1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 PetscScalar diag; 6 } Mat_ConstantDiagonal; 7 8 static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 9 { 10 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)A->data; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 if (ncols) *ncols = 1; 15 if (cols) { 16 ierr = PetscMalloc1(1,cols);CHKERRQ(ierr); 17 (*cols)[0] = row; 18 } 19 if (vals) { 20 ierr = PetscMalloc1(1,vals);CHKERRQ(ierr); 21 (*vals)[0] = ctx->diag; 22 } 23 PetscFunctionReturn(0); 24 } 25 26 static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[]) 27 { 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 if (ncols) *ncols = 0; 32 if (cols) { 33 ierr = PetscFree(*cols);CHKERRQ(ierr); 34 } 35 if (vals) { 36 ierr = PetscFree(*vals);CHKERRQ(ierr); 37 } 38 PetscFunctionReturn(0); 39 } 40 41 static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y) 42 { 43 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)A->data; 44 PetscErrorCode ierr; 45 46 PetscFunctionBegin; 47 ierr = VecAXPBY(y,ctx->diag,0.0,x);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3) 52 { 53 PetscErrorCode ierr; 54 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)mat->data; 55 56 PetscFunctionBegin; 57 if (v2 == v3) { 58 ierr = VecAXPBY(v3,ctx->diag,1.0,v1);CHKERRQ(ierr); 59 } else { 60 ierr = VecAXPBYPCZ(v3,ctx->diag,1.0,0.0,v1,v2);CHKERRQ(ierr); 61 } 62 PetscFunctionReturn(0); 63 } 64 65 static PetscErrorCode MatMultTransposeAdd_ConstantDiagonal(Mat mat,Vec v1,Vec v2,Vec v3) 66 { 67 PetscErrorCode ierr; 68 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)mat->data; 69 70 PetscFunctionBegin; 71 if (v2 == v3) { 72 ierr = VecAXPBY(v3,ctx->diag,1.0,v1);CHKERRQ(ierr); 73 } else { 74 ierr = VecAXPBYPCZ(v3,ctx->diag,1.0,0.0,v1,v2);CHKERRQ(ierr); 75 } 76 PetscFunctionReturn(0); 77 } 78 79 static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B) 80 { 81 PetscErrorCode ierr; 82 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data; 83 84 PetscFunctionBegin; 85 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 86 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 87 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 88 ierr = MatSetType(*B,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 89 ierr = PetscLayoutReference(A->rmap,&(*B)->rmap);CHKERRQ(ierr); 90 ierr = PetscLayoutReference(A->cmap,&(*B)->cmap);CHKERRQ(ierr); 91 if (op == MAT_COPY_VALUES) { 92 Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal*)(*B)->data; 93 bctx->diag = actx->diag; 94 } 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat,PetscBool *missing,PetscInt *dd) 99 { 100 PetscFunctionBegin; 101 *missing = PETSC_FALSE; 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat) 106 { 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = PetscFree(mat->data);CHKERRQ(ierr); 111 PetscFunctionReturn(0); 112 } 113 114 static PetscErrorCode MatView_ConstantDiagonal(Mat J,PetscViewer viewer) 115 { 116 PetscErrorCode ierr; 117 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 118 PetscBool iascii; 119 120 PetscFunctionBegin; 121 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 122 if (iascii) { 123 PetscViewerFormat format; 124 125 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 126 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(0); 127 #if !defined(PETSC_USE_COMPLEX) 128 ierr = PetscViewerASCIIPrintf(viewer,"Diagonal value: %g\n",ctx->diag);CHKERRQ(ierr); 129 #else 130 ierr = PetscViewerASCIIPrintf(viewer,"Diagonal value: %g + i %g\n",PetscRealPart(ctx->diag),PetscImaginaryPart(ctx->diag));CHKERRQ(ierr); 131 #endif 132 } 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J,MatAssemblyType mt) 137 { 138 PetscFunctionBegin; 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode MatMult_ConstantDiagonal(Mat J,Vec x,Vec y) 143 { 144 PetscErrorCode ierr; 145 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 146 147 PetscFunctionBegin; 148 ierr = VecAXPBY(y,ctx->diag,0.0,x);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J,Vec x) 153 { 154 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)J->data; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = VecSet(x,ctx->diag);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode MatShift_ConstantDiagonal(Mat Y,PetscScalar a) 163 { 164 PetscErrorCode ierr; 165 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 166 167 PetscFunctionBegin; 168 if (a != 0.) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 169 ctx->diag += a; 170 PetscFunctionReturn(0); 171 } 172 173 static PetscErrorCode MatScale_ConstantDiagonal(Mat Y,PetscScalar a) 174 { 175 PetscErrorCode ierr; 176 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 177 178 PetscFunctionBegin; 179 if (a != 1.) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 180 ctx->diag *= a; 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y) 185 { 186 PetscErrorCode ierr; 187 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)Y->data; 188 189 PetscFunctionBegin; 190 if (ctx->diag != 0.0) {ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr);} 191 ctx->diag = 0.0; 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode MatSOR_ConstantDiagonal(Mat matin,Vec x,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec y) 196 { 197 PetscErrorCode ierr; 198 Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal*)matin->data; 199 200 PetscFunctionBegin; 201 if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 202 else matin->factorerrortype = MAT_FACTOR_NOERROR; 203 ierr = VecAXPBY(y,1.0/ctx->diag,0.0,x);CHKERRQ(ierr); 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A,MatInfoType flag,MatInfo *info) 208 { 209 PetscFunctionBegin; 210 info->block_size = 1.0; 211 info->nz_allocated = 1.0; 212 info->nz_used = 1.0; 213 info->nz_unneeded = 0.0; 214 info->assemblies = A->num_ass; 215 info->mallocs = 0.0; 216 info->memory = ((PetscObject)A)->mem; 217 if (A->factortype) { 218 info->fill_ratio_given = 1.0; 219 info->fill_ratio_needed = 1.0; 220 info->factor_mallocs = 0.0; 221 } else { 222 info->fill_ratio_given = 0; 223 info->fill_ratio_needed = 0; 224 info->factor_mallocs = 0; 225 } 226 PetscFunctionReturn(0); 227 } 228 229 /*@ 230 MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal 231 232 Collective 233 234 Input Parameters: 235 + comm - MPI communicator 236 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 237 This value should be the same as the local size used in creating the 238 y vector for the matrix-vector product y = Ax. 239 . n - This value should be the same as the local size used in creating the 240 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 241 calculated if N is given) For square matrices n is almost always m. 242 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 243 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 244 - diag - the diagonal value 245 246 Output Parameter: 247 . J - the diagonal matrix 248 249 Level: advanced 250 251 Notes: 252 Only supports square matrices with the same number of local rows and columns 253 254 .seealso: MatDestroy(), MATCONSTANTDIAGONAL, MatScale(), MatShift(), MatMult(), MatGetDiagonal(), MatGetFactor(), MatSolve() 255 256 @*/ 257 PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar diag,Mat *J) 258 { 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 ierr = MatCreate(comm,J);CHKERRQ(ierr); 263 ierr = MatSetSizes(*J,m,n,M,N);CHKERRQ(ierr); 264 ierr = MatSetType(*J,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 265 ierr = MatShift(*J,diag);CHKERRQ(ierr); 266 ierr = MatSetUp(*J);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A) 271 { 272 PetscErrorCode ierr; 273 Mat_ConstantDiagonal *ctx; 274 275 PetscFunctionBegin; 276 ierr = PetscNew(&ctx);CHKERRQ(ierr); 277 ctx->diag = 0.0; 278 A->data = (void*)ctx; 279 280 A->assembled = PETSC_TRUE; 281 A->preallocated = PETSC_TRUE; 282 283 A->ops->mult = MatMult_ConstantDiagonal; 284 A->ops->multadd = MatMultAdd_ConstantDiagonal; 285 A->ops->multtranspose = MatMultTranspose_ConstantDiagonal; 286 A->ops->multtransposeadd = MatMultTransposeAdd_ConstantDiagonal; 287 A->ops->duplicate = MatDuplicate_ConstantDiagonal; 288 A->ops->missingdiagonal = MatMissingDiagonal_ConstantDiagonal; 289 A->ops->getrow = MatGetRow_ConstantDiagonal; 290 A->ops->restorerow = MatRestoreRow_ConstantDiagonal; 291 A->ops->sor = MatSOR_ConstantDiagonal; 292 A->ops->shift = MatShift_ConstantDiagonal; 293 A->ops->scale = MatScale_ConstantDiagonal; 294 A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal; 295 A->ops->view = MatView_ConstantDiagonal; 296 A->ops->zeroentries = MatZeroEntries_ConstantDiagonal; 297 A->ops->assemblyend = MatAssemblyEnd_ConstantDiagonal; 298 A->ops->destroy = MatDestroy_ConstantDiagonal; 299 A->ops->getinfo = MatGetInfo_ConstantDiagonal; 300 301 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCONSTANTDIAGONAL);CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact,Mat A,const MatFactorInfo *info) 306 { 307 PetscErrorCode ierr; 308 Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal*)A->data,*fctx = (Mat_ConstantDiagonal*)fact->data; 309 310 PetscFunctionBegin; 311 if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 312 else fact->factorerrortype = MAT_FACTOR_NOERROR; 313 fctx->diag = 1.0/actx->diag; 314 ierr = PetscObjectStateIncrease((PetscObject)fact);CHKERRQ(ierr); 315 fact->ops->solve = MatMult_ConstantDiagonal; 316 PetscFunctionReturn(0); 317 } 318 319 static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 320 { 321 PetscFunctionBegin; 322 fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal; 323 PetscFunctionReturn(0); 324 } 325 326 static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact,Mat A,IS isrow,const MatFactorInfo *info) 327 { 328 PetscFunctionBegin; 329 fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal; PetscFunctionReturn(0); 330 PetscFunctionReturn(0); 331 } 332 333 PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A,MatFactorType ftype,Mat *B) 334 { 335 PetscInt n = A->rmap->n, N = A->rmap->N; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 ierr = MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A),n,n,N,N,0,B);CHKERRQ(ierr); 340 341 (*B)->factortype = ftype; 342 (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 343 (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal; 344 (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 345 (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal; 346 347 (*B)->ops->shift = NULL; 348 (*B)->ops->scale = NULL; 349 (*B)->ops->mult = NULL; 350 (*B)->ops->sor = NULL; 351 (*B)->ops->zeroentries = NULL; 352 353 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 354 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 355 PetscFunctionReturn(0); 356 } 357