Actual source code: cdiagonal.c

  1: #include <petsc/private/matimpl.h>

  3: typedef struct {
  4:   PetscScalar diag;
  5: } Mat_ConstantDiagonal;

  7: static PetscErrorCode MatAXPY_ConstantDiagonal(Mat Y, PetscScalar a, Mat X, MatStructure str)
  8: {
  9:   Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data;
 10:   Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data;

 12:   PetscFunctionBegin;
 13:   yctx->diag += a * xctx->diag;
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

 17: static PetscErrorCode MatEqual_ConstantDiagonal(Mat Y, Mat X, PetscBool *equal)
 18: {
 19:   Mat_ConstantDiagonal *yctx = (Mat_ConstantDiagonal *)Y->data;
 20:   Mat_ConstantDiagonal *xctx = (Mat_ConstantDiagonal *)X->data;

 22:   PetscFunctionBegin;
 23:   *equal = (yctx->diag == xctx->diag) ? PETSC_TRUE : PETSC_FALSE;
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: static PetscErrorCode MatGetRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
 28: {
 29:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;

 31:   PetscFunctionBegin;
 32:   if (ncols) *ncols = 1;
 33:   if (cols) {
 34:     PetscCall(PetscMalloc1(1, cols));
 35:     (*cols)[0] = row;
 36:   }
 37:   if (vals) {
 38:     PetscCall(PetscMalloc1(1, vals));
 39:     (*vals)[0] = ctx->diag;
 40:   }
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
 45: {
 46:   PetscFunctionBegin;
 47:   if (cols) PetscCall(PetscFree(*cols));
 48:   if (vals) PetscCall(PetscFree(*vals));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
 53: {
 54:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;

 56:   PetscFunctionBegin;
 57:   if (v2 == v3) {
 58:     PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1));
 59:   } else {
 60:     PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2));
 61:   }
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
 66: {
 67:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;

 69:   PetscFunctionBegin;
 70:   if (v2 == v3) {
 71:     PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1));
 72:   } else {
 73:     PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2));
 74:   }
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
 79: {
 80:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;

 82:   PetscFunctionBegin;
 83:   if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag);
 84:   else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
 85:   PetscFunctionReturn(PETSC_SUCCESS);
 86: }

 88: static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])

 90: {
 91:   Mat B;

 93:   PetscFunctionBegin;
 94:   PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
 95:   PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat));
 96:   PetscCall(MatDestroy(&B));
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
101: {
102:   Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;

104:   PetscFunctionBegin;
105:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
106:   PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
107:   PetscCall(MatSetBlockSizesFromMats(*B, A, A));
108:   PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL));
109:   PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
110:   PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
111:   if (op == MAT_COPY_VALUES) {
112:     Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
113:     bctx->diag                 = actx->diag;
114:   }
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd)
119: {
120:   PetscFunctionBegin;
121:   *missing = PETSC_FALSE;
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
126: {
127:   PetscFunctionBegin;
128:   PetscCall(PetscFree(mat->data));
129:   mat->structural_symmetry_eternal = PETSC_FALSE;
130:   mat->symmetry_eternal            = PETSC_FALSE;
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
135: {
136:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
137:   PetscBool             iascii;

139:   PetscFunctionBegin;
140:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
141:   if (iascii) {
142:     PetscViewerFormat format;

144:     PetscCall(PetscViewerGetFormat(viewer, &format));
145:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
146: #if defined(PETSC_USE_COMPLEX)
147:     PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
148: #else
149:     PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag)));
150: #endif
151:   }
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt)
156: {
157:   PetscFunctionBegin;
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
162: {
163:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

165:   PetscFunctionBegin;
166:   PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
167:   PetscFunctionReturn(PETSC_SUCCESS);
168: }

170: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
171: {
172:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

174:   PetscFunctionBegin;
175:   PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
180: {
181:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

183:   PetscFunctionBegin;
184:   PetscCall(VecSet(x, ctx->diag));
185:   PetscFunctionReturn(PETSC_SUCCESS);
186: }

188: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
189: {
190:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

192:   PetscFunctionBegin;
193:   ctx->diag += a;
194:   PetscFunctionReturn(PETSC_SUCCESS);
195: }

197: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
198: {
199:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

201:   PetscFunctionBegin;
202:   ctx->diag *= a;
203:   PetscFunctionReturn(PETSC_SUCCESS);
204: }

206: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
207: {
208:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

210:   PetscFunctionBegin;
211:   ctx->diag = 0.0;
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }

215: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
216: {
217:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;

219:   PetscFunctionBegin;
220:   if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
221:   else matin->factorerrortype = MAT_FACTOR_NOERROR;
222:   PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
227: {
228:   PetscFunctionBegin;
229:   PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
230:   PetscFunctionReturn(PETSC_SUCCESS);
231: }

233: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
234: {
235:   PetscFunctionBegin;
236:   info->block_size   = 1.0;
237:   info->nz_allocated = 1.0;
238:   info->nz_used      = 1.0;
239:   info->nz_unneeded  = 0.0;
240:   info->assemblies   = A->num_ass;
241:   info->mallocs      = 0.0;
242:   info->memory       = 0; /* REVIEW ME */
243:   if (A->factortype) {
244:     info->fill_ratio_given  = 1.0;
245:     info->fill_ratio_needed = 1.0;
246:     info->factor_mallocs    = 0.0;
247:   } else {
248:     info->fill_ratio_given  = 0;
249:     info->fill_ratio_needed = 0;
250:     info->factor_mallocs    = 0;
251:   }
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: /*@
256:   MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal

258:   Collective

260:   Input Parameters:
261: + comm - MPI communicator
262: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
263:            This value should be the same as the local size used in creating the
264:            y vector for the matrix-vector product y = Ax.
265: . n    - This value should be the same as the local size used in creating the
266:        x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
267:        calculated if `N` is given) For square matrices n is almost always `m`.
268: . M    - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
269: . N    - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
270: - diag - the diagonal value

272:   Output Parameter:
273: . J - the diagonal matrix

275:   Level: advanced

277:   Notes:
278:   Only supports square matrices with the same number of local rows and columns

280: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
281: @*/
282: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
283: {
284:   PetscFunctionBegin;
285:   PetscCall(MatCreate(comm, J));
286:   PetscCall(MatSetSizes(*J, m, n, M, N));
287:   PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
288:   PetscCall(MatShift(*J, diag));
289:   PetscCall(MatSetUp(*J));
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
294: {
295:   Mat_ConstantDiagonal *ctx;

297:   PetscFunctionBegin;
298:   PetscCall(PetscNew(&ctx));
299:   ctx->diag = 0.0;
300:   A->data   = (void *)ctx;

302:   A->assembled                   = PETSC_TRUE;
303:   A->preallocated                = PETSC_TRUE;
304:   A->structurally_symmetric      = PETSC_BOOL3_TRUE;
305:   A->structural_symmetry_eternal = PETSC_TRUE;
306:   A->symmetric                   = PETSC_BOOL3_TRUE;
307:   if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
308:   A->symmetry_eternal = PETSC_TRUE;

310:   A->ops->mult                      = MatMult_ConstantDiagonal;
311:   A->ops->multadd                   = MatMultAdd_ConstantDiagonal;
312:   A->ops->multtranspose             = MatMult_ConstantDiagonal;
313:   A->ops->multtransposeadd          = MatMultAdd_ConstantDiagonal;
314:   A->ops->multhermitiantranspose    = MatMultHermitianTranspose_ConstantDiagonal;
315:   A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
316:   A->ops->solve                     = MatSolve_ConstantDiagonal;
317:   A->ops->solvetranspose            = MatSolve_ConstantDiagonal;
318:   A->ops->norm                      = MatNorm_ConstantDiagonal;
319:   A->ops->createsubmatrices         = MatCreateSubMatrices_ConstantDiagonal;
320:   A->ops->duplicate                 = MatDuplicate_ConstantDiagonal;
321:   A->ops->missingdiagonal           = MatMissingDiagonal_ConstantDiagonal;
322:   A->ops->getrow                    = MatGetRow_ConstantDiagonal;
323:   A->ops->restorerow                = MatRestoreRow_ConstantDiagonal;
324:   A->ops->sor                       = MatSOR_ConstantDiagonal;
325:   A->ops->shift                     = MatShift_ConstantDiagonal;
326:   A->ops->scale                     = MatScale_ConstantDiagonal;
327:   A->ops->getdiagonal               = MatGetDiagonal_ConstantDiagonal;
328:   A->ops->view                      = MatView_ConstantDiagonal;
329:   A->ops->zeroentries               = MatZeroEntries_ConstantDiagonal;
330:   A->ops->assemblyend               = MatAssemblyEnd_ConstantDiagonal;
331:   A->ops->destroy                   = MatDestroy_ConstantDiagonal;
332:   A->ops->getinfo                   = MatGetInfo_ConstantDiagonal;
333:   A->ops->equal                     = MatEqual_ConstantDiagonal;
334:   A->ops->axpy                      = MatAXPY_ConstantDiagonal;

336:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
341: {
342:   Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;

344:   PetscFunctionBegin;
345:   if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
346:   else fact->factorerrortype = MAT_FACTOR_NOERROR;
347:   fctx->diag       = 1.0 / actx->diag;
348:   fact->ops->solve = MatMult_ConstantDiagonal;
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

352: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
353: {
354:   PetscFunctionBegin;
355:   fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
360: {
361:   PetscFunctionBegin;
362:   fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
363:   PetscFunctionReturn(PETSC_SUCCESS);
364: }

366: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
367: {
368:   PetscInt n = A->rmap->n, N = A->rmap->N;

370:   PetscFunctionBegin;
371:   PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));

373:   (*B)->factortype                  = ftype;
374:   (*B)->ops->ilufactorsymbolic      = MatFactorSymbolic_LU_ConstantDiagonal;
375:   (*B)->ops->lufactorsymbolic       = MatFactorSymbolic_LU_ConstantDiagonal;
376:   (*B)->ops->iccfactorsymbolic      = MatFactorSymbolic_Cholesky_ConstantDiagonal;
377:   (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;

379:   (*B)->ops->shift       = NULL;
380:   (*B)->ops->scale       = NULL;
381:   (*B)->ops->mult        = NULL;
382:   (*B)->ops->sor         = NULL;
383:   (*B)->ops->zeroentries = NULL;

385:   PetscCall(PetscFree((*B)->solvertype));
386:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }