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: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_constantdiagonal_C", NULL));
117: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_diagonal_C", NULL));
118: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_constantdiagonal_C", NULL));
119: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_anytype_C", NULL));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
124: {
125: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
126: PetscBool isascii;
128: PetscFunctionBegin;
129: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
130: if (isascii) {
131: PetscViewerFormat format;
133: PetscCall(PetscViewerGetFormat(viewer, &format));
134: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
135: if (PetscImaginaryPart(ctx->diag) == 0) {
136: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag)));
137: } else {
138: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
139: }
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
145: {
146: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
148: PetscFunctionBegin;
149: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
154: {
155: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
157: PetscFunctionBegin;
158: PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
163: {
164: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
166: PetscFunctionBegin;
167: PetscCall(VecSet(x, ctx->diag));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
172: {
173: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
175: PetscFunctionBegin;
176: ctx->diag += a;
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
181: {
182: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
184: PetscFunctionBegin;
185: ctx->diag *= a;
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
190: {
191: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
193: PetscFunctionBegin;
194: ctx->diag = 0.0;
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y)
199: {
200: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
202: PetscFunctionBegin;
203: ctx->diag = PetscConj(ctx->diag);
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout)
208: {
209: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
211: PetscFunctionBegin;
212: if (reuse == MAT_INPLACE_MATRIX) {
213: PetscLayout tmplayout = A->rmap;
215: A->rmap = A->cmap;
216: A->cmap = tmplayout;
217: } else {
218: if (reuse == MAT_INITIAL_MATRIX) {
219: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout));
220: } else {
221: PetscCall(MatZeroEntries(*matout));
222: PetscCall(MatShift(*matout, ctx->diag));
223: }
224: }
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand)
229: {
230: PetscMPIInt rank;
231: MPI_Comm comm;
232: PetscScalar v = 0.0;
233: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
235: PetscFunctionBegin;
236: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
237: PetscCallMPI(MPI_Comm_rank(comm, &rank));
238: if (!rank) PetscCall(PetscRandomGetValue(rand, &v));
239: PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm));
240: ctx->diag = v;
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
245: {
246: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
248: PetscFunctionBegin;
249: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
250: else matin->factorerrortype = MAT_FACTOR_NOERROR;
251: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
256: {
257: PetscFunctionBegin;
258: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
263: {
264: PetscFunctionBegin;
265: info->block_size = 1.0;
266: info->nz_allocated = 1.0;
267: info->nz_used = 1.0;
268: info->nz_unneeded = 0.0;
269: info->assemblies = A->num_ass;
270: info->mallocs = 0.0;
271: info->memory = 0; /* REVIEW ME */
272: if (A->factortype) {
273: info->fill_ratio_given = 1.0;
274: info->fill_ratio_needed = 1.0;
275: info->factor_mallocs = 0.0;
276: } else {
277: info->fill_ratio_given = 0;
278: info->fill_ratio_needed = 0;
279: info->factor_mallocs = 0;
280: }
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
287: Collective
289: Input Parameters:
290: + comm - MPI communicator
291: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
292: This value should be the same as the local size used in creating the
293: y vector for the matrix-vector product y = Ax.
294: . n - This value should be the same as the local size used in creating the
295: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
296: calculated if `N` is given) For square matrices n is almost always `m`.
297: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
298: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
299: - diag - the diagonal value
301: Output Parameter:
302: . J - the diagonal matrix
304: Level: advanced
306: Notes:
307: Only supports square matrices with the same number of local rows and columns
309: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
310: @*/
311: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
312: {
313: PetscFunctionBegin;
314: PetscCall(MatCreate(comm, J));
315: PetscCall(MatSetSizes(*J, m, n, M, N));
316: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
317: PetscCall(MatShift(*J, diag));
318: PetscCall(MatSetUp(*J));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
325: Not collective
327: Input Parameter:
328: . mat - a `MATCONSTANTDIAGONAL`
330: Output Parameter:
331: . value - the scalar value
333: Level: developer
335: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
336: @*/
337: PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
338: {
339: PetscFunctionBegin;
340: PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
341: PetscFunctionReturn(PETSC_SUCCESS);
342: }
344: static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
345: {
346: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
348: PetscFunctionBegin;
349: *value = ctx->diag;
350: PetscFunctionReturn(PETSC_SUCCESS);
351: }
353: /* PtAP for constantdiagonal * constantdiagonal: C = alpha * beta^2 * I */
354: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_ConstDiag(Mat C)
355: {
356: Mat A = C->product->A, P = C->product->B;
357: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *p = (Mat_ConstantDiagonal *)P->data, *c = (Mat_ConstantDiagonal *)C->data;
359: PetscFunctionBegin;
360: MatCheckProduct(C, 1);
361: c->diag = a->diag * p->diag * p->diag;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_ConstDiag(Mat C)
366: {
367: Mat P = C->product->B;
369: PetscFunctionBegin;
370: MatCheckProduct(C, 1);
371: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
372: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
373: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
374: C->assembled = PETSC_TRUE;
375: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_ConstDiag;
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: /* PtAP for constantdiagonal A and diagonal P: C_i = alpha * p_i^2 */
380: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Diagonal(Mat C)
381: {
382: Mat A = C->product->A, P = C->product->B;
383: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
384: Vec pdiag, cdiag;
386: PetscFunctionBegin;
387: MatCheckProduct(C, 1);
388: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
389: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
390: PetscCall(VecPointwiseMult(cdiag, pdiag, pdiag));
391: PetscCall(VecScale(cdiag, a->diag));
392: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
393: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Diagonal(Mat C)
398: {
399: Mat P = C->product->B;
400: Vec pdiag, cdiag;
402: PetscFunctionBegin;
403: MatCheckProduct(C, 1);
404: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
405: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
406: PetscCall(MatSetType(C, MATDIAGONAL));
407: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
408: PetscCall(VecDuplicate(pdiag, &cdiag));
409: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
410: PetscCall(VecDestroy(&cdiag));
411: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
412: C->assembled = PETSC_TRUE;
413: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Diagonal;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /* PtAP for diagonal A and constantdiagonal P: C_i = beta^2 * a_i */
418: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_ConstDiag(Mat C)
419: {
420: Mat A = C->product->A, P = C->product->B;
421: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
422: Vec adiag, cdiag;
424: PetscFunctionBegin;
425: MatCheckProduct(C, 1);
426: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
427: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
428: PetscCall(VecCopy(adiag, cdiag));
429: PetscCall(VecScale(cdiag, p->diag * p->diag));
430: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
431: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_ConstDiag(Mat C)
436: {
437: Mat A = C->product->A, P = C->product->B;
438: Vec adiag, cdiag;
440: PetscFunctionBegin;
441: MatCheckProduct(C, 1);
442: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
443: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
444: PetscCall(MatSetType(C, MATDIAGONAL));
445: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
446: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
447: PetscCall(VecDuplicate(adiag, &cdiag));
448: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
449: PetscCall(VecDestroy(&cdiag));
450: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
451: C->assembled = PETSC_TRUE;
452: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_ConstDiag;
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: /* PtAP for any (non-diagonal) A and constantdiagonal P: C = beta^2 * A */
457: static PetscErrorCode MatProductNumeric_PtAP_Anytype_ConstDiag(Mat C)
458: {
459: Mat A = C->product->A, P = C->product->B;
460: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
462: PetscFunctionBegin;
463: MatCheckProduct(C, 1);
464: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
465: PetscCall(MatScale(C, p->diag * p->diag));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: static PetscErrorCode MatProductSymbolic_PtAP_Anytype_ConstDiag(Mat C)
470: {
471: Mat A = C->product->A;
472: Mat_Product *product = C->product;
473: Mat Cwork;
475: PetscFunctionBegin;
476: MatCheckProduct(C, 1);
477: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
478: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
479: C->product = NULL;
480: PetscCall(MatHeaderReplace(C, &Cwork));
481: C->product = product;
482: C->ops->productnumeric = MatProductNumeric_PtAP_Anytype_ConstDiag;
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: /* PtAP for constantdiagonal A and any non-diagonal P: C = alpha * P^T * P */
487: typedef struct {
488: Mat PtP; /* P^T * P result via MatProduct AtB */
489: PetscObjectState pnnzstate; /* P's nonzero state when inner symbolic was last built */
490: } MatProductCtx_PtAP_ConstDiag_Anytype;
492: static PetscErrorCode MatProductCtxDestroy_PtAP_ConstDiag_Anytype(PetscCtxRt data)
493: {
494: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = *(MatProductCtx_PtAP_ConstDiag_Anytype **)data;
496: PetscFunctionBegin;
497: PetscCall(MatDestroy(&ctx->PtP));
498: PetscCall(PetscFree(ctx));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Anytype(Mat C)
503: {
504: Mat_Product *product = C->product;
505: Mat A = product->A, P = product->B;
506: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = (MatProductCtx_PtAP_ConstDiag_Anytype *)product->data;
507: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
508: PetscObjectState pnnzstate;
510: PetscFunctionBegin;
511: MatCheckProduct(C, 1);
512: /* Rebuild inner symbolic if P's nonzero structure has changed */
513: PetscCall(MatGetNonzeroState(P, &pnnzstate));
514: if (pnnzstate != ctx->pnnzstate) {
515: PetscCall(MatDestroy(&ctx->PtP));
516: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
517: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
518: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
519: PetscCall(MatProductSetFromOptions(ctx->PtP));
520: PetscCall(MatProductSymbolic(ctx->PtP));
521: ctx->pnnzstate = pnnzstate;
522: }
523: /* Compute P^T * P */
524: PetscCall(MatProductNumeric(ctx->PtP));
525: PetscCall(MatCopy(ctx->PtP, C, SAME_NONZERO_PATTERN));
526: PetscCall(MatScale(C, a->diag));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Anytype(Mat C)
531: {
532: Mat_Product *product = C->product;
533: Mat P = product->B;
534: MatProductCtx_PtAP_ConstDiag_Anytype *ctx;
535: Mat Cwork;
537: PetscFunctionBegin;
538: MatCheckProduct(C, 1);
539: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
540: PetscCall(PetscNew(&ctx));
542: /* PtP = P^T * P (symbolic) */
543: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
544: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
545: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
546: PetscCall(MatProductSetFromOptions(ctx->PtP));
547: PetscCall(MatProductSymbolic(ctx->PtP));
549: /* Record P's nonzero state so numeric phase can detect structural changes */
550: PetscCall(MatGetNonzeroState(P, &ctx->pnnzstate));
552: /* Set up C with the same structure as PtP */
553: PetscCall(MatDuplicate(ctx->PtP, MAT_DO_NOT_COPY_VALUES, &Cwork));
554: C->product = NULL;
555: PetscCall(MatHeaderReplace(C, &Cwork));
556: C->product = product;
557: C->assembled = PETSC_TRUE;
558: product->data = ctx;
559: product->destroy = MatProductCtxDestroy_PtAP_ConstDiag_Anytype;
560: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Anytype;
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode MatProductSetFromOptions_ConstDiag_ConstDiag(Mat C)
565: {
566: Mat_Product *product = C->product;
568: PetscFunctionBegin;
569: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_ConstDiag;
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Diagonal(Mat C)
574: {
575: Mat_Product *product = C->product;
577: PetscFunctionBegin;
578: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Diagonal;
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: static PetscErrorCode MatProductSetFromOptions_Diagonal_ConstDiag(Mat C)
583: {
584: Mat_Product *product = C->product;
586: PetscFunctionBegin;
587: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_ConstDiag;
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Anytype(Mat C)
592: {
593: Mat_Product *product = C->product;
594: PetscBool Acdiag, Bcdiag;
596: PetscFunctionBegin;
597: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATCONSTANTDIAGONAL, &Acdiag));
598: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, MATCONSTANTDIAGONAL, &Bcdiag));
599: if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
600: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Anytype_ConstDiag;
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*MC
605: MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
606: along the diagonal.
608: Level: advanced
610: .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
611: M*/
612: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
613: {
614: Mat_ConstantDiagonal *ctx;
616: PetscFunctionBegin;
617: PetscCall(PetscNew(&ctx));
618: ctx->diag = 0.0;
619: A->data = (void *)ctx;
621: A->assembled = PETSC_TRUE;
622: A->preallocated = PETSC_TRUE;
623: A->structurally_symmetric = PETSC_BOOL3_TRUE;
624: A->structural_symmetry_eternal = PETSC_TRUE;
625: A->symmetric = PETSC_BOOL3_TRUE;
626: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
627: A->symmetry_eternal = PETSC_TRUE;
629: A->ops->mult = MatMult_ConstantDiagonal;
630: A->ops->multadd = MatMultAdd_ConstantDiagonal;
631: A->ops->multtranspose = MatMult_ConstantDiagonal;
632: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
633: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
634: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
635: A->ops->solve = MatSolve_ConstantDiagonal;
636: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
637: A->ops->norm = MatNorm_ConstantDiagonal;
638: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
639: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
640: A->ops->getrow = MatGetRow_ConstantDiagonal;
641: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
642: A->ops->sor = MatSOR_ConstantDiagonal;
643: A->ops->shift = MatShift_ConstantDiagonal;
644: A->ops->scale = MatScale_ConstantDiagonal;
645: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
646: A->ops->view = MatView_ConstantDiagonal;
647: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
648: A->ops->destroy = MatDestroy_ConstantDiagonal;
649: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
650: A->ops->equal = MatEqual_ConstantDiagonal;
651: A->ops->axpy = MatAXPY_ConstantDiagonal;
652: A->ops->setrandom = MatSetRandom_ConstantDiagonal;
653: A->ops->conjugate = MatConjugate_ConstantDiagonal;
654: A->ops->transpose = MatTranspose_ConstantDiagonal;
656: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
657: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
658: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_constantdiagonal_C", MatProductSetFromOptions_ConstDiag_ConstDiag));
659: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_diagonal_C", MatProductSetFromOptions_ConstDiag_Diagonal));
660: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_constantdiagonal_C", MatProductSetFromOptions_Diagonal_ConstDiag));
661: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_ConstDiag_Anytype));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
666: {
667: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
669: PetscFunctionBegin;
670: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
671: else fact->factorerrortype = MAT_FACTOR_NOERROR;
672: fctx->diag = 1.0 / actx->diag;
673: fact->ops->solve = MatMult_ConstantDiagonal;
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
678: {
679: PetscFunctionBegin;
680: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
685: {
686: PetscFunctionBegin;
687: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
692: {
693: PetscInt n = A->rmap->n, N = A->rmap->N;
695: PetscFunctionBegin;
696: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
698: (*B)->factortype = ftype;
699: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
700: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
701: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
702: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
704: (*B)->ops->shift = NULL;
705: (*B)->ops->scale = NULL;
706: (*B)->ops->mult = NULL;
707: (*B)->ops->sor = NULL;
708: (*B)->ops->zeroentries = NULL;
710: PetscCall(PetscFree((*B)->solvertype));
711: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }