Actual source code: cdiagonal.c


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

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

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

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

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

 22:   PetscFunctionBegin;
 23:   if (ncols) *ncols = 1;
 24:   if (cols) {
 25:     PetscCall(PetscMalloc1(1, cols));
 26:     (*cols)[0] = row;
 27:   }
 28:   if (vals) {
 29:     PetscCall(PetscMalloc1(1, vals));
 30:     (*vals)[0] = ctx->diag;
 31:   }
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
 36: {
 37:   PetscFunctionBegin;
 38:   if (ncols) *ncols = 0;
 39:   if (cols) PetscCall(PetscFree(*cols));
 40:   if (vals) PetscCall(PetscFree(*vals));
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y)
 45: {
 46:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;

 48:   PetscFunctionBegin;
 49:   PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

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

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

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

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

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

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

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

 91: {
 92:   Mat B;

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

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

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

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

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

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

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

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

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

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

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

169: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
170: {
171:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

173:   PetscFunctionBegin;
174:   PetscCall(VecSet(x, ctx->diag));
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
179: {
180:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

182:   PetscFunctionBegin;
183:   ctx->diag += a;
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

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

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

196: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
197: {
198:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

200:   PetscFunctionBegin;
201:   ctx->diag = 0.0;
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
206: {
207:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;

209:   PetscFunctionBegin;
210:   if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
211:   else matin->factorerrortype = MAT_FACTOR_NOERROR;
212:   PetscCall(VecAXPBY(y, 1.0 / ctx->diag, 0.0, x));
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
217: {
218:   PetscFunctionBegin;
219:   info->block_size   = 1.0;
220:   info->nz_allocated = 1.0;
221:   info->nz_used      = 1.0;
222:   info->nz_unneeded  = 0.0;
223:   info->assemblies   = A->num_ass;
224:   info->mallocs      = 0.0;
225:   info->memory       = 0; /* REVIEW ME */
226:   if (A->factortype) {
227:     info->fill_ratio_given  = 1.0;
228:     info->fill_ratio_needed = 1.0;
229:     info->factor_mallocs    = 0.0;
230:   } else {
231:     info->fill_ratio_given  = 0;
232:     info->fill_ratio_needed = 0;
233:     info->factor_mallocs    = 0;
234:   }
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: /*@
239:    MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal

241:    Collective

243:    Input Parameters:
244: +  comm - MPI communicator
245: .  m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
246:            This value should be the same as the local size used in creating the
247:            y vector for the matrix-vector product y = Ax.
248: .  n - This value should be the same as the local size used in creating the
249:        x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
250:        calculated if `N` is given) For square matrices n is almost always `m`.
251: .  M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
252: .  N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
253: -  diag - the diagonal value

255:    Output Parameter:
256: .  J - the diagonal matrix

258:    Level: advanced

260:    Notes:
261:     Only supports square matrices with the same number of local rows and columns

263: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
264: @*/
265: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
266: {
267:   PetscFunctionBegin;
268:   PetscCall(MatCreate(comm, J));
269:   PetscCall(MatSetSizes(*J, m, n, M, N));
270:   PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
271:   PetscCall(MatShift(*J, diag));
272:   PetscCall(MatSetUp(*J));
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
277: {
278:   Mat_ConstantDiagonal *ctx;

280:   PetscFunctionBegin;
281:   PetscCall(PetscNew(&ctx));
282:   ctx->diag = 0.0;
283:   A->data   = (void *)ctx;

285:   A->assembled    = PETSC_TRUE;
286:   A->preallocated = PETSC_TRUE;

288:   A->ops->mult              = MatMult_ConstantDiagonal;
289:   A->ops->multadd           = MatMultAdd_ConstantDiagonal;
290:   A->ops->multtranspose     = MatMultTranspose_ConstantDiagonal;
291:   A->ops->multtransposeadd  = MatMultTransposeAdd_ConstantDiagonal;
292:   A->ops->norm              = MatNorm_ConstantDiagonal;
293:   A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
294:   A->ops->duplicate         = MatDuplicate_ConstantDiagonal;
295:   A->ops->missingdiagonal   = MatMissingDiagonal_ConstantDiagonal;
296:   A->ops->getrow            = MatGetRow_ConstantDiagonal;
297:   A->ops->restorerow        = MatRestoreRow_ConstantDiagonal;
298:   A->ops->sor               = MatSOR_ConstantDiagonal;
299:   A->ops->shift             = MatShift_ConstantDiagonal;
300:   A->ops->scale             = MatScale_ConstantDiagonal;
301:   A->ops->getdiagonal       = MatGetDiagonal_ConstantDiagonal;
302:   A->ops->view              = MatView_ConstantDiagonal;
303:   A->ops->zeroentries       = MatZeroEntries_ConstantDiagonal;
304:   A->ops->assemblyend       = MatAssemblyEnd_ConstantDiagonal;
305:   A->ops->destroy           = MatDestroy_ConstantDiagonal;
306:   A->ops->getinfo           = MatGetInfo_ConstantDiagonal;
307:   A->ops->axpy              = MatAXPY_ConstantDiagonal;

309:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
314: {
315:   Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;

317:   PetscFunctionBegin;
318:   if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
319:   else fact->factorerrortype = MAT_FACTOR_NOERROR;
320:   fctx->diag       = 1.0 / actx->diag;
321:   fact->ops->solve = MatMult_ConstantDiagonal;
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

325: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
326: {
327:   PetscFunctionBegin;
328:   fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
329:   PetscFunctionReturn(PETSC_SUCCESS);
330: }

332: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
333: {
334:   PetscFunctionBegin;
335:   fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
340: {
341:   PetscInt n = A->rmap->n, N = A->rmap->N;

343:   PetscFunctionBegin;
344:   PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));

346:   (*B)->factortype                  = ftype;
347:   (*B)->ops->ilufactorsymbolic      = MatFactorSymbolic_LU_ConstantDiagonal;
348:   (*B)->ops->lufactorsymbolic       = MatFactorSymbolic_LU_ConstantDiagonal;
349:   (*B)->ops->iccfactorsymbolic      = MatFactorSymbolic_Cholesky_ConstantDiagonal;
350:   (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;

352:   (*B)->ops->shift       = NULL;
353:   (*B)->ops->scale       = NULL;
354:   (*B)->ops->mult        = NULL;
355:   (*B)->ops->sor         = NULL;
356:   (*B)->ops->zeroentries = NULL;

358:   PetscCall(PetscFree((*B)->solvertype));
359:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
360:   PetscFunctionReturn(PETSC_SUCCESS);
361: }