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_constantdiagonal_seqdense_C", NULL));
120: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_mpidense_C", NULL));
121: #if PetscDefined(HAVE_CUDA)
122: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_seqdensecuda_C", NULL));
123: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_mpidensecuda_C", NULL));
124: #endif
125: #if PetscDefined(HAVE_HIP)
126: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_seqdensehip_C", NULL));
127: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_constantdiagonal_mpidensehip_C", NULL));
128: #endif
129: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_anytype_C", NULL));
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 isascii;
138: PetscFunctionBegin;
139: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
140: if (isascii) {
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 (PetscImaginaryPart(ctx->diag) == 0) {
146: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag)));
147: } else {
148: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
149: }
150: }
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
155: {
156: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
158: PetscFunctionBegin;
159: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
164: {
165: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
167: PetscFunctionBegin;
168: PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
173: {
174: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
176: PetscFunctionBegin;
177: PetscCall(VecSet(x, ctx->diag));
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
182: {
183: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
185: PetscFunctionBegin;
186: ctx->diag += a;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
191: {
192: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
194: PetscFunctionBegin;
195: ctx->diag *= a;
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
200: {
201: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
203: PetscFunctionBegin;
204: ctx->diag = 0.0;
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y)
209: {
210: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
212: PetscFunctionBegin;
213: ctx->diag = PetscConj(ctx->diag);
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout)
218: {
219: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
221: PetscFunctionBegin;
222: if (reuse == MAT_INPLACE_MATRIX) {
223: PetscLayout tmplayout = A->rmap;
225: A->rmap = A->cmap;
226: A->cmap = tmplayout;
227: } else {
228: if (reuse == MAT_INITIAL_MATRIX) {
229: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout));
230: } else {
231: PetscCall(MatZeroEntries(*matout));
232: PetscCall(MatShift(*matout, ctx->diag));
233: }
234: }
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand)
239: {
240: PetscMPIInt rank;
241: MPI_Comm comm;
242: PetscScalar v = 0.0;
243: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
245: PetscFunctionBegin;
246: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
247: PetscCallMPI(MPI_Comm_rank(comm, &rank));
248: if (!rank) PetscCall(PetscRandomGetValue(rand, &v));
249: PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm));
250: ctx->diag = v;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
255: {
256: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
258: PetscFunctionBegin;
259: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
260: else matin->factorerrortype = MAT_FACTOR_NOERROR;
261: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
266: {
267: PetscFunctionBegin;
268: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
273: {
274: PetscFunctionBegin;
275: info->block_size = 1.0;
276: info->nz_allocated = 1.0;
277: info->nz_used = 1.0;
278: info->nz_unneeded = 0.0;
279: info->assemblies = A->num_ass;
280: info->mallocs = 0.0;
281: info->memory = 0; /* REVIEW ME */
282: if (A->factortype) {
283: info->fill_ratio_given = 1.0;
284: info->fill_ratio_needed = 1.0;
285: info->factor_mallocs = 0.0;
286: } else {
287: info->fill_ratio_given = 0;
288: info->fill_ratio_needed = 0;
289: info->factor_mallocs = 0;
290: }
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
297: Collective
299: Input Parameters:
300: + comm - MPI communicator
301: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
302: This value should be the same as the local size used in creating the
303: y vector for the matrix-vector product y = Ax.
304: . n - This value should be the same as the local size used in creating the
305: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
306: calculated if `N` is given) For square matrices n is almost always `m`.
307: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
308: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
309: - diag - the diagonal value
311: Output Parameter:
312: . J - the diagonal matrix
314: Level: advanced
316: Notes:
317: Only supports square matrices with the same number of local rows and columns
319: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
320: @*/
321: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
322: {
323: PetscFunctionBegin;
324: PetscCall(MatCreate(comm, J));
325: PetscCall(MatSetSizes(*J, m, n, M, N));
326: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
327: PetscCall(MatShift(*J, diag));
328: PetscCall(MatSetUp(*J));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*@
333: MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
335: Not collective
337: Input Parameter:
338: . mat - a `MATCONSTANTDIAGONAL`
340: Output Parameter:
341: . value - the scalar value
343: Level: developer
345: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
346: @*/
347: PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
348: {
349: PetscFunctionBegin;
350: PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
355: {
356: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
358: PetscFunctionBegin;
359: *value = ctx->diag;
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode MatProductNumeric_ConstDiag_Dense(Mat C)
364: {
365: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)C->product->A->data;
367: PetscFunctionBegin;
368: MatCheckProduct(C, 1);
369: PetscCall(MatCopy(C->product->B, C, SAME_NONZERO_PATTERN));
370: PetscCall(MatScale(C, p->diag));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode MatProductSymbolic_ConstDiag_Dense(Mat C)
375: {
376: Mat A, B;
377: PetscInt n, N, m, M;
379: PetscFunctionBegin;
380: MatCheckProduct(C, 1);
381: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
382: A = C->product->A;
383: B = C->product->B;
384: PetscCall(MatGetLocalSize(B, NULL, &n));
385: PetscCall(MatGetSize(B, NULL, &N));
386: PetscCall(MatGetLocalSize(A, &m, NULL));
387: PetscCall(MatGetSize(A, &M, NULL));
388: PetscCall(MatSetSizes(C, m, n, M, N));
389: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
390: PetscCall(MatSetVecType(C, B->defaultvectype));
391: PetscCall(MatSetUp(C));
392: C->ops->productnumeric = MatProductNumeric_ConstDiag_Dense;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /* PtAP for constantdiagonal * constantdiagonal: C = alpha * beta^2 * I */
397: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_ConstDiag(Mat C)
398: {
399: Mat A = C->product->A, P = C->product->B;
400: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *p = (Mat_ConstantDiagonal *)P->data, *c = (Mat_ConstantDiagonal *)C->data;
402: PetscFunctionBegin;
403: MatCheckProduct(C, 1);
404: c->diag = a->diag * p->diag * p->diag;
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_ConstDiag(Mat C)
409: {
410: Mat P = C->product->B;
412: PetscFunctionBegin;
413: MatCheckProduct(C, 1);
414: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
415: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
416: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
417: C->assembled = PETSC_TRUE;
418: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_ConstDiag;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /* PtAP for constantdiagonal A and diagonal P: C_i = alpha * p_i^2 */
423: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Diagonal(Mat C)
424: {
425: Mat A = C->product->A, P = C->product->B;
426: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
427: Vec pdiag, cdiag;
429: PetscFunctionBegin;
430: MatCheckProduct(C, 1);
431: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
432: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
433: PetscCall(VecPointwiseMult(cdiag, pdiag, pdiag));
434: PetscCall(VecScale(cdiag, a->diag));
435: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
436: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Diagonal(Mat C)
441: {
442: Mat P = C->product->B;
443: Vec pdiag, cdiag;
445: PetscFunctionBegin;
446: MatCheckProduct(C, 1);
447: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
448: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
449: PetscCall(MatSetType(C, MATDIAGONAL));
450: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
451: PetscCall(VecDuplicate(pdiag, &cdiag));
452: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
453: PetscCall(VecDestroy(&cdiag));
454: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
455: C->assembled = PETSC_TRUE;
456: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Diagonal;
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /* PtAP for diagonal A and constantdiagonal P: C_i = beta^2 * a_i */
461: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_ConstDiag(Mat C)
462: {
463: Mat A = C->product->A, P = C->product->B;
464: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
465: Vec adiag, cdiag;
467: PetscFunctionBegin;
468: MatCheckProduct(C, 1);
469: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
470: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
471: PetscCall(VecCopy(adiag, cdiag));
472: PetscCall(VecScale(cdiag, p->diag * p->diag));
473: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
474: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_ConstDiag(Mat C)
479: {
480: Mat A = C->product->A, P = C->product->B;
481: Vec adiag, cdiag;
483: PetscFunctionBegin;
484: MatCheckProduct(C, 1);
485: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
486: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
487: PetscCall(MatSetType(C, MATDIAGONAL));
488: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
489: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
490: PetscCall(VecDuplicate(adiag, &cdiag));
491: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
492: PetscCall(VecDestroy(&cdiag));
493: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
494: C->assembled = PETSC_TRUE;
495: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_ConstDiag;
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /* PtAP for any (non-diagonal) A and constantdiagonal P: C = beta^2 * A */
500: static PetscErrorCode MatProductNumeric_PtAP_Anytype_ConstDiag(Mat C)
501: {
502: Mat A = C->product->A, P = C->product->B;
503: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
505: PetscFunctionBegin;
506: MatCheckProduct(C, 1);
507: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
508: PetscCall(MatScale(C, p->diag * p->diag));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode MatProductSymbolic_PtAP_Anytype_ConstDiag(Mat C)
513: {
514: Mat A = C->product->A;
515: Mat_Product *product = C->product;
516: Mat Cwork;
518: PetscFunctionBegin;
519: MatCheckProduct(C, 1);
520: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
521: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
522: C->product = NULL;
523: PetscCall(MatHeaderReplace(C, &Cwork));
524: C->product = product;
525: C->ops->productnumeric = MatProductNumeric_PtAP_Anytype_ConstDiag;
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /* PtAP for constantdiagonal A and any non-diagonal P: C = alpha * P^T * P */
530: typedef struct {
531: Mat PtP; /* P^T * P result via MatProduct AtB */
532: PetscObjectState pnnzstate; /* P's nonzero state when inner symbolic was last built */
533: } MatProductCtx_PtAP_ConstDiag_Anytype;
535: static PetscErrorCode MatProductCtxDestroy_PtAP_ConstDiag_Anytype(PetscCtxRt data)
536: {
537: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = *(MatProductCtx_PtAP_ConstDiag_Anytype **)data;
539: PetscFunctionBegin;
540: PetscCall(MatDestroy(&ctx->PtP));
541: PetscCall(PetscFree(ctx));
542: PetscFunctionReturn(PETSC_SUCCESS);
543: }
545: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Anytype(Mat C)
546: {
547: Mat_Product *product = C->product;
548: Mat A = product->A, P = product->B;
549: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = (MatProductCtx_PtAP_ConstDiag_Anytype *)product->data;
550: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
551: PetscObjectState pnnzstate;
553: PetscFunctionBegin;
554: MatCheckProduct(C, 1);
555: /* Rebuild inner symbolic if P's nonzero structure has changed */
556: PetscCall(MatGetNonzeroState(P, &pnnzstate));
557: if (pnnzstate != ctx->pnnzstate) {
558: PetscCall(MatDestroy(&ctx->PtP));
559: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
560: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
561: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
562: PetscCall(MatProductSetFromOptions(ctx->PtP));
563: PetscCall(MatProductSymbolic(ctx->PtP));
564: ctx->pnnzstate = pnnzstate;
565: }
566: /* Compute P^T * P */
567: PetscCall(MatProductNumeric(ctx->PtP));
568: PetscCall(MatCopy(ctx->PtP, C, SAME_NONZERO_PATTERN));
569: PetscCall(MatScale(C, a->diag));
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Anytype(Mat C)
574: {
575: Mat_Product *product = C->product;
576: Mat P = product->B;
577: MatProductCtx_PtAP_ConstDiag_Anytype *ctx;
578: Mat Cwork;
580: PetscFunctionBegin;
581: MatCheckProduct(C, 1);
582: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
583: PetscCall(PetscNew(&ctx));
585: /* PtP = P^T * P (symbolic) */
586: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
587: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
588: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
589: PetscCall(MatProductSetFromOptions(ctx->PtP));
590: PetscCall(MatProductSymbolic(ctx->PtP));
592: /* Record P's nonzero state so numeric phase can detect structural changes */
593: PetscCall(MatGetNonzeroState(P, &ctx->pnnzstate));
595: /* Set up C with the same structure as PtP */
596: PetscCall(MatDuplicate(ctx->PtP, MAT_DO_NOT_COPY_VALUES, &Cwork));
597: C->product = NULL;
598: PetscCall(MatHeaderReplace(C, &Cwork));
599: C->product = product;
600: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
601: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
602: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
603: product->data = ctx;
604: product->destroy = MatProductCtxDestroy_PtAP_ConstDiag_Anytype;
605: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Anytype;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /* AB for MATCONSTANTDIAGONAL A and any non-diagonal B: C = alpha * B */
610: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Anytype(Mat C)
611: {
612: Mat A = C->product->A, B = C->product->B;
613: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
615: PetscFunctionBegin;
616: MatCheckProduct(C, 1);
617: PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
618: PetscCall(MatScale(C, a->diag));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Anytype(Mat C)
623: {
624: Mat B = C->product->B;
625: Mat_Product *product = C->product;
626: Mat Cwork;
628: PetscFunctionBegin;
629: MatCheckProduct(C, 1);
630: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
631: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &Cwork));
632: C->product = NULL;
633: PetscCall(MatHeaderReplace(C, &Cwork));
634: C->product = product;
635: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Anytype;
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /* AB for any non-diagonal A and MATCONSTANTDIAGONAL B: C = beta * A */
640: static PetscErrorCode MatProductNumeric_AB_Anytype_ConstDiag(Mat C)
641: {
642: Mat A = C->product->A, B = C->product->B;
643: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
645: PetscFunctionBegin;
646: MatCheckProduct(C, 1);
647: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
648: PetscCall(MatScale(C, b->diag));
649: PetscFunctionReturn(PETSC_SUCCESS);
650: }
652: static PetscErrorCode MatProductSymbolic_AB_Anytype_ConstDiag(Mat C)
653: {
654: Mat A = C->product->A;
655: Mat_Product *product = C->product;
656: Mat Cwork;
658: PetscFunctionBegin;
659: MatCheckProduct(C, 1);
660: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
661: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
662: C->product = NULL;
663: PetscCall(MatHeaderReplace(C, &Cwork));
664: C->product = product;
665: C->ops->productnumeric = MatProductNumeric_AB_Anytype_ConstDiag;
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /* AB for MATCONSTANTDIAGONAL * MATCONSTANTDIAGONAL: C = alpha * beta * I */
670: static PetscErrorCode MatProductNumeric_AB_ConstDiag_ConstDiag(Mat C)
671: {
672: Mat A = C->product->A, B = C->product->B;
673: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *b = (Mat_ConstantDiagonal *)B->data, *c = (Mat_ConstantDiagonal *)C->data;
675: PetscFunctionBegin;
676: MatCheckProduct(C, 1);
677: c->diag = a->diag * b->diag;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_ConstDiag(Mat C)
682: {
683: Mat A = C->product->A, B = C->product->B;
685: PetscFunctionBegin;
686: MatCheckProduct(C, 1);
687: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
688: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
689: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
690: C->assembled = PETSC_TRUE;
691: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_ConstDiag;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /* AB for MATCONSTANTDIAGONAL A and MATDIAGONAL B: C_i = alpha * b_i */
696: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Diagonal(Mat C)
697: {
698: Mat A = C->product->A, B = C->product->B;
699: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
700: Vec bdiag, cdiag;
702: PetscFunctionBegin;
703: MatCheckProduct(C, 1);
704: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
705: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
706: PetscCall(VecCopy(bdiag, cdiag));
707: PetscCall(VecScale(cdiag, a->diag));
708: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
709: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Diagonal(Mat C)
714: {
715: Mat A = C->product->A, B = C->product->B;
716: Vec bdiag, cdiag;
718: PetscFunctionBegin;
719: MatCheckProduct(C, 1);
720: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
721: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
722: PetscCall(MatSetType(C, MATDIAGONAL));
723: /* Duplicate B's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
724: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
725: PetscCall(VecDuplicate(bdiag, &cdiag));
726: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
727: PetscCall(VecDestroy(&cdiag));
728: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
729: C->assembled = PETSC_TRUE;
730: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Diagonal;
731: PetscFunctionReturn(PETSC_SUCCESS);
732: }
734: /* AB for MATDIAGONAL A and MATCONSTANTDIAGONAL B: C_i = beta * a_i */
735: static PetscErrorCode MatProductNumeric_AB_Diagonal_ConstDiag(Mat C)
736: {
737: Mat A = C->product->A, B = C->product->B;
738: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
739: Vec adiag, cdiag;
741: PetscFunctionBegin;
742: MatCheckProduct(C, 1);
743: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
744: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
745: PetscCall(VecCopy(adiag, cdiag));
746: PetscCall(VecScale(cdiag, b->diag));
747: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
748: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode MatProductSymbolic_AB_Diagonal_ConstDiag(Mat C)
753: {
754: Mat A = C->product->A, B = C->product->B;
755: Vec adiag, cdiag;
757: PetscFunctionBegin;
758: MatCheckProduct(C, 1);
759: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
760: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
761: PetscCall(MatSetType(C, MATDIAGONAL));
762: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
763: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
764: PetscCall(VecDuplicate(adiag, &cdiag));
765: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
766: PetscCall(VecDestroy(&cdiag));
767: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
768: C->assembled = PETSC_TRUE;
769: C->ops->productnumeric = MatProductNumeric_AB_Diagonal_ConstDiag;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: static PetscErrorCode MatProductSetFromOptions_ConstDiag_ConstDiag(Mat C)
774: {
775: Mat_Product *product = C->product;
777: PetscFunctionBegin;
778: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_ConstDiag;
779: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_ConstDiag;
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Diagonal(Mat C)
784: {
785: Mat_Product *product = C->product;
787: PetscFunctionBegin;
788: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Diagonal;
789: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Diagonal;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: static PetscErrorCode MatProductSetFromOptions_Diagonal_ConstDiag(Mat C)
794: {
795: Mat_Product *product = C->product;
797: PetscFunctionBegin;
798: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_ConstDiag;
799: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_Diagonal_ConstDiag;
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense_AB(Mat C)
804: {
805: PetscFunctionBegin;
806: C->ops->productsymbolic = MatProductSymbolic_ConstDiag_Dense;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense(Mat C)
811: {
812: Mat_Product *product = C->product;
814: PetscFunctionBegin;
815: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_ConstDiag_Dense_AB(C));
816: else if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Anytype(Mat C)
821: {
822: Mat_Product *product = C->product;
823: PetscBool Acdiag, Bcdiag;
825: PetscFunctionBegin;
826: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATCONSTANTDIAGONAL, &Acdiag));
827: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, MATCONSTANTDIAGONAL, &Bcdiag));
828: if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
829: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Anytype_ConstDiag;
830: else if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Anytype;
831: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_Anytype_ConstDiag;
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /*MC
836: MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
837: along the diagonal.
839: Level: advanced
841: .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
842: M*/
843: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
844: {
845: Mat_ConstantDiagonal *ctx;
847: PetscFunctionBegin;
848: PetscCall(PetscNew(&ctx));
849: ctx->diag = 0.0;
850: A->data = (void *)ctx;
852: A->assembled = PETSC_TRUE;
853: A->preallocated = PETSC_TRUE;
854: A->structurally_symmetric = PETSC_BOOL3_TRUE;
855: A->structural_symmetry_eternal = PETSC_TRUE;
856: A->symmetric = PETSC_BOOL3_TRUE;
857: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
858: A->symmetry_eternal = PETSC_TRUE;
860: A->ops->mult = MatMult_ConstantDiagonal;
861: A->ops->multadd = MatMultAdd_ConstantDiagonal;
862: A->ops->multtranspose = MatMult_ConstantDiagonal;
863: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
864: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
865: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
866: A->ops->solve = MatSolve_ConstantDiagonal;
867: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
868: A->ops->norm = MatNorm_ConstantDiagonal;
869: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
870: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
871: A->ops->getrow = MatGetRow_ConstantDiagonal;
872: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
873: A->ops->sor = MatSOR_ConstantDiagonal;
874: A->ops->shift = MatShift_ConstantDiagonal;
875: A->ops->scale = MatScale_ConstantDiagonal;
876: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
877: A->ops->view = MatView_ConstantDiagonal;
878: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
879: A->ops->destroy = MatDestroy_ConstantDiagonal;
880: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
881: A->ops->equal = MatEqual_ConstantDiagonal;
882: A->ops->axpy = MatAXPY_ConstantDiagonal;
883: A->ops->setrandom = MatSetRandom_ConstantDiagonal;
884: A->ops->conjugate = MatConjugate_ConstantDiagonal;
885: A->ops->transpose = MatTranspose_ConstantDiagonal;
887: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
888: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
889: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_constantdiagonal_C", MatProductSetFromOptions_ConstDiag_ConstDiag));
890: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_diagonal_C", MatProductSetFromOptions_ConstDiag_Diagonal));
891: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_constantdiagonal_C", MatProductSetFromOptions_Diagonal_ConstDiag));
892: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdense_C", MatProductSetFromOptions_ConstDiag_Dense));
893: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidense_C", MatProductSetFromOptions_ConstDiag_Dense));
894: #if PetscDefined(HAVE_CUDA)
895: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdensecuda_C", MatProductSetFromOptions_ConstDiag_Dense));
896: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidensecuda_C", MatProductSetFromOptions_ConstDiag_Dense));
897: #endif
898: #if PetscDefined(HAVE_HIP)
899: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdensehip_C", MatProductSetFromOptions_ConstDiag_Dense));
900: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidensehip_C", MatProductSetFromOptions_ConstDiag_Dense));
901: #endif
902: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_ConstDiag_Anytype));
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
907: {
908: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
910: PetscFunctionBegin;
911: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
912: else fact->factorerrortype = MAT_FACTOR_NOERROR;
913: fctx->diag = 1.0 / actx->diag;
914: fact->ops->solve = MatMult_ConstantDiagonal;
915: PetscFunctionReturn(PETSC_SUCCESS);
916: }
918: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
919: {
920: PetscFunctionBegin;
921: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
926: {
927: PetscFunctionBegin;
928: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
929: PetscFunctionReturn(PETSC_SUCCESS);
930: }
932: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
933: {
934: PetscInt n = A->rmap->n, N = A->rmap->N;
936: PetscFunctionBegin;
937: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
939: (*B)->factortype = ftype;
940: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
941: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
942: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
943: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
945: (*B)->ops->shift = NULL;
946: (*B)->ops->scale = NULL;
947: (*B)->ops->mult = NULL;
948: (*B)->ops->sor = NULL;
949: (*B)->ops->zeroentries = NULL;
951: PetscCall(PetscFree((*B)->solvertype));
952: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
953: PetscFunctionReturn(PETSC_SUCCESS);
954: }