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:   yctx->diag += a * xctx->diag;
 14:   return 0;
 15: }

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

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

 33: static PetscErrorCode MatRestoreRow_ConstantDiagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt *cols[], PetscScalar *vals[])
 34: {
 35:   if (ncols) *ncols = 0;
 36:   if (cols) PetscFree(*cols);
 37:   if (vals) PetscFree(*vals);
 38:   return 0;
 39: }

 41: static PetscErrorCode MatMultTranspose_ConstantDiagonal(Mat A, Vec x, Vec y)
 42: {
 43:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;

 45:   VecAXPBY(y, ctx->diag, 0.0, x);
 46:   return 0;
 47: }

 49: static PetscErrorCode MatMultAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
 50: {
 51:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;

 53:   if (v2 == v3) {
 54:     VecAXPBY(v3, ctx->diag, 1.0, v1);
 55:   } else {
 56:     VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2);
 57:   }
 58:   return 0;
 59: }

 61: static PetscErrorCode MatMultTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
 62: {
 63:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;

 65:   if (v2 == v3) {
 66:     VecAXPBY(v3, ctx->diag, 1.0, v1);
 67:   } else {
 68:     VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2);
 69:   }
 70:   return 0;
 71: }

 73: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
 74: {
 75:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;

 77:   if (type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY) *nrm = PetscAbsScalar(ctx->diag);
 78:   else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
 79:   return 0;
 80: }

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

 84: {
 85:   Mat B;

 87:   MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B);
 88:   MatCreateSubMatrices(B, n, irow, icol, scall, submat);
 89:   MatDestroy(&B);
 90:   return 0;
 91: }

 93: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
 94: {
 95:   Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;

 97:   MatCreate(PetscObjectComm((PetscObject)A), B);
 98:   MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
 99:   MatSetBlockSizesFromMats(*B, A, A);
100:   MatSetType(*B, MATCONSTANTDIAGONAL);
101:   PetscLayoutReference(A->rmap, &(*B)->rmap);
102:   PetscLayoutReference(A->cmap, &(*B)->cmap);
103:   if (op == MAT_COPY_VALUES) {
104:     Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
105:     bctx->diag                 = actx->diag;
106:   }
107:   return 0;
108: }

110: static PetscErrorCode MatMissingDiagonal_ConstantDiagonal(Mat mat, PetscBool *missing, PetscInt *dd)
111: {
112:   *missing = PETSC_FALSE;
113:   return 0;
114: }

116: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
117: {
118:   PetscFree(mat->data);
119:   return 0;
120: }

122: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
123: {
124:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
125:   PetscBool             iascii;

127:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
128:   if (iascii) {
129:     PetscViewerFormat format;

131:     PetscViewerGetFormat(viewer, &format);
132:     if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) return 0;
133: #if defined(PETSC_USE_COMPLEX)
134:     PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag));
135: #else
136:     PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)(ctx->diag));
137: #endif
138:   }
139:   return 0;
140: }

142: static PetscErrorCode MatAssemblyEnd_ConstantDiagonal(Mat J, MatAssemblyType mt)
143: {
144:   return 0;
145: }

147: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
148: {
149:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

151:   VecAXPBY(y, ctx->diag, 0.0, x);
152:   return 0;
153: }

155: PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
156: {
157:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;

159:   VecSet(x, ctx->diag);
160:   return 0;
161: }

163: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
164: {
165:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

167:   ctx->diag += a;
168:   return 0;
169: }

171: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
172: {
173:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

175:   ctx->diag *= a;
176:   return 0;
177: }

179: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
180: {
181:   Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;

183:   ctx->diag = 0.0;
184:   return 0;
185: }

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

191:   if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
192:   else matin->factorerrortype = MAT_FACTOR_NOERROR;
193:   VecAXPBY(y, 1.0 / ctx->diag, 0.0, x);
194:   return 0;
195: }

197: PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
198: {
199:   info->block_size   = 1.0;
200:   info->nz_allocated = 1.0;
201:   info->nz_used      = 1.0;
202:   info->nz_unneeded  = 0.0;
203:   info->assemblies   = A->num_ass;
204:   info->mallocs      = 0.0;
205:   info->memory       = 0; /* REVIEW ME */
206:   if (A->factortype) {
207:     info->fill_ratio_given  = 1.0;
208:     info->fill_ratio_needed = 1.0;
209:     info->factor_mallocs    = 0.0;
210:   } else {
211:     info->fill_ratio_given  = 0;
212:     info->fill_ratio_needed = 0;
213:     info->factor_mallocs    = 0;
214:   }
215:   return 0;
216: }

218: /*@
219:    MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal

221:    Collective

223:    Input Parameters:
224: +  comm - MPI communicator
225: .  m - number of local rows (or `PETSC_DECIDE` to have calculated if M is given)
226:            This value should be the same as the local size used in creating the
227:            y vector for the matrix-vector product y = Ax.
228: .  n - This value should be the same as the local size used in creating the
229:        x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have
230:        calculated if N is given) For square matrices n is almost always m.
231: .  M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
232: .  N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
233: -  diag - the diagonal value

235:    Output Parameter:
236: .  J - the diagonal matrix

238:    Level: advanced

240:    Notes:
241:     Only supports square matrices with the same number of local rows and columns

243: .seealso: `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
244: @*/
245: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
246: {
247:   MatCreate(comm, J);
248:   MatSetSizes(*J, m, n, M, N);
249:   MatSetType(*J, MATCONSTANTDIAGONAL);
250:   MatShift(*J, diag);
251:   MatSetUp(*J);
252:   return 0;
253: }

255: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
256: {
257:   Mat_ConstantDiagonal *ctx;

259:   PetscNew(&ctx);
260:   ctx->diag = 0.0;
261:   A->data   = (void *)ctx;

263:   A->assembled    = PETSC_TRUE;
264:   A->preallocated = PETSC_TRUE;

266:   A->ops->mult              = MatMult_ConstantDiagonal;
267:   A->ops->multadd           = MatMultAdd_ConstantDiagonal;
268:   A->ops->multtranspose     = MatMultTranspose_ConstantDiagonal;
269:   A->ops->multtransposeadd  = MatMultTransposeAdd_ConstantDiagonal;
270:   A->ops->norm              = MatNorm_ConstantDiagonal;
271:   A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
272:   A->ops->duplicate         = MatDuplicate_ConstantDiagonal;
273:   A->ops->missingdiagonal   = MatMissingDiagonal_ConstantDiagonal;
274:   A->ops->getrow            = MatGetRow_ConstantDiagonal;
275:   A->ops->restorerow        = MatRestoreRow_ConstantDiagonal;
276:   A->ops->sor               = MatSOR_ConstantDiagonal;
277:   A->ops->shift             = MatShift_ConstantDiagonal;
278:   A->ops->scale             = MatScale_ConstantDiagonal;
279:   A->ops->getdiagonal       = MatGetDiagonal_ConstantDiagonal;
280:   A->ops->view              = MatView_ConstantDiagonal;
281:   A->ops->zeroentries       = MatZeroEntries_ConstantDiagonal;
282:   A->ops->assemblyend       = MatAssemblyEnd_ConstantDiagonal;
283:   A->ops->destroy           = MatDestroy_ConstantDiagonal;
284:   A->ops->getinfo           = MatGetInfo_ConstantDiagonal;
285:   A->ops->axpy              = MatAXPY_ConstantDiagonal;

287:   PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL);
288:   return 0;
289: }

291: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
292: {
293:   Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;

295:   if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
296:   else fact->factorerrortype = MAT_FACTOR_NOERROR;
297:   fctx->diag       = 1.0 / actx->diag;
298:   fact->ops->solve = MatMult_ConstantDiagonal;
299:   return 0;
300: }

302: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
303: {
304:   fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
305:   return 0;
306: }

308: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
309: {
310:   fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
311:   return 0;
312: }

314: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
315: {
316:   PetscInt n = A->rmap->n, N = A->rmap->N;

318:   MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B);

320:   (*B)->factortype                  = ftype;
321:   (*B)->ops->ilufactorsymbolic      = MatFactorSymbolic_LU_ConstantDiagonal;
322:   (*B)->ops->lufactorsymbolic       = MatFactorSymbolic_LU_ConstantDiagonal;
323:   (*B)->ops->iccfactorsymbolic      = MatFactorSymbolic_Cholesky_ConstantDiagonal;
324:   (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;

326:   (*B)->ops->shift       = NULL;
327:   (*B)->ops->scale       = NULL;
328:   (*B)->ops->mult        = NULL;
329:   (*B)->ops->sor         = NULL;
330:   (*B)->ops->zeroentries = NULL;

332:   PetscFree((*B)->solvertype);
333:   PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype);
334:   return 0;
335: }