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: PetscCall(PetscObjectStateIncrease((PetscObject)A));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
256: {
257: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
259: PetscFunctionBegin;
260: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
261: else matin->factorerrortype = MAT_FACTOR_NOERROR;
262: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
267: {
268: PetscFunctionBegin;
269: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
274: {
275: PetscFunctionBegin;
276: info->block_size = 1.0;
277: info->nz_allocated = 1.0;
278: info->nz_used = 1.0;
279: info->nz_unneeded = 0.0;
280: info->assemblies = A->num_ass;
281: info->mallocs = 0.0;
282: info->memory = 0; /* REVIEW ME */
283: if (A->factortype) {
284: info->fill_ratio_given = 1.0;
285: info->fill_ratio_needed = 1.0;
286: info->factor_mallocs = 0.0;
287: } else {
288: info->fill_ratio_given = 0;
289: info->fill_ratio_needed = 0;
290: info->factor_mallocs = 0;
291: }
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@
296: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
298: Collective
300: Input Parameters:
301: + comm - MPI communicator
302: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
303: This value should be the same as the local size used in creating the
304: y vector for the matrix-vector product y = Ax.
305: . n - This value should be the same as the local size used in creating the
306: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
307: calculated if `N` is given) For square matrices n is almost always `m`.
308: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
309: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
310: - diag - the diagonal value
312: Output Parameter:
313: . J - the diagonal matrix
315: Level: advanced
317: Notes:
318: Only supports square matrices with the same number of local rows and columns
320: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
321: @*/
322: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
323: {
324: PetscFunctionBegin;
325: PetscCall(MatCreate(comm, J));
326: PetscCall(MatSetSizes(*J, m, n, M, N));
327: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
328: PetscCall(MatShift(*J, diag));
329: PetscCall(MatSetUp(*J));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*@
334: MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
336: Not collective
338: Input Parameter:
339: . mat - a `MATCONSTANTDIAGONAL`
341: Output Parameter:
342: . value - the scalar value
344: Level: developer
346: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
347: @*/
348: PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
349: {
350: PetscFunctionBegin;
351: PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
356: {
357: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
359: PetscFunctionBegin;
360: *value = ctx->diag;
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode MatProductNumeric_ConstDiag_Dense(Mat C)
365: {
366: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)C->product->A->data;
368: PetscFunctionBegin;
369: MatCheckProduct(C, 1);
370: PetscCall(MatCopy(C->product->B, C, SAME_NONZERO_PATTERN));
371: PetscCall(MatScale(C, p->diag));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode MatProductSymbolic_ConstDiag_Dense(Mat C)
376: {
377: Mat A, B;
378: PetscInt n, N, m, M;
380: PetscFunctionBegin;
381: MatCheckProduct(C, 1);
382: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
383: A = C->product->A;
384: B = C->product->B;
385: PetscCall(MatGetLocalSize(B, NULL, &n));
386: PetscCall(MatGetSize(B, NULL, &N));
387: PetscCall(MatGetLocalSize(A, &m, NULL));
388: PetscCall(MatGetSize(A, &M, NULL));
389: PetscCall(MatSetSizes(C, m, n, M, N));
390: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
391: PetscCall(MatSetVecType(C, B->defaultvectype));
392: PetscCall(MatSetUp(C));
393: C->ops->productnumeric = MatProductNumeric_ConstDiag_Dense;
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /* PtAP for constantdiagonal * constantdiagonal: C = alpha * beta^2 * I */
398: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_ConstDiag(Mat C)
399: {
400: Mat A = C->product->A, P = C->product->B;
401: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *p = (Mat_ConstantDiagonal *)P->data, *c = (Mat_ConstantDiagonal *)C->data;
403: PetscFunctionBegin;
404: MatCheckProduct(C, 1);
405: c->diag = a->diag * p->diag * p->diag;
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_ConstDiag(Mat C)
410: {
411: Mat P = C->product->B;
413: PetscFunctionBegin;
414: MatCheckProduct(C, 1);
415: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
416: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
417: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
418: C->assembled = PETSC_TRUE;
419: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_ConstDiag;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /* PtAP for constantdiagonal A and diagonal P: C_i = alpha * p_i^2 */
424: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Diagonal(Mat C)
425: {
426: Mat A = C->product->A, P = C->product->B;
427: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
428: Vec pdiag, cdiag;
430: PetscFunctionBegin;
431: MatCheckProduct(C, 1);
432: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
433: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
434: PetscCall(VecPointwiseMult(cdiag, pdiag, pdiag));
435: PetscCall(VecScale(cdiag, a->diag));
436: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
437: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Diagonal(Mat C)
442: {
443: Mat P = C->product->B;
444: Vec pdiag, cdiag;
446: PetscFunctionBegin;
447: MatCheckProduct(C, 1);
448: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
449: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
450: PetscCall(MatSetType(C, MATDIAGONAL));
451: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
452: PetscCall(VecDuplicate(pdiag, &cdiag));
453: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
454: PetscCall(VecDestroy(&cdiag));
455: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
456: C->assembled = PETSC_TRUE;
457: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Diagonal;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: /* PtAP for diagonal A and constantdiagonal P: C_i = beta^2 * a_i */
462: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_ConstDiag(Mat C)
463: {
464: Mat A = C->product->A, P = C->product->B;
465: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
466: Vec adiag, cdiag;
468: PetscFunctionBegin;
469: MatCheckProduct(C, 1);
470: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
471: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
472: PetscCall(VecCopy(adiag, cdiag));
473: PetscCall(VecScale(cdiag, p->diag * p->diag));
474: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
475: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_ConstDiag(Mat C)
480: {
481: Mat A = C->product->A, P = C->product->B;
482: Vec adiag, cdiag;
484: PetscFunctionBegin;
485: MatCheckProduct(C, 1);
486: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
487: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
488: PetscCall(MatSetType(C, MATDIAGONAL));
489: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
490: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
491: PetscCall(VecDuplicate(adiag, &cdiag));
492: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
493: PetscCall(VecDestroy(&cdiag));
494: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
495: C->assembled = PETSC_TRUE;
496: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_ConstDiag;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /* PtAP for any (non-diagonal) A and constantdiagonal P: C = beta^2 * A */
501: static PetscErrorCode MatProductNumeric_PtAP_Anytype_ConstDiag(Mat C)
502: {
503: Mat A = C->product->A, P = C->product->B;
504: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
506: PetscFunctionBegin;
507: MatCheckProduct(C, 1);
508: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
509: PetscCall(MatScale(C, p->diag * p->diag));
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: static PetscErrorCode MatProductSymbolic_PtAP_Anytype_ConstDiag(Mat C)
514: {
515: Mat A = C->product->A;
516: Mat_Product *product = C->product;
517: Mat Cwork;
519: PetscFunctionBegin;
520: MatCheckProduct(C, 1);
521: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
522: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
523: C->product = NULL;
524: PetscCall(MatHeaderReplace(C, &Cwork));
525: C->product = product;
526: C->ops->productnumeric = MatProductNumeric_PtAP_Anytype_ConstDiag;
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /* PtAP for constantdiagonal A and any non-diagonal P: C = alpha * P^T * P */
531: typedef struct {
532: Mat PtP; /* P^T * P result via MatProduct AtB */
533: PetscObjectState pnnzstate; /* P's nonzero state when inner symbolic was last built */
534: } MatProductCtx_PtAP_ConstDiag_Anytype;
536: static PetscErrorCode MatProductCtxDestroy_PtAP_ConstDiag_Anytype(PetscCtxRt data)
537: {
538: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = *(MatProductCtx_PtAP_ConstDiag_Anytype **)data;
540: PetscFunctionBegin;
541: PetscCall(MatDestroy(&ctx->PtP));
542: PetscCall(PetscFree(ctx));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Anytype(Mat C)
547: {
548: Mat_Product *product = C->product;
549: Mat A = product->A, P = product->B;
550: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = (MatProductCtx_PtAP_ConstDiag_Anytype *)product->data;
551: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
552: PetscObjectState pnnzstate;
554: PetscFunctionBegin;
555: MatCheckProduct(C, 1);
556: /* Rebuild inner symbolic if P's nonzero structure has changed */
557: PetscCall(MatGetNonzeroState(P, &pnnzstate));
558: if (pnnzstate != ctx->pnnzstate) {
559: PetscCall(MatDestroy(&ctx->PtP));
560: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
561: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
562: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
563: PetscCall(MatProductSetFromOptions(ctx->PtP));
564: PetscCall(MatProductSymbolic(ctx->PtP));
565: ctx->pnnzstate = pnnzstate;
566: }
567: /* Compute P^T * P */
568: PetscCall(MatProductNumeric(ctx->PtP));
569: PetscCall(MatCopy(ctx->PtP, C, SAME_NONZERO_PATTERN));
570: PetscCall(MatScale(C, a->diag));
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }
574: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Anytype(Mat C)
575: {
576: Mat_Product *product = C->product;
577: Mat P = product->B;
578: MatProductCtx_PtAP_ConstDiag_Anytype *ctx;
579: Mat Cwork;
581: PetscFunctionBegin;
582: MatCheckProduct(C, 1);
583: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
584: PetscCall(PetscNew(&ctx));
586: /* PtP = P^T * P (symbolic) */
587: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
588: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
589: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
590: PetscCall(MatProductSetFromOptions(ctx->PtP));
591: PetscCall(MatProductSymbolic(ctx->PtP));
593: /* Record P's nonzero state so numeric phase can detect structural changes */
594: PetscCall(MatGetNonzeroState(P, &ctx->pnnzstate));
596: /* Set up C with the same structure as PtP */
597: PetscCall(MatDuplicate(ctx->PtP, MAT_DO_NOT_COPY_VALUES, &Cwork));
598: C->product = NULL;
599: PetscCall(MatHeaderReplace(C, &Cwork));
600: C->product = product;
601: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
602: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
603: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
604: product->data = ctx;
605: product->destroy = MatProductCtxDestroy_PtAP_ConstDiag_Anytype;
606: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Anytype;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /* AB for MATCONSTANTDIAGONAL A and any non-diagonal B: C = alpha * B */
611: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Anytype(Mat C)
612: {
613: Mat A = C->product->A, B = C->product->B;
614: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
616: PetscFunctionBegin;
617: MatCheckProduct(C, 1);
618: PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
619: PetscCall(MatScale(C, a->diag));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Anytype(Mat C)
624: {
625: Mat B = C->product->B;
626: Mat_Product *product = C->product;
627: Mat Cwork;
629: PetscFunctionBegin;
630: MatCheckProduct(C, 1);
631: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
632: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &Cwork));
633: C->product = NULL;
634: PetscCall(MatHeaderReplace(C, &Cwork));
635: C->product = product;
636: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Anytype;
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /* AB for any non-diagonal A and MATCONSTANTDIAGONAL B: C = beta * A */
641: static PetscErrorCode MatProductNumeric_AB_Anytype_ConstDiag(Mat C)
642: {
643: Mat A = C->product->A, B = C->product->B;
644: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
646: PetscFunctionBegin;
647: MatCheckProduct(C, 1);
648: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
649: PetscCall(MatScale(C, b->diag));
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode MatProductSymbolic_AB_Anytype_ConstDiag(Mat C)
654: {
655: Mat A = C->product->A;
656: Mat_Product *product = C->product;
657: Mat Cwork;
659: PetscFunctionBegin;
660: MatCheckProduct(C, 1);
661: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
662: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
663: C->product = NULL;
664: PetscCall(MatHeaderReplace(C, &Cwork));
665: C->product = product;
666: C->ops->productnumeric = MatProductNumeric_AB_Anytype_ConstDiag;
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: /* AB for MATCONSTANTDIAGONAL * MATCONSTANTDIAGONAL: C = alpha * beta * I */
671: static PetscErrorCode MatProductNumeric_AB_ConstDiag_ConstDiag(Mat C)
672: {
673: Mat A = C->product->A, B = C->product->B;
674: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *b = (Mat_ConstantDiagonal *)B->data, *c = (Mat_ConstantDiagonal *)C->data;
676: PetscFunctionBegin;
677: MatCheckProduct(C, 1);
678: c->diag = a->diag * b->diag;
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_ConstDiag(Mat C)
683: {
684: Mat A = C->product->A, B = C->product->B;
686: PetscFunctionBegin;
687: MatCheckProduct(C, 1);
688: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
689: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
690: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
691: C->assembled = PETSC_TRUE;
692: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_ConstDiag;
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: /* AB for MATCONSTANTDIAGONAL A and MATDIAGONAL B: C_i = alpha * b_i */
697: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Diagonal(Mat C)
698: {
699: Mat A = C->product->A, B = C->product->B;
700: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
701: Vec bdiag, cdiag;
703: PetscFunctionBegin;
704: MatCheckProduct(C, 1);
705: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
706: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
707: PetscCall(VecCopy(bdiag, cdiag));
708: PetscCall(VecScale(cdiag, a->diag));
709: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
710: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
711: PetscFunctionReturn(PETSC_SUCCESS);
712: }
714: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Diagonal(Mat C)
715: {
716: Mat A = C->product->A, B = C->product->B;
717: Vec bdiag, cdiag;
719: PetscFunctionBegin;
720: MatCheckProduct(C, 1);
721: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
722: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
723: PetscCall(MatSetType(C, MATDIAGONAL));
724: /* Duplicate B's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
725: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
726: PetscCall(VecDuplicate(bdiag, &cdiag));
727: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
728: PetscCall(VecDestroy(&cdiag));
729: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
730: C->assembled = PETSC_TRUE;
731: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Diagonal;
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: /* AB for MATDIAGONAL A and MATCONSTANTDIAGONAL B: C_i = beta * a_i */
736: static PetscErrorCode MatProductNumeric_AB_Diagonal_ConstDiag(Mat C)
737: {
738: Mat A = C->product->A, B = C->product->B;
739: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
740: Vec adiag, cdiag;
742: PetscFunctionBegin;
743: MatCheckProduct(C, 1);
744: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
745: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
746: PetscCall(VecCopy(adiag, cdiag));
747: PetscCall(VecScale(cdiag, b->diag));
748: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
749: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: static PetscErrorCode MatProductSymbolic_AB_Diagonal_ConstDiag(Mat C)
754: {
755: Mat A = C->product->A, B = C->product->B;
756: Vec adiag, cdiag;
758: PetscFunctionBegin;
759: MatCheckProduct(C, 1);
760: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
761: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
762: PetscCall(MatSetType(C, MATDIAGONAL));
763: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
764: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
765: PetscCall(VecDuplicate(adiag, &cdiag));
766: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
767: PetscCall(VecDestroy(&cdiag));
768: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
769: C->assembled = PETSC_TRUE;
770: C->ops->productnumeric = MatProductNumeric_AB_Diagonal_ConstDiag;
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: static PetscErrorCode MatProductSetFromOptions_ConstDiag_ConstDiag(Mat C)
775: {
776: Mat_Product *product = C->product;
778: PetscFunctionBegin;
779: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_ConstDiag;
780: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_ConstDiag;
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Diagonal(Mat C)
785: {
786: Mat_Product *product = C->product;
788: PetscFunctionBegin;
789: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Diagonal;
790: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Diagonal;
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode MatProductSetFromOptions_Diagonal_ConstDiag(Mat C)
795: {
796: Mat_Product *product = C->product;
798: PetscFunctionBegin;
799: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_ConstDiag;
800: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_Diagonal_ConstDiag;
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense_AB(Mat C)
805: {
806: PetscFunctionBegin;
807: C->ops->productsymbolic = MatProductSymbolic_ConstDiag_Dense;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense(Mat C)
812: {
813: Mat_Product *product = C->product;
815: PetscFunctionBegin;
816: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_ConstDiag_Dense_AB(C));
817: else if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Anytype(Mat C)
822: {
823: Mat_Product *product = C->product;
824: PetscBool Acdiag, Bcdiag;
826: PetscFunctionBegin;
827: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATCONSTANTDIAGONAL, &Acdiag));
828: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, MATCONSTANTDIAGONAL, &Bcdiag));
829: if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
830: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Anytype_ConstDiag;
831: else if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Anytype;
832: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_Anytype_ConstDiag;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /*MC
837: MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
838: along the diagonal.
840: Level: advanced
842: .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
843: M*/
844: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
845: {
846: Mat_ConstantDiagonal *ctx;
848: PetscFunctionBegin;
849: PetscCall(PetscNew(&ctx));
850: ctx->diag = 0.0;
851: A->data = (void *)ctx;
853: A->assembled = PETSC_TRUE;
854: A->preallocated = PETSC_TRUE;
855: A->structurally_symmetric = PETSC_BOOL3_TRUE;
856: A->structural_symmetry_eternal = PETSC_TRUE;
857: A->symmetric = PETSC_BOOL3_TRUE;
858: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
859: A->symmetry_eternal = PETSC_TRUE;
861: A->ops->mult = MatMult_ConstantDiagonal;
862: A->ops->multadd = MatMultAdd_ConstantDiagonal;
863: A->ops->multtranspose = MatMult_ConstantDiagonal;
864: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
865: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
866: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
867: A->ops->solve = MatSolve_ConstantDiagonal;
868: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
869: A->ops->norm = MatNorm_ConstantDiagonal;
870: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
871: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
872: A->ops->getrow = MatGetRow_ConstantDiagonal;
873: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
874: A->ops->sor = MatSOR_ConstantDiagonal;
875: A->ops->shift = MatShift_ConstantDiagonal;
876: A->ops->scale = MatScale_ConstantDiagonal;
877: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
878: A->ops->view = MatView_ConstantDiagonal;
879: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
880: A->ops->destroy = MatDestroy_ConstantDiagonal;
881: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
882: A->ops->equal = MatEqual_ConstantDiagonal;
883: A->ops->axpy = MatAXPY_ConstantDiagonal;
884: A->ops->setrandom = MatSetRandom_ConstantDiagonal;
885: A->ops->conjugate = MatConjugate_ConstantDiagonal;
886: A->ops->transpose = MatTranspose_ConstantDiagonal;
888: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
889: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
890: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_constantdiagonal_C", MatProductSetFromOptions_ConstDiag_ConstDiag));
891: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_diagonal_C", MatProductSetFromOptions_ConstDiag_Diagonal));
892: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_constantdiagonal_C", MatProductSetFromOptions_Diagonal_ConstDiag));
893: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdense_C", MatProductSetFromOptions_ConstDiag_Dense));
894: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidense_C", MatProductSetFromOptions_ConstDiag_Dense));
895: #if PetscDefined(HAVE_CUDA)
896: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdensecuda_C", MatProductSetFromOptions_ConstDiag_Dense));
897: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidensecuda_C", MatProductSetFromOptions_ConstDiag_Dense));
898: #endif
899: #if PetscDefined(HAVE_HIP)
900: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdensehip_C", MatProductSetFromOptions_ConstDiag_Dense));
901: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidensehip_C", MatProductSetFromOptions_ConstDiag_Dense));
902: #endif
903: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_ConstDiag_Anytype));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
908: {
909: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
911: PetscFunctionBegin;
912: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
913: else fact->factorerrortype = MAT_FACTOR_NOERROR;
914: fctx->diag = 1.0 / actx->diag;
915: fact->ops->solve = MatMult_ConstantDiagonal;
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
920: {
921: PetscFunctionBegin;
922: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
927: {
928: PetscFunctionBegin;
929: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
934: {
935: PetscInt n = A->rmap->n, N = A->rmap->N;
937: PetscFunctionBegin;
938: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
940: (*B)->factortype = ftype;
941: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
942: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
943: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
944: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
946: (*B)->ops->shift = NULL;
947: (*B)->ops->scale = NULL;
948: (*B)->ops->mult = NULL;
949: (*B)->ops->sor = NULL;
950: (*B)->ops->zeroentries = NULL;
952: PetscCall(PetscFree((*B)->solvertype));
953: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }