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) PetscCall(VecAXPBY(v3, ctx->diag, 1.0, v1));
58: else PetscCall(VecAXPBYPCZ(v3, ctx->diag, 1.0, 0.0, v1, v2));
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode MatMultHermitianTransposeAdd_ConstantDiagonal(Mat mat, Vec v1, Vec v2, Vec v3)
63: {
64: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
66: PetscFunctionBegin;
67: if (v2 == v3) PetscCall(VecAXPBY(v3, PetscConj(ctx->diag), 1.0, v1));
68: else PetscCall(VecAXPBYPCZ(v3, PetscConj(ctx->diag), 1.0, 0.0, v1, v2));
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: static PetscErrorCode MatNorm_ConstantDiagonal(Mat A, NormType type, PetscReal *nrm)
73: {
74: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
76: PetscFunctionBegin;
77: PetscCheck(type == NORM_FROBENIUS || type == NORM_2 || type == NORM_1 || type == NORM_INFINITY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm");
78: *nrm = PetscAbsScalar(ctx->diag);
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode MatCreateSubMatrices_ConstantDiagonal(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
83: {
84: Mat B;
86: PetscFunctionBegin;
87: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B));
88: PetscCall(MatCreateSubMatrices(B, n, irow, icol, scall, submat));
89: PetscCall(MatDestroy(&B));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode MatDuplicate_ConstantDiagonal(Mat A, MatDuplicateOption op, Mat *B)
94: {
95: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data;
97: PetscFunctionBegin;
98: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
99: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
100: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
101: PetscCall(MatSetType(*B, MATCONSTANTDIAGONAL));
102: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
103: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
104: if (op == MAT_COPY_VALUES) {
105: Mat_ConstantDiagonal *bctx = (Mat_ConstantDiagonal *)(*B)->data;
106: bctx->diag = actx->diag;
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode MatDestroy_ConstantDiagonal(Mat mat)
112: {
113: PetscFunctionBegin;
114: PetscCall(PetscFree(mat->data));
115: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConstantDiagonalGetConstant_C", NULL));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
120: {
121: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
122: PetscBool isascii;
124: PetscFunctionBegin;
125: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
126: if (isascii) {
127: PetscViewerFormat format;
129: PetscCall(PetscViewerGetFormat(viewer, &format));
130: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
131: if (PetscImaginaryPart(ctx->diag) == 0) {
132: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag)));
133: } else {
134: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
135: }
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
141: {
142: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
144: PetscFunctionBegin;
145: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
150: {
151: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
153: PetscFunctionBegin;
154: PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
159: {
160: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
162: PetscFunctionBegin;
163: PetscCall(VecSet(x, ctx->diag));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
168: {
169: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
171: PetscFunctionBegin;
172: ctx->diag += a;
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
177: {
178: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
180: PetscFunctionBegin;
181: ctx->diag *= a;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
186: {
187: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
189: PetscFunctionBegin;
190: ctx->diag = 0.0;
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y)
195: {
196: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
198: PetscFunctionBegin;
199: ctx->diag = PetscConj(ctx->diag);
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout)
204: {
205: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
207: PetscFunctionBegin;
208: if (reuse == MAT_INPLACE_MATRIX) {
209: PetscLayout tmplayout = A->rmap;
211: A->rmap = A->cmap;
212: A->cmap = tmplayout;
213: } else {
214: if (reuse == MAT_INITIAL_MATRIX) {
215: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout));
216: } else {
217: PetscCall(MatZeroEntries(*matout));
218: PetscCall(MatShift(*matout, ctx->diag));
219: }
220: }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand)
225: {
226: PetscMPIInt rank;
227: MPI_Comm comm;
228: PetscScalar v = 0.0;
229: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
231: PetscFunctionBegin;
232: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
233: PetscCallMPI(MPI_Comm_rank(comm, &rank));
234: if (!rank) PetscCall(PetscRandomGetValue(rand, &v));
235: PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm));
236: ctx->diag = v;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
241: {
242: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
244: PetscFunctionBegin;
245: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
246: else matin->factorerrortype = MAT_FACTOR_NOERROR;
247: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
252: {
253: PetscFunctionBegin;
254: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
259: {
260: PetscFunctionBegin;
261: info->block_size = 1.0;
262: info->nz_allocated = 1.0;
263: info->nz_used = 1.0;
264: info->nz_unneeded = 0.0;
265: info->assemblies = A->num_ass;
266: info->mallocs = 0.0;
267: info->memory = 0; /* REVIEW ME */
268: if (A->factortype) {
269: info->fill_ratio_given = 1.0;
270: info->fill_ratio_needed = 1.0;
271: info->factor_mallocs = 0.0;
272: } else {
273: info->fill_ratio_given = 0;
274: info->fill_ratio_needed = 0;
275: info->factor_mallocs = 0;
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@
281: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
283: Collective
285: Input Parameters:
286: + comm - MPI communicator
287: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
288: This value should be the same as the local size used in creating the
289: y vector for the matrix-vector product y = Ax.
290: . n - This value should be the same as the local size used in creating the
291: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
292: calculated if `N` is given) For square matrices n is almost always `m`.
293: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
294: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
295: - diag - the diagonal value
297: Output Parameter:
298: . J - the diagonal matrix
300: Level: advanced
302: Notes:
303: Only supports square matrices with the same number of local rows and columns
305: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
306: @*/
307: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
308: {
309: PetscFunctionBegin;
310: PetscCall(MatCreate(comm, J));
311: PetscCall(MatSetSizes(*J, m, n, M, N));
312: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
313: PetscCall(MatShift(*J, diag));
314: PetscCall(MatSetUp(*J));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@
319: MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
321: Not collective
323: Input Parameter:
324: . mat - a `MATCONSTANTDIAGONAL`
326: Output Parameter:
327: . value - the scalar value
329: Level: developer
331: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
332: @*/
333: PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
334: {
335: PetscFunctionBegin;
336: PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
341: {
342: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
344: PetscFunctionBegin;
345: *value = ctx->diag;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*MC
350: MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
351: along the diagonal.
353: Level: advanced
355: .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
356: M*/
357: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
358: {
359: Mat_ConstantDiagonal *ctx;
361: PetscFunctionBegin;
362: PetscCall(PetscNew(&ctx));
363: ctx->diag = 0.0;
364: A->data = (void *)ctx;
366: A->assembled = PETSC_TRUE;
367: A->preallocated = PETSC_TRUE;
368: A->structurally_symmetric = PETSC_BOOL3_TRUE;
369: A->structural_symmetry_eternal = PETSC_TRUE;
370: A->symmetric = PETSC_BOOL3_TRUE;
371: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
372: A->symmetry_eternal = PETSC_TRUE;
374: A->ops->mult = MatMult_ConstantDiagonal;
375: A->ops->multadd = MatMultAdd_ConstantDiagonal;
376: A->ops->multtranspose = MatMult_ConstantDiagonal;
377: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
378: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
379: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
380: A->ops->solve = MatSolve_ConstantDiagonal;
381: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
382: A->ops->norm = MatNorm_ConstantDiagonal;
383: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
384: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
385: A->ops->getrow = MatGetRow_ConstantDiagonal;
386: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
387: A->ops->sor = MatSOR_ConstantDiagonal;
388: A->ops->shift = MatShift_ConstantDiagonal;
389: A->ops->scale = MatScale_ConstantDiagonal;
390: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
391: A->ops->view = MatView_ConstantDiagonal;
392: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
393: A->ops->destroy = MatDestroy_ConstantDiagonal;
394: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
395: A->ops->equal = MatEqual_ConstantDiagonal;
396: A->ops->axpy = MatAXPY_ConstantDiagonal;
397: A->ops->setrandom = MatSetRandom_ConstantDiagonal;
398: A->ops->conjugate = MatConjugate_ConstantDiagonal;
399: A->ops->transpose = MatTranspose_ConstantDiagonal;
401: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
402: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
407: {
408: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
410: PetscFunctionBegin;
411: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
412: else fact->factorerrortype = MAT_FACTOR_NOERROR;
413: fctx->diag = 1.0 / actx->diag;
414: fact->ops->solve = MatMult_ConstantDiagonal;
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
419: {
420: PetscFunctionBegin;
421: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
426: {
427: PetscFunctionBegin;
428: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }
432: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
433: {
434: PetscInt n = A->rmap->n, N = A->rmap->N;
436: PetscFunctionBegin;
437: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
439: (*B)->factortype = ftype;
440: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
441: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
442: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
443: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
445: (*B)->ops->shift = NULL;
446: (*B)->ops->scale = NULL;
447: (*B)->ops->mult = NULL;
448: (*B)->ops->sor = NULL;
449: (*B)->ops->zeroentries = NULL;
451: PetscCall(PetscFree((*B)->solvertype));
452: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }