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: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_anytype_C", NULL));
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: static PetscErrorCode MatView_ConstantDiagonal(Mat J, PetscViewer viewer)
126: {
127: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
128: PetscBool isascii;
130: PetscFunctionBegin;
131: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
132: if (isascii) {
133: PetscViewerFormat format;
135: PetscCall(PetscViewerGetFormat(viewer, &format));
136: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
137: if (PetscImaginaryPart(ctx->diag) == 0) {
138: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g\n", (double)PetscRealPart(ctx->diag)));
139: } else {
140: PetscCall(PetscViewerASCIIPrintf(viewer, "Diagonal value: %g + i %g\n", (double)PetscRealPart(ctx->diag), (double)PetscImaginaryPart(ctx->diag)));
141: }
142: }
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: static PetscErrorCode MatMult_ConstantDiagonal(Mat J, Vec x, Vec y)
147: {
148: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
150: PetscFunctionBegin;
151: PetscCall(VecAXPBY(y, ctx->diag, 0.0, x));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode MatMultHermitianTranspose_ConstantDiagonal(Mat J, Vec x, Vec y)
156: {
157: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
159: PetscFunctionBegin;
160: PetscCall(VecAXPBY(y, PetscConj(ctx->diag), 0.0, x));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode MatGetDiagonal_ConstantDiagonal(Mat J, Vec x)
165: {
166: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)J->data;
168: PetscFunctionBegin;
169: PetscCall(VecSet(x, ctx->diag));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode MatShift_ConstantDiagonal(Mat Y, PetscScalar a)
174: {
175: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
177: PetscFunctionBegin;
178: ctx->diag += a;
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: static PetscErrorCode MatScale_ConstantDiagonal(Mat Y, PetscScalar a)
183: {
184: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
186: PetscFunctionBegin;
187: ctx->diag *= a;
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: static PetscErrorCode MatZeroEntries_ConstantDiagonal(Mat Y)
192: {
193: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
195: PetscFunctionBegin;
196: ctx->diag = 0.0;
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode MatConjugate_ConstantDiagonal(Mat Y)
201: {
202: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)Y->data;
204: PetscFunctionBegin;
205: ctx->diag = PetscConj(ctx->diag);
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode MatTranspose_ConstantDiagonal(Mat A, MatReuse reuse, Mat *matout)
210: {
211: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
213: PetscFunctionBegin;
214: if (reuse == MAT_INPLACE_MATRIX) {
215: PetscLayout tmplayout = A->rmap;
217: A->rmap = A->cmap;
218: A->cmap = tmplayout;
219: } else {
220: if (reuse == MAT_INITIAL_MATRIX) {
221: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), A->cmap->n, A->rmap->n, A->cmap->N, A->rmap->N, ctx->diag, matout));
222: } else {
223: PetscCall(MatZeroEntries(*matout));
224: PetscCall(MatShift(*matout, ctx->diag));
225: }
226: }
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: static PetscErrorCode MatSetRandom_ConstantDiagonal(Mat A, PetscRandom rand)
231: {
232: PetscMPIInt rank;
233: MPI_Comm comm;
234: PetscScalar v = 0.0;
235: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)A->data;
237: PetscFunctionBegin;
238: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
239: PetscCallMPI(MPI_Comm_rank(comm, &rank));
240: if (!rank) PetscCall(PetscRandomGetValue(rand, &v));
241: PetscCallMPI(MPI_Bcast(&v, 1, MPIU_SCALAR, 0, comm));
242: ctx->diag = v;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static PetscErrorCode MatSolve_ConstantDiagonal(Mat matin, Vec b, Vec x)
247: {
248: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)matin->data;
250: PetscFunctionBegin;
251: if (ctx->diag == 0.0) matin->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
252: else matin->factorerrortype = MAT_FACTOR_NOERROR;
253: PetscCall(VecAXPBY(x, 1.0 / ctx->diag, 0.0, b));
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: static PetscErrorCode MatSOR_ConstantDiagonal(Mat matin, Vec x, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec y)
258: {
259: PetscFunctionBegin;
260: PetscCall(MatSolve_ConstantDiagonal(matin, x, y));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: static PetscErrorCode MatGetInfo_ConstantDiagonal(Mat A, MatInfoType flag, MatInfo *info)
265: {
266: PetscFunctionBegin;
267: info->block_size = 1.0;
268: info->nz_allocated = 1.0;
269: info->nz_used = 1.0;
270: info->nz_unneeded = 0.0;
271: info->assemblies = A->num_ass;
272: info->mallocs = 0.0;
273: info->memory = 0; /* REVIEW ME */
274: if (A->factortype) {
275: info->fill_ratio_given = 1.0;
276: info->fill_ratio_needed = 1.0;
277: info->factor_mallocs = 0.0;
278: } else {
279: info->fill_ratio_given = 0;
280: info->fill_ratio_needed = 0;
281: info->factor_mallocs = 0;
282: }
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: MatCreateConstantDiagonal - Creates a matrix with a uniform value along the diagonal
289: Collective
291: Input Parameters:
292: + comm - MPI communicator
293: . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
294: This value should be the same as the local size used in creating the
295: y vector for the matrix-vector product y = Ax.
296: . n - This value should be the same as the local size used in creating the
297: x vector for the matrix-vector product y = Ax. (or `PETSC_DECIDE` to have
298: calculated if `N` is given) For square matrices n is almost always `m`.
299: . M - number of global rows (or `PETSC_DETERMINE` to have calculated if m is given)
300: . N - number of global columns (or `PETSC_DETERMINE` to have calculated if n is given)
301: - diag - the diagonal value
303: Output Parameter:
304: . J - the diagonal matrix
306: Level: advanced
308: Notes:
309: Only supports square matrices with the same number of local rows and columns
311: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatGetFactor()`, `MatSolve()`
312: @*/
313: PetscErrorCode MatCreateConstantDiagonal(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar diag, Mat *J)
314: {
315: PetscFunctionBegin;
316: PetscCall(MatCreate(comm, J));
317: PetscCall(MatSetSizes(*J, m, n, M, N));
318: PetscCall(MatSetType(*J, MATCONSTANTDIAGONAL));
319: PetscCall(MatShift(*J, diag));
320: PetscCall(MatSetUp(*J));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: MatConstantDiagonalGetConstant - Get the scalar constant of a constant diagonal matrix
327: Not collective
329: Input Parameter:
330: . mat - a `MATCONSTANTDIAGONAL`
332: Output Parameter:
333: . value - the scalar value
335: Level: developer
337: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`
338: @*/
339: PetscErrorCode MatConstantDiagonalGetConstant(Mat mat, PetscScalar *value)
340: {
341: PetscFunctionBegin;
342: PetscUseMethod(mat, "MatConstantDiagonalGetConstant_C", (Mat, PetscScalar *), (mat, value));
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode MatConstantDiagonalGetConstant_ConstantDiagonal(Mat mat, PetscScalar *value)
347: {
348: Mat_ConstantDiagonal *ctx = (Mat_ConstantDiagonal *)mat->data;
350: PetscFunctionBegin;
351: *value = ctx->diag;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode MatProductNumeric_ConstDiag_Dense(Mat C)
356: {
357: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)C->product->A->data;
359: PetscFunctionBegin;
360: MatCheckProduct(C, 1);
361: PetscCall(MatCopy(C->product->B, C, SAME_NONZERO_PATTERN));
362: PetscCall(MatScale(C, p->diag));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static PetscErrorCode MatProductSymbolic_ConstDiag_Dense(Mat C)
367: {
368: Mat A, B;
369: PetscInt n, N, m, M;
371: PetscFunctionBegin;
372: MatCheckProduct(C, 1);
373: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
374: A = C->product->A;
375: B = C->product->B;
376: PetscCall(MatGetLocalSize(B, NULL, &n));
377: PetscCall(MatGetSize(B, NULL, &N));
378: PetscCall(MatGetLocalSize(A, &m, NULL));
379: PetscCall(MatGetSize(A, &M, NULL));
380: PetscCall(MatSetSizes(C, m, n, M, N));
381: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
382: PetscCall(MatSetUp(C));
383: C->ops->productnumeric = MatProductNumeric_ConstDiag_Dense;
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /* PtAP for constantdiagonal * constantdiagonal: C = alpha * beta^2 * I */
388: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_ConstDiag(Mat C)
389: {
390: Mat A = C->product->A, P = C->product->B;
391: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *p = (Mat_ConstantDiagonal *)P->data, *c = (Mat_ConstantDiagonal *)C->data;
393: PetscFunctionBegin;
394: MatCheckProduct(C, 1);
395: c->diag = a->diag * p->diag * p->diag;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_ConstDiag(Mat C)
400: {
401: Mat P = C->product->B;
403: PetscFunctionBegin;
404: MatCheckProduct(C, 1);
405: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
406: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
407: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
408: C->assembled = PETSC_TRUE;
409: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_ConstDiag;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /* PtAP for constantdiagonal A and diagonal P: C_i = alpha * p_i^2 */
414: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Diagonal(Mat C)
415: {
416: Mat A = C->product->A, P = C->product->B;
417: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
418: Vec pdiag, cdiag;
420: PetscFunctionBegin;
421: MatCheckProduct(C, 1);
422: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
423: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
424: PetscCall(VecPointwiseMult(cdiag, pdiag, pdiag));
425: PetscCall(VecScale(cdiag, a->diag));
426: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
427: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Diagonal(Mat C)
432: {
433: Mat P = C->product->B;
434: Vec pdiag, cdiag;
436: PetscFunctionBegin;
437: MatCheckProduct(C, 1);
438: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
439: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
440: PetscCall(MatSetType(C, MATDIAGONAL));
441: PetscCall(MatDiagonalGetDiagonal(P, &pdiag));
442: PetscCall(VecDuplicate(pdiag, &cdiag));
443: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
444: PetscCall(VecDestroy(&cdiag));
445: PetscCall(MatDiagonalRestoreDiagonal(P, &pdiag));
446: C->assembled = PETSC_TRUE;
447: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Diagonal;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /* PtAP for diagonal A and constantdiagonal P: C_i = beta^2 * a_i */
452: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_ConstDiag(Mat C)
453: {
454: Mat A = C->product->A, P = C->product->B;
455: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
456: Vec adiag, cdiag;
458: PetscFunctionBegin;
459: MatCheckProduct(C, 1);
460: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
461: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
462: PetscCall(VecCopy(adiag, cdiag));
463: PetscCall(VecScale(cdiag, p->diag * p->diag));
464: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
465: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_ConstDiag(Mat C)
470: {
471: Mat A = C->product->A, P = C->product->B;
472: Vec adiag, cdiag;
474: PetscFunctionBegin;
475: MatCheckProduct(C, 1);
476: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
477: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
478: PetscCall(MatSetType(C, MATDIAGONAL));
479: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
480: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
481: PetscCall(VecDuplicate(adiag, &cdiag));
482: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
483: PetscCall(VecDestroy(&cdiag));
484: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
485: C->assembled = PETSC_TRUE;
486: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_ConstDiag;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /* PtAP for any (non-diagonal) A and constantdiagonal P: C = beta^2 * A */
491: static PetscErrorCode MatProductNumeric_PtAP_Anytype_ConstDiag(Mat C)
492: {
493: Mat A = C->product->A, P = C->product->B;
494: Mat_ConstantDiagonal *p = (Mat_ConstantDiagonal *)P->data;
496: PetscFunctionBegin;
497: MatCheckProduct(C, 1);
498: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
499: PetscCall(MatScale(C, p->diag * p->diag));
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }
503: static PetscErrorCode MatProductSymbolic_PtAP_Anytype_ConstDiag(Mat C)
504: {
505: Mat A = C->product->A;
506: Mat_Product *product = C->product;
507: Mat Cwork;
509: PetscFunctionBegin;
510: MatCheckProduct(C, 1);
511: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
512: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
513: C->product = NULL;
514: PetscCall(MatHeaderReplace(C, &Cwork));
515: C->product = product;
516: C->ops->productnumeric = MatProductNumeric_PtAP_Anytype_ConstDiag;
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /* PtAP for constantdiagonal A and any non-diagonal P: C = alpha * P^T * P */
521: typedef struct {
522: Mat PtP; /* P^T * P result via MatProduct AtB */
523: PetscObjectState pnnzstate; /* P's nonzero state when inner symbolic was last built */
524: } MatProductCtx_PtAP_ConstDiag_Anytype;
526: static PetscErrorCode MatProductCtxDestroy_PtAP_ConstDiag_Anytype(PetscCtxRt data)
527: {
528: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = *(MatProductCtx_PtAP_ConstDiag_Anytype **)data;
530: PetscFunctionBegin;
531: PetscCall(MatDestroy(&ctx->PtP));
532: PetscCall(PetscFree(ctx));
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: static PetscErrorCode MatProductNumeric_PtAP_ConstDiag_Anytype(Mat C)
537: {
538: Mat_Product *product = C->product;
539: Mat A = product->A, P = product->B;
540: MatProductCtx_PtAP_ConstDiag_Anytype *ctx = (MatProductCtx_PtAP_ConstDiag_Anytype *)product->data;
541: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
542: PetscObjectState pnnzstate;
544: PetscFunctionBegin;
545: MatCheckProduct(C, 1);
546: /* Rebuild inner symbolic if P's nonzero structure has changed */
547: PetscCall(MatGetNonzeroState(P, &pnnzstate));
548: if (pnnzstate != ctx->pnnzstate) {
549: PetscCall(MatDestroy(&ctx->PtP));
550: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
551: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
552: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
553: PetscCall(MatProductSetFromOptions(ctx->PtP));
554: PetscCall(MatProductSymbolic(ctx->PtP));
555: ctx->pnnzstate = pnnzstate;
556: }
557: /* Compute P^T * P */
558: PetscCall(MatProductNumeric(ctx->PtP));
559: PetscCall(MatCopy(ctx->PtP, C, SAME_NONZERO_PATTERN));
560: PetscCall(MatScale(C, a->diag));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode MatProductSymbolic_PtAP_ConstDiag_Anytype(Mat C)
565: {
566: Mat_Product *product = C->product;
567: Mat P = product->B;
568: MatProductCtx_PtAP_ConstDiag_Anytype *ctx;
569: Mat Cwork;
571: PetscFunctionBegin;
572: MatCheckProduct(C, 1);
573: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
574: PetscCall(PetscNew(&ctx));
576: /* PtP = P^T * P (symbolic) */
577: PetscCall(MatProductCreate(P, P, NULL, &ctx->PtP));
578: PetscCall(MatProductSetType(ctx->PtP, MATPRODUCT_AtB));
579: PetscCall(MatProductSetFill(ctx->PtP, product->fill));
580: PetscCall(MatProductSetFromOptions(ctx->PtP));
581: PetscCall(MatProductSymbolic(ctx->PtP));
583: /* Record P's nonzero state so numeric phase can detect structural changes */
584: PetscCall(MatGetNonzeroState(P, &ctx->pnnzstate));
586: /* Set up C with the same structure as PtP */
587: PetscCall(MatDuplicate(ctx->PtP, MAT_DO_NOT_COPY_VALUES, &Cwork));
588: C->product = NULL;
589: PetscCall(MatHeaderReplace(C, &Cwork));
590: C->product = product;
591: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
592: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
593: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
594: product->data = ctx;
595: product->destroy = MatProductCtxDestroy_PtAP_ConstDiag_Anytype;
596: C->ops->productnumeric = MatProductNumeric_PtAP_ConstDiag_Anytype;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /* AB for MATCONSTANTDIAGONAL A and any non-diagonal B: C = alpha * B */
601: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Anytype(Mat C)
602: {
603: Mat A = C->product->A, B = C->product->B;
604: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
606: PetscFunctionBegin;
607: MatCheckProduct(C, 1);
608: PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
609: PetscCall(MatScale(C, a->diag));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Anytype(Mat C)
614: {
615: Mat B = C->product->B;
616: Mat_Product *product = C->product;
617: Mat Cwork;
619: PetscFunctionBegin;
620: MatCheckProduct(C, 1);
621: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
622: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &Cwork));
623: C->product = NULL;
624: PetscCall(MatHeaderReplace(C, &Cwork));
625: C->product = product;
626: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Anytype;
627: PetscFunctionReturn(PETSC_SUCCESS);
628: }
630: /* AB for any non-diagonal A and MATCONSTANTDIAGONAL B: C = beta * A */
631: static PetscErrorCode MatProductNumeric_AB_Anytype_ConstDiag(Mat C)
632: {
633: Mat A = C->product->A, B = C->product->B;
634: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
636: PetscFunctionBegin;
637: MatCheckProduct(C, 1);
638: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
639: PetscCall(MatScale(C, b->diag));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode MatProductSymbolic_AB_Anytype_ConstDiag(Mat C)
644: {
645: Mat A = C->product->A;
646: Mat_Product *product = C->product;
647: Mat Cwork;
649: PetscFunctionBegin;
650: MatCheckProduct(C, 1);
651: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
652: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
653: C->product = NULL;
654: PetscCall(MatHeaderReplace(C, &Cwork));
655: C->product = product;
656: C->ops->productnumeric = MatProductNumeric_AB_Anytype_ConstDiag;
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /* AB for MATCONSTANTDIAGONAL * MATCONSTANTDIAGONAL: C = alpha * beta * I */
661: static PetscErrorCode MatProductNumeric_AB_ConstDiag_ConstDiag(Mat C)
662: {
663: Mat A = C->product->A, B = C->product->B;
664: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data, *b = (Mat_ConstantDiagonal *)B->data, *c = (Mat_ConstantDiagonal *)C->data;
666: PetscFunctionBegin;
667: MatCheckProduct(C, 1);
668: c->diag = a->diag * b->diag;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_ConstDiag(Mat C)
673: {
674: Mat A = C->product->A, B = C->product->B;
676: PetscFunctionBegin;
677: MatCheckProduct(C, 1);
678: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
679: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
680: PetscCall(MatSetType(C, MATCONSTANTDIAGONAL));
681: C->assembled = PETSC_TRUE;
682: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_ConstDiag;
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /* AB for MATCONSTANTDIAGONAL A and MATDIAGONAL B: C_i = alpha * b_i */
687: static PetscErrorCode MatProductNumeric_AB_ConstDiag_Diagonal(Mat C)
688: {
689: Mat A = C->product->A, B = C->product->B;
690: Mat_ConstantDiagonal *a = (Mat_ConstantDiagonal *)A->data;
691: Vec bdiag, cdiag;
693: PetscFunctionBegin;
694: MatCheckProduct(C, 1);
695: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
696: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
697: PetscCall(VecCopy(bdiag, cdiag));
698: PetscCall(VecScale(cdiag, a->diag));
699: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
700: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: static PetscErrorCode MatProductSymbolic_AB_ConstDiag_Diagonal(Mat C)
705: {
706: Mat A = C->product->A, B = C->product->B;
707: Vec bdiag, cdiag;
709: PetscFunctionBegin;
710: MatCheckProduct(C, 1);
711: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
712: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
713: PetscCall(MatSetType(C, MATDIAGONAL));
714: /* Duplicate B's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
715: PetscCall(MatDiagonalGetDiagonal(B, &bdiag));
716: PetscCall(VecDuplicate(bdiag, &cdiag));
717: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
718: PetscCall(VecDestroy(&cdiag));
719: PetscCall(MatDiagonalRestoreDiagonal(B, &bdiag));
720: C->assembled = PETSC_TRUE;
721: C->ops->productnumeric = MatProductNumeric_AB_ConstDiag_Diagonal;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /* AB for MATDIAGONAL A and MATCONSTANTDIAGONAL B: C_i = beta * a_i */
726: static PetscErrorCode MatProductNumeric_AB_Diagonal_ConstDiag(Mat C)
727: {
728: Mat A = C->product->A, B = C->product->B;
729: Mat_ConstantDiagonal *b = (Mat_ConstantDiagonal *)B->data;
730: Vec adiag, cdiag;
732: PetscFunctionBegin;
733: MatCheckProduct(C, 1);
734: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
735: PetscCall(MatDiagonalGetDiagonal(C, &cdiag));
736: PetscCall(VecCopy(adiag, cdiag));
737: PetscCall(VecScale(cdiag, b->diag));
738: PetscCall(MatDiagonalRestoreDiagonal(C, &cdiag));
739: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
740: PetscFunctionReturn(PETSC_SUCCESS);
741: }
743: static PetscErrorCode MatProductSymbolic_AB_Diagonal_ConstDiag(Mat C)
744: {
745: Mat A = C->product->A, B = C->product->B;
746: Vec adiag, cdiag;
748: PetscFunctionBegin;
749: MatCheckProduct(C, 1);
750: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
751: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
752: PetscCall(MatSetType(C, MATDIAGONAL));
753: /* Duplicate A's diagonal Vec so C inherits the correct VecType (e.g., Kokkos, CUDA, HIP) */
754: PetscCall(MatDiagonalGetDiagonal(A, &adiag));
755: PetscCall(VecDuplicate(adiag, &cdiag));
756: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
757: PetscCall(VecDestroy(&cdiag));
758: PetscCall(MatDiagonalRestoreDiagonal(A, &adiag));
759: C->assembled = PETSC_TRUE;
760: C->ops->productnumeric = MatProductNumeric_AB_Diagonal_ConstDiag;
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: static PetscErrorCode MatProductSetFromOptions_ConstDiag_ConstDiag(Mat C)
765: {
766: Mat_Product *product = C->product;
768: PetscFunctionBegin;
769: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_ConstDiag;
770: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_ConstDiag;
771: PetscFunctionReturn(PETSC_SUCCESS);
772: }
774: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Diagonal(Mat C)
775: {
776: Mat_Product *product = C->product;
778: PetscFunctionBegin;
779: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Diagonal;
780: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Diagonal;
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode MatProductSetFromOptions_Diagonal_ConstDiag(Mat C)
785: {
786: Mat_Product *product = C->product;
788: PetscFunctionBegin;
789: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_ConstDiag;
790: else if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_AB_Diagonal_ConstDiag;
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense_AB(Mat C)
795: {
796: PetscFunctionBegin;
797: C->ops->productsymbolic = MatProductSymbolic_ConstDiag_Dense;
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Dense(Mat C)
802: {
803: Mat_Product *product = C->product;
805: PetscFunctionBegin;
806: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_ConstDiag_Dense_AB(C));
807: else if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode MatProductSetFromOptions_ConstDiag_Anytype(Mat C)
812: {
813: Mat_Product *product = C->product;
814: PetscBool Acdiag, Bcdiag;
816: PetscFunctionBegin;
817: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATCONSTANTDIAGONAL, &Acdiag));
818: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, MATCONSTANTDIAGONAL, &Bcdiag));
819: if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_ConstDiag_Anytype;
820: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Anytype_ConstDiag;
821: else if (Acdiag && !Bcdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_ConstDiag_Anytype;
822: else if (Bcdiag && !Acdiag && (product->type == MATPRODUCT_AB)) C->ops->productsymbolic = MatProductSymbolic_AB_Anytype_ConstDiag;
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*MC
827: MATCONSTANTDIAGONAL - "constant-diagonal" - A diagonal matrix type with a uniform value
828: along the diagonal.
830: Level: advanced
832: .seealso: [](ch_matrices), `Mat`, `MatCreateConstantDiagonal()`
833: M*/
834: PETSC_EXTERN PetscErrorCode MatCreate_ConstantDiagonal(Mat A)
835: {
836: Mat_ConstantDiagonal *ctx;
838: PetscFunctionBegin;
839: PetscCall(PetscNew(&ctx));
840: ctx->diag = 0.0;
841: A->data = (void *)ctx;
843: A->assembled = PETSC_TRUE;
844: A->preallocated = PETSC_TRUE;
845: A->structurally_symmetric = PETSC_BOOL3_TRUE;
846: A->structural_symmetry_eternal = PETSC_TRUE;
847: A->symmetric = PETSC_BOOL3_TRUE;
848: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
849: A->symmetry_eternal = PETSC_TRUE;
851: A->ops->mult = MatMult_ConstantDiagonal;
852: A->ops->multadd = MatMultAdd_ConstantDiagonal;
853: A->ops->multtranspose = MatMult_ConstantDiagonal;
854: A->ops->multtransposeadd = MatMultAdd_ConstantDiagonal;
855: A->ops->multhermitiantranspose = MatMultHermitianTranspose_ConstantDiagonal;
856: A->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_ConstantDiagonal;
857: A->ops->solve = MatSolve_ConstantDiagonal;
858: A->ops->solvetranspose = MatSolve_ConstantDiagonal;
859: A->ops->norm = MatNorm_ConstantDiagonal;
860: A->ops->createsubmatrices = MatCreateSubMatrices_ConstantDiagonal;
861: A->ops->duplicate = MatDuplicate_ConstantDiagonal;
862: A->ops->getrow = MatGetRow_ConstantDiagonal;
863: A->ops->restorerow = MatRestoreRow_ConstantDiagonal;
864: A->ops->sor = MatSOR_ConstantDiagonal;
865: A->ops->shift = MatShift_ConstantDiagonal;
866: A->ops->scale = MatScale_ConstantDiagonal;
867: A->ops->getdiagonal = MatGetDiagonal_ConstantDiagonal;
868: A->ops->view = MatView_ConstantDiagonal;
869: A->ops->zeroentries = MatZeroEntries_ConstantDiagonal;
870: A->ops->destroy = MatDestroy_ConstantDiagonal;
871: A->ops->getinfo = MatGetInfo_ConstantDiagonal;
872: A->ops->equal = MatEqual_ConstantDiagonal;
873: A->ops->axpy = MatAXPY_ConstantDiagonal;
874: A->ops->setrandom = MatSetRandom_ConstantDiagonal;
875: A->ops->conjugate = MatConjugate_ConstantDiagonal;
876: A->ops->transpose = MatTranspose_ConstantDiagonal;
878: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCONSTANTDIAGONAL));
879: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConstantDiagonalGetConstant_C", MatConstantDiagonalGetConstant_ConstantDiagonal));
880: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_constantdiagonal_C", MatProductSetFromOptions_ConstDiag_ConstDiag));
881: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_diagonal_C", MatProductSetFromOptions_ConstDiag_Diagonal));
882: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_constantdiagonal_C", MatProductSetFromOptions_Diagonal_ConstDiag));
883: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_seqdense_C", MatProductSetFromOptions_ConstDiag_Dense));
884: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_constantdiagonal_mpidense_C", MatProductSetFromOptions_ConstDiag_Dense));
885: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_ConstDiag_Anytype));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: static PetscErrorCode MatFactorNumeric_ConstantDiagonal(Mat fact, Mat A, const MatFactorInfo *info)
890: {
891: Mat_ConstantDiagonal *actx = (Mat_ConstantDiagonal *)A->data, *fctx = (Mat_ConstantDiagonal *)fact->data;
893: PetscFunctionBegin;
894: if (actx->diag == 0.0) fact->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
895: else fact->factorerrortype = MAT_FACTOR_NOERROR;
896: fctx->diag = 1.0 / actx->diag;
897: fact->ops->solve = MatMult_ConstantDiagonal;
898: PetscFunctionReturn(PETSC_SUCCESS);
899: }
901: static PetscErrorCode MatFactorSymbolic_LU_ConstantDiagonal(Mat fact, Mat A, IS isrow, IS iscol, const MatFactorInfo *info)
902: {
903: PetscFunctionBegin;
904: fact->ops->lufactornumeric = MatFactorNumeric_ConstantDiagonal;
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: static PetscErrorCode MatFactorSymbolic_Cholesky_ConstantDiagonal(Mat fact, Mat A, IS isrow, const MatFactorInfo *info)
909: {
910: PetscFunctionBegin;
911: fact->ops->choleskyfactornumeric = MatFactorNumeric_ConstantDiagonal;
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: PETSC_INTERN PetscErrorCode MatGetFactor_constantdiagonal_petsc(Mat A, MatFactorType ftype, Mat *B)
916: {
917: PetscInt n = A->rmap->n, N = A->rmap->N;
919: PetscFunctionBegin;
920: PetscCall(MatCreateConstantDiagonal(PetscObjectComm((PetscObject)A), n, n, N, N, 0, B));
922: (*B)->factortype = ftype;
923: (*B)->ops->ilufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
924: (*B)->ops->lufactorsymbolic = MatFactorSymbolic_LU_ConstantDiagonal;
925: (*B)->ops->iccfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
926: (*B)->ops->choleskyfactorsymbolic = MatFactorSymbolic_Cholesky_ConstantDiagonal;
928: (*B)->ops->shift = NULL;
929: (*B)->ops->scale = NULL;
930: (*B)->ops->mult = NULL;
931: (*B)->ops->sor = NULL;
932: (*B)->ops->zeroentries = NULL;
934: PetscCall(PetscFree((*B)->solvertype));
935: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*B)->solvertype));
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }