Actual source code: diagonal.c
1: #include <petsc/private/matimpl.h>
2: #include <petsc/private/vecimpl.h>
4: static PetscErrorCode MatDiagonalSetUpDiagonal(Mat A)
5: {
6: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
8: PetscFunctionBegin;
9: if (!ctx->diag_valid) {
10: PetscAssert(ctx->inv_diag_valid, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Neither diagonal nor inverse diagonal is in a valid state");
11: PetscCall(VecCopy(ctx->inv_diag, ctx->diag));
12: PetscCall(VecReciprocal(ctx->diag));
13: ctx->diag_valid = PETSC_TRUE;
14: }
15: PetscFunctionReturn(PETSC_SUCCESS);
16: }
18: static PetscErrorCode MatDiagonalSetUpInverseDiagonal(Mat A)
19: {
20: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
22: PetscFunctionBegin;
23: if (!ctx->inv_diag_valid) {
24: PetscAssert(ctx->diag_valid, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Neither diagonal nor inverse diagonal is in a valid state");
25: PetscCall(VecCopy(ctx->diag, ctx->inv_diag));
26: PetscCall(VecReciprocal(ctx->inv_diag));
27: ctx->inv_diag_valid = PETSC_TRUE;
28: }
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: static PetscErrorCode MatAXPY_Diagonal(Mat Y, PetscScalar a, Mat X, MatStructure str)
33: {
34: Mat_Diagonal *yctx = (Mat_Diagonal *)Y->data;
35: Mat_Diagonal *xctx = (Mat_Diagonal *)X->data;
37: PetscFunctionBegin;
38: PetscCall(MatDiagonalSetUpDiagonal(Y));
39: PetscCall(MatDiagonalSetUpDiagonal(X));
40: PetscCall(VecAXPY(yctx->diag, a, xctx->diag));
41: yctx->inv_diag_valid = PETSC_FALSE;
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: static PetscErrorCode MatGetRow_Diagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
46: {
47: Mat_Diagonal *mat = (Mat_Diagonal *)A->data;
48: PetscInt rstart = A->rmap->rstart, rend = A->rmap->rend;
50: PetscFunctionBegin;
51: PetscCheck(row >= rstart && row < rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only local rows");
52: if (ncols) *ncols = 1;
53: if (cols) {
54: if (!mat->col) PetscCall(PetscMalloc1(1, &mat->col));
55: *mat->col = row;
56: *cols = mat->col;
57: }
58: if (vals) {
59: const PetscScalar *v;
61: if (!mat->val) PetscCall(PetscMalloc1(1, &mat->val));
62: PetscCall(VecGetArrayRead(mat->diag, &v));
63: *mat->val = v[row - rstart];
64: *vals = mat->val;
65: PetscCall(VecRestoreArrayRead(mat->diag, &v));
66: }
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode MatMult_Diagonal(Mat A, Vec x, Vec y)
71: {
72: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
74: PetscFunctionBegin;
75: PetscCall(MatDiagonalSetUpDiagonal(A));
76: PetscCall(VecPointwiseMult(y, ctx->diag, x));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode MatMultAdd_Diagonal(Mat mat, Vec v1, Vec v2, Vec v3)
81: {
82: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
84: PetscFunctionBegin;
85: PetscCall(MatDiagonalSetUpDiagonal(mat));
86: if (v2 != v3) {
87: PetscCall(VecPointwiseMult(v3, ctx->diag, v1));
88: PetscCall(VecAXPY(v3, 1.0, v2));
89: } else {
90: Vec w;
91: PetscCall(VecDuplicate(v3, &w));
92: PetscCall(VecPointwiseMult(w, ctx->diag, v1));
93: PetscCall(VecAXPY(v3, 1.0, w));
94: PetscCall(VecDestroy(&w));
95: }
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode MatNorm_Diagonal(Mat A, NormType type, PetscReal *nrm)
100: {
101: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
103: PetscFunctionBegin;
104: PetscCall(MatDiagonalSetUpDiagonal(A));
105: type = (type == NORM_FROBENIUS) ? NORM_2 : NORM_INFINITY;
106: PetscCall(VecNorm(ctx->diag, type, nrm));
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }
110: static PetscErrorCode MatDuplicate_Diagonal(Mat A, MatDuplicateOption op, Mat *B)
111: {
112: Mat_Diagonal *actx = (Mat_Diagonal *)A->data;
114: PetscFunctionBegin;
115: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
116: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
117: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
118: PetscCall(MatSetType(*B, MATDIAGONAL));
119: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
120: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
121: if (op == MAT_COPY_VALUES) {
122: Mat_Diagonal *bctx = (Mat_Diagonal *)(*B)->data;
124: PetscCall(MatSetUp(A));
125: PetscCall(MatSetUp(*B));
126: PetscCall(VecCopy(actx->diag, bctx->diag));
127: PetscCall(VecCopy(actx->inv_diag, bctx->inv_diag));
128: bctx->diag_valid = actx->diag_valid;
129: bctx->inv_diag_valid = actx->inv_diag_valid;
130: }
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /*@
135: MatDiagonalGetDiagonal - Get the diagonal of a `MATDIAGONAL`
137: Input Parameter:
138: . A - the `MATDIAGONAL`
140: Output Parameter:
141: . diag - the `Vec` that defines the diagonal
143: Level: developer
145: Note:
146: The user must call
147: `MatDiagonalRestoreDiagonal()` before using the matrix again.
149: For a copy of the diagonal values, rather than a reference, use `MatGetDiagonal()`
151: Any changes to the obtained vector immediately change the action of the `Mat`.
152: The matrix can be changed more efficiently by accessing this vector and changing its values, instead of filling a work vector and using `MatDiagonalSet()`
154: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`, `MatGetDiagonal()`
155: @*/
156: PetscErrorCode MatDiagonalGetDiagonal(Mat A, Vec *diag)
157: {
158: PetscFunctionBegin;
160: PetscAssertPointer(diag, 2);
161: *diag = NULL;
162: PetscUseMethod(A, "MatDiagonalGetDiagonal_C", (Mat, Vec *), (A, diag));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: static PetscErrorCode MatDiagonalGetDiagonal_Diagonal(Mat A, Vec *diag)
167: {
168: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
170: PetscFunctionBegin;
171: PetscCall(MatSetUp(A));
172: PetscCall(MatDiagonalSetUpDiagonal(A));
173: *diag = ctx->diag;
174: PetscCall(PetscObjectStateGet((PetscObject)*diag, &ctx->diag_state));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /*@
179: MatDiagonalRestoreDiagonal - Restore the diagonal of a `MATDIAGONAL`
181: Input Parameters:
182: + A - the `MATDIAGONAL`
183: - diag - the `Vec` obtained from `MatDiagonalGetDiagonal()`
185: Level: developer
187: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetDiagonal()`
188: @*/
189: PetscErrorCode MatDiagonalRestoreDiagonal(Mat A, Vec *diag)
190: {
191: PetscFunctionBegin;
193: PetscAssertPointer(diag, 2);
194: PetscUseMethod(A, "MatDiagonalRestoreDiagonal_C", (Mat, Vec *), (A, diag));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode MatDiagonalRestoreDiagonal_Diagonal(Mat A, Vec *diag)
199: {
200: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
201: PetscObjectState diag_state;
203: PetscFunctionBegin;
204: PetscCheck(ctx->diag == *diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
205: ctx->diag_valid = PETSC_TRUE;
206: PetscCall(PetscObjectStateGet((PetscObject)*diag, &diag_state));
207: if (ctx->diag_state != diag_state) {
208: PetscCall(PetscObjectStateIncrease((PetscObject)A));
209: ctx->inv_diag_valid = PETSC_FALSE;
210: }
211: *diag = NULL;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: MatDiagonalGetInverseDiagonal - Get the inverse diagonal of a `MATDIAGONAL`
218: Input Parameter:
219: . A - the `MATDIAGONAL`
221: Output Parameter:
222: . inv_diag - the `Vec` that defines the inverse diagonal
224: Level: developer
226: Note:
227: The user must call
228: `MatDiagonalRestoreInverseDiagonal()` before using the matrix again.
230: If a matrix is created only to call `MatSolve()` (which happens for `MATLMVMDIAGBROYDEN`),
231: using `MatDiagonalGetInverseDiagonal()` and `MatDiagonalRestoreInverseDiagonal()` avoids copies
232: and avoids any call to `VecReciprocal()`.
234: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MATLMVMBROYDEN`, `MatSolve()`
235: @*/
236: PetscErrorCode MatDiagonalGetInverseDiagonal(Mat A, Vec *inv_diag)
237: {
238: PetscFunctionBegin;
240: PetscAssertPointer(inv_diag, 2);
241: *inv_diag = NULL;
242: PetscUseMethod(A, "MatDiagonalGetInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: static PetscErrorCode MatDiagonalGetInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
247: {
248: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
250: PetscFunctionBegin;
251: PetscCall(MatSetUp(A));
252: PetscCall(MatDiagonalSetUpInverseDiagonal(A));
253: *inv_diag = ctx->inv_diag;
254: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &ctx->inv_diag_state));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@
259: MatDiagonalRestoreInverseDiagonal - Restore the inverse diagonal of a `MATDIAGONAL`
261: Input Parameters:
262: + A - the `MATDIAGONAL`
263: - inv_diag - the `Vec` obtained from `MatDiagonalGetInverseDiagonal()`
265: Level: developer
267: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetInverseDiagonal()`
268: @*/
269: PetscErrorCode MatDiagonalRestoreInverseDiagonal(Mat A, Vec *inv_diag)
270: {
271: PetscFunctionBegin;
273: PetscAssertPointer(inv_diag, 2);
274: PetscUseMethod(A, "MatDiagonalRestoreInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode MatDiagonalRestoreInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
279: {
280: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
281: PetscObjectState inv_diag_state;
283: PetscFunctionBegin;
284: PetscCheck(ctx->inv_diag == *inv_diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
285: ctx->inv_diag_valid = PETSC_TRUE;
286: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &inv_diag_state));
287: if (ctx->inv_diag_state != inv_diag_state) {
288: PetscCall(PetscObjectStateIncrease((PetscObject)A));
289: ctx->diag_valid = PETSC_FALSE;
290: }
291: *inv_diag = NULL;
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: static PetscErrorCode MatPermute_Diagonal(Mat A, IS rowp, IS colp, Mat *B)
296: {
297: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
298: Vec v;
300: PetscFunctionBegin;
301: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
302: PetscCall(VecDuplicate(ctx->diag, &v));
303: PetscCall(VecCopy(ctx->diag, v));
304: PetscCall(VecPermute(v, rowp, PETSC_FALSE));
305: PetscCall(MatCreateDiagonal(v, B));
306: PetscCall(VecDestroy(&v));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode MatSetRandom_Diagonal(Mat A, PetscRandom rand)
311: {
312: Vec d;
314: PetscFunctionBegin;
315: PetscCall(MatDiagonalGetDiagonal(A, &d));
316: PetscCall(VecSetRandom(d, rand));
317: PetscCall(MatDiagonalRestoreDiagonal(A, &d));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode MatDestroy_Diagonal(Mat mat)
322: {
323: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
325: PetscFunctionBegin;
326: PetscCall(VecDestroy(&ctx->diag));
327: PetscCall(VecDestroy(&ctx->inv_diag));
328: PetscCall(PetscFree(ctx->col));
329: PetscCall(PetscFree(ctx->val));
330: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetDiagonal_C", NULL));
331: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreDiagonal_C", NULL));
332: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetInverseDiagonal_C", NULL));
333: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreInverseDiagonal_C", NULL));
334: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_seqdense_C", NULL));
335: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_mpidense_C", NULL));
336: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_diagonal_C", NULL));
337: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_anytype_C", NULL));
338: PetscCall(PetscFree(mat->data));
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode MatView_Diagonal(Mat J, PetscViewer viewer)
343: {
344: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
345: PetscBool isascii;
347: PetscFunctionBegin;
348: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
349: if (isascii) {
350: PetscViewerFormat format;
352: PetscCall(MatDiagonalSetUpDiagonal(J));
353: PetscCall(PetscViewerGetFormat(viewer, &format));
354: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
355: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ctx->diag, viewer));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
358: PetscCall(VecView(ctx->diag, viewer));
359: }
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: static PetscErrorCode MatGetDiagonal_Diagonal(Mat J, Vec x)
364: {
365: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
367: PetscFunctionBegin;
368: PetscCall(MatDiagonalSetUpDiagonal(J));
369: PetscCall(VecCopy(ctx->diag, x));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode MatADot_Diagonal_Local(Mat A, Vec x, Vec y, PetscScalar *val)
374: {
375: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
376: PetscInt n = x->map->n;
377: const PetscScalar *ya, *xa, *wa;
378: PetscScalar sum = 0;
380: PetscFunctionBegin;
381: PetscCheckSameTypeAndComm(x, 2, ctx->diag, 1);
382: PetscCheckSameTypeAndComm(y, 3, ctx->diag, 1);
383: PetscCall(VecGetArrayRead(x, &xa));
384: PetscCall(VecGetArrayRead(y, &ya));
385: PetscCall(VecGetArrayRead(ctx->diag, &wa));
386: for (PetscInt i = 0; i < n; i++) {
387: sum += PetscConj(ya[i]) * wa[i] * xa[i];
388: }
389: if (n > 0) PetscCall(PetscLogFlops(3.0 * n));
390: PetscCall(VecRestoreArrayRead(x, &xa));
391: PetscCall(VecRestoreArrayRead(y, &ya));
392: PetscCall(VecRestoreArrayRead(ctx->diag, &wa));
393: *val = sum;
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode MatADot_Diagonal_MPI(Mat A, Vec x, Vec y, PetscScalar *val)
398: {
399: PetscFunctionBegin;
400: PetscCall(MatDiagonalSetUpDiagonal(A));
401: PetscUseTypeMethod(A, adot_local, x, y, val);
402: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, val, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)A)));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode MatANorm_Diagonal_Local(Mat A, Vec x, PetscReal *val)
407: {
408: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
409: PetscInt n = x->map->n;
410: const PetscScalar *xa, *wa;
411: PetscScalar sum = 0;
413: PetscFunctionBegin;
414: PetscCheckSameTypeAndComm(x, 2, ctx->diag, 1);
415: PetscCall(VecGetArrayRead(x, &xa));
416: PetscCall(VecGetArrayRead(ctx->diag, &wa));
417: for (PetscInt i = 0; i < n; i++) {
418: sum += PetscConj(xa[i]) * wa[i] * xa[i];
419: }
420: if (n > 0) PetscCall(PetscLogFlops(3.0 * n));
421: PetscCall(VecRestoreArrayRead(x, &xa));
422: PetscCall(VecRestoreArrayRead(ctx->diag, &wa));
423: PetscCheck(PetscAbsReal(PetscImaginaryPart(sum)) < 100 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix argument is not Hermitian (diagonal has nonzero imaginary parts)");
424: *val = PetscRealPart(sum);
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode MatANorm_Diagonal_MPI(Mat A, Vec x, PetscReal *val)
429: {
430: PetscFunctionBegin;
431: PetscCall(MatDiagonalSetUpDiagonal(A));
432: PetscUseTypeMethod(A, anorm_local, x, val);
433: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, val, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)A)));
434: *val = PetscSqrtReal(*val);
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode MatANorm_Diagonal_Seq(Mat A, Vec x, PetscReal *val)
439: {
440: PetscFunctionBegin;
441: PetscCall(MatDiagonalSetUpDiagonal(A));
442: PetscUseTypeMethod(A, anorm_local, x, val);
443: PetscCheck(*val >= 0.0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix argument is not positive definite (diagonal has negative entries)");
444: *val = PetscSqrtReal(*val);
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /*@
449: MatDiagonalSetDiagonal - Sets the diagonal for a `MATDIAGONAL`
451: Collective
453: Input Parameter:
454: + J - the `MATDIAGONAL` matrix
455: - diag - the vector for the diagonal
457: Level: advanced
459: Notes:
460: The input vector `diag` will be referenced internally: any changes to `diag`
461: will affect the matrix `J`.
463: This routine can only be called once for the given `J` matrix.
465: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatCreateDiagonal()`, `MATDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatSolve()`,
466: `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`,
467: `MATCONSTANTDIAGONAL`
468: @*/
469: PetscErrorCode MatDiagonalSetDiagonal(Mat J, Vec diag)
470: {
471: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
472: PetscMPIInt comm_size;
473: PetscLayout map;
474: VecType type;
475: PetscBool iskokkos, ismpi, iscuda, iship;
477: PetscFunctionBegin;
480: PetscCheckSameComm(J, 1, diag, 2);
481: PetscCheck(!ctx->diag, PetscObjectComm((PetscObject)J), PETSC_ERR_ARG_WRONGSTATE, "MATDIAGONAL already has a diagonal");
482: ctx->diag = diag;
483: PetscCall(PetscObjectReference((PetscObject)ctx->diag));
484: PetscCall(VecDuplicate(ctx->diag, &ctx->inv_diag));
485: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)J), &comm_size));
486: PetscCall(VecGetLayout(ctx->diag, &map));
487: PetscCall(MatSetLayouts(J, map, map));
488: PetscCall(PetscFree(J->defaultvectype));
489: PetscCall(VecGetType(ctx->diag, &type));
490: PetscCall(PetscStrallocpy(type, &J->defaultvectype));
492: PetscCall(PetscObjectTypeCompareAny((PetscObject)ctx->diag, &iscuda, VECCUDA, VECMPICUDA, VECSEQCUDA, ""));
493: PetscCall(PetscObjectTypeCompareAny((PetscObject)ctx->diag, &iship, VECHIP, VECMPIHIP, VECSEQHIP, ""));
494: PetscCall(PetscObjectTypeCompareAny((PetscObject)ctx->diag, &iskokkos, VECKOKKOS, VECMPIKOKKOS, VECSEQKOKKOS, ""));
495: PetscCall(PetscObjectBaseTypeCompare((PetscObject)ctx->diag, VECMPI, &ismpi));
496: ismpi = ismpi || comm_size > 1;
497: J->ops->adot_local = MatADot_Diagonal_Local;
498: J->ops->adot = MatADot_Diagonal_Local;
499: J->ops->anorm_local = MatANorm_Diagonal_Local;
500: if (iskokkos) {
501: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)ctx->diag), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
502: #if PetscDefined(HAVE_KOKKOS_KERNELS)
503: J->ops->adot_local = MatADot_Diagonal_SeqKokkos;
504: J->ops->adot = MatADot_Diagonal_SeqKokkos;
505: J->ops->anorm_local = MatANormSq_Diagonal_SeqKokkos;
506: #endif
507: }
508: if (iscuda) {
509: PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)ctx->diag), PETSC_ERR_SUP, "Reconfigure using CUDA support");
510: #if PetscDefined(HAVE_CUDA)
511: J->ops->adot_local = MatADot_Diagonal_SeqCUDA;
512: J->ops->adot = MatADot_Diagonal_SeqCUDA;
513: J->ops->anorm_local = MatANormSq_Diagonal_SeqCUDA;
514: #endif
515: }
516: if (iship) {
517: PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)ctx->diag), PETSC_ERR_SUP, "Reconfigure using HIP support");
518: #if PetscDefined(HAVE_HIP)
519: J->ops->adot_local = MatADot_Diagonal_SeqHIP;
520: J->ops->adot = MatADot_Diagonal_SeqHIP;
521: J->ops->anorm_local = MatANormSq_Diagonal_SeqHIP;
522: #endif
523: }
524: if (ismpi) {
525: J->ops->adot = MatADot_Diagonal_MPI;
526: J->ops->anorm = MatANorm_Diagonal_MPI;
527: } else {
528: J->ops->anorm = MatANorm_Diagonal_Seq;
529: }
530: ctx->col = NULL;
531: ctx->val = NULL;
532: ctx->diag_valid = PETSC_TRUE;
533: J->assembled = PETSC_TRUE;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: static PetscErrorCode MatDiagonalSet_Diagonal(Mat J, Vec D, InsertMode is)
538: {
539: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
541: PetscFunctionBegin;
542: switch (is) {
543: case ADD_VALUES:
544: case ADD_ALL_VALUES:
545: case ADD_BC_VALUES:
546: PetscCall(MatDiagonalSetUpDiagonal(J));
547: PetscCall(VecAXPY(ctx->diag, 1.0, D));
548: ctx->inv_diag_valid = PETSC_FALSE;
549: break;
550: case INSERT_VALUES:
551: case INSERT_BC_VALUES:
552: case INSERT_ALL_VALUES:
553: PetscCall(VecCopy(D, ctx->diag));
554: ctx->diag_valid = PETSC_TRUE;
555: ctx->inv_diag_valid = PETSC_FALSE;
556: break;
557: case MAX_VALUES:
558: PetscCall(MatDiagonalSetUpDiagonal(J));
559: PetscCall(VecPointwiseMax(ctx->diag, D, ctx->diag));
560: ctx->inv_diag_valid = PETSC_FALSE;
561: break;
562: case MIN_VALUES:
563: PetscCall(MatDiagonalSetUpDiagonal(J));
564: PetscCall(VecPointwiseMin(ctx->diag, D, ctx->diag));
565: ctx->inv_diag_valid = PETSC_FALSE;
566: break;
567: case NOT_SET_VALUES:
568: break;
569: }
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: static PetscErrorCode MatShift_Diagonal(Mat Y, PetscScalar a)
574: {
575: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
577: PetscFunctionBegin;
578: PetscCall(MatDiagonalSetUpDiagonal(Y));
579: PetscCall(VecShift(ctx->diag, a));
580: ctx->inv_diag_valid = PETSC_FALSE;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode MatScale_Diagonal(Mat Y, PetscScalar a)
585: {
586: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
588: PetscFunctionBegin;
589: PetscCall(MatDiagonalSetUpDiagonal(Y));
590: PetscCall(VecScale(ctx->diag, a));
591: ctx->inv_diag_valid = PETSC_FALSE;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: static PetscErrorCode MatDiagonalScale_Diagonal(Mat Y, Vec l, Vec r)
596: {
597: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
599: PetscFunctionBegin;
600: PetscCall(MatDiagonalSetUpDiagonal(Y));
601: if (l) {
602: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, l));
603: ctx->inv_diag_valid = PETSC_FALSE;
604: }
605: if (r) {
606: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, r));
607: ctx->inv_diag_valid = PETSC_FALSE;
608: }
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode MatConjugate_Diagonal(Mat Y)
613: {
614: PetscFunctionBegin;
615: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
617: if (ctx->diag_valid) {
618: PetscCall(VecConjugate(ctx->diag));
619: PetscCall(PetscObjectStateGet((PetscObject)ctx->diag, &ctx->diag_state));
620: }
621: if (ctx->inv_diag_valid) {
622: PetscCall(VecConjugate(ctx->inv_diag));
623: PetscCall(PetscObjectStateGet((PetscObject)ctx->inv_diag, &ctx->inv_diag_state));
624: }
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: static PetscErrorCode MatTranspose_Diagonal(Mat A, MatReuse reuse, Mat *matout)
629: {
630: PetscFunctionBegin;
631: if (reuse == MAT_INPLACE_MATRIX) {
632: PetscLayout tmplayout = A->rmap;
634: A->rmap = A->cmap;
635: A->cmap = tmplayout;
636: } else {
637: Vec diag, newdiag;
638: if (reuse == MAT_INITIAL_MATRIX) {
639: PetscCall(MatDiagonalGetDiagonal(A, &diag));
640: PetscCall(VecDuplicate(diag, &newdiag));
641: PetscCall(VecCopy(diag, newdiag));
642: PetscCall(MatDiagonalRestoreDiagonal(A, &diag));
643: PetscCall(MatCreateDiagonal(newdiag, matout));
644: PetscCall(VecDestroy(&newdiag));
645: } else {
646: PetscCall(MatDiagonalGetDiagonal(A, &diag));
647: PetscCall(MatDiagonalGetDiagonal(*matout, &newdiag));
648: PetscCall(VecCopy(diag, newdiag));
649: PetscCall(MatDiagonalRestoreDiagonal(*matout, &newdiag));
650: PetscCall(MatDiagonalRestoreDiagonal(A, &diag));
651: }
652: }
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode MatSetUp_Diagonal(Mat A)
657: {
658: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
660: PetscFunctionBegin;
661: if (!ctx->diag) {
662: Vec diag;
664: PetscCall(PetscLayoutSetUp(A->cmap));
665: PetscCall(PetscLayoutSetUp(A->rmap));
666: PetscCall(MatCreateVecs(A, &diag, NULL));
667: PetscCall(MatDiagonalSetDiagonal(A, diag));
668: PetscCall(VecDestroy(&diag));
669: }
670: A->assembled = PETSC_TRUE;
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: static PetscErrorCode MatZeroEntries_Diagonal(Mat Y)
675: {
676: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
678: PetscFunctionBegin;
679: PetscCall(VecZeroEntries(ctx->diag));
680: ctx->diag_valid = PETSC_TRUE;
681: ctx->inv_diag_valid = PETSC_FALSE;
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode MatSolve_Diagonal(Mat matin, Vec b, Vec x)
686: {
687: Mat_Diagonal *ctx = (Mat_Diagonal *)matin->data;
689: PetscFunctionBegin;
690: PetscCall(MatDiagonalSetUpInverseDiagonal(matin));
691: PetscCall(VecPointwiseMult(x, b, ctx->inv_diag));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: static PetscErrorCode MatGetInfo_Diagonal(Mat A, MatInfoType flag, MatInfo *info)
696: {
697: PetscFunctionBegin;
698: info->block_size = 1.0;
699: info->nz_allocated = A->cmap->N;
700: info->nz_used = A->cmap->N;
701: info->nz_unneeded = 0.0;
702: info->assemblies = A->num_ass;
703: info->mallocs = 0.0;
704: info->memory = 0;
705: info->fill_ratio_given = 0;
706: info->fill_ratio_needed = 0;
707: info->factor_mallocs = 0;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: /*@
712: MatCreateDiagonal - Creates a matrix defined by a given vector along its diagonal.
714: Collective
716: Input Parameter:
717: . diag - vector for the diagonal
719: Output Parameter:
720: . J - the diagonal matrix
722: Level: advanced
724: Notes:
725: Only supports square matrices with the same number of local rows and columns.
727: The input vector `diag` will be referenced internally: any changes to `diag`
728: will affect the matrix `J`.
730: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatSolve()`,
731: `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`,
732: `MATCONSTANTDIAGONAL`
733: @*/
734: PetscErrorCode MatCreateDiagonal(Vec diag, Mat *J)
735: {
736: PetscFunctionBegin;
738: PetscCall(MatCreate(PetscObjectComm((PetscObject)diag), J));
739: PetscCall(MatSetType(*J, MATDIAGONAL));
740: PetscCall(MatDiagonalSetDiagonal(*J, diag));
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode MatProductNumeric_Diagonal_Dense(Mat C)
745: {
746: Mat A, B;
747: Mat_Diagonal *a;
748: PetscScalar *c;
749: const PetscScalar *b, *alpha;
750: PetscInt ldb, ldc;
752: PetscFunctionBegin;
753: MatCheckProduct(C, 1);
754: A = C->product->A;
755: B = C->product->B;
756: PetscCall(MatDiagonalSetUpDiagonal(A));
757: a = (Mat_Diagonal *)A->data;
758: PetscCall(VecGetArrayRead(a->diag, &alpha));
759: PetscCall(MatDenseGetLDA(B, &ldb));
760: PetscCall(MatDenseGetLDA(C, &ldc));
761: PetscCall(MatDenseGetArrayRead(B, &b));
762: PetscCall(MatDenseGetArrayWrite(C, &c));
763: for (PetscInt j = 0; j < B->cmap->N; j++)
764: for (PetscInt i = 0; i < B->rmap->n; i++) c[i + j * ldc] = alpha[i] * b[i + j * ldb];
765: PetscCall(MatDenseRestoreArrayWrite(C, &c));
766: PetscCall(MatDenseRestoreArrayRead(B, &b));
767: PetscCall(VecRestoreArrayRead(a->diag, &alpha));
768: PetscCall(PetscLogFlops(B->cmap->N * B->rmap->n));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: static PetscErrorCode MatProductSymbolic_Diagonal_Dense(Mat C)
773: {
774: Mat A, B;
775: PetscInt n, N, m, M;
777: PetscFunctionBegin;
778: MatCheckProduct(C, 1);
779: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
780: A = C->product->A;
781: B = C->product->B;
782: PetscCall(MatGetLocalSize(C, &m, &n));
783: PetscCall(MatGetSize(C, &M, &N));
784: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
785: PetscCall(MatGetLocalSize(B, NULL, &n));
786: PetscCall(MatGetSize(B, NULL, &N));
787: PetscCall(MatGetLocalSize(A, &m, NULL));
788: PetscCall(MatGetSize(A, &M, NULL));
789: PetscCall(MatSetSizes(C, m, n, M, N));
790: }
791: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
792: PetscCall(MatSetUp(C));
793: C->ops->productnumeric = MatProductNumeric_Diagonal_Dense;
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /* PtAP for diagonal * diagonal: C_ii = d_A_i * d_P_i^2 */
798: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_Diagonal(Mat C)
799: {
800: Mat A = C->product->A, P = C->product->B;
801: Mat_Diagonal *a = (Mat_Diagonal *)A->data, *p = (Mat_Diagonal *)P->data, *c = (Mat_Diagonal *)C->data;
803: PetscFunctionBegin;
804: MatCheckProduct(C, 1);
805: PetscCall(MatDiagonalSetUpDiagonal(A));
806: PetscCall(MatDiagonalSetUpDiagonal(P));
807: PetscCall(VecPointwiseMult(c->diag, a->diag, p->diag));
808: PetscCall(VecPointwiseMult(c->diag, c->diag, p->diag));
809: c->diag_valid = PETSC_TRUE;
810: c->inv_diag_valid = PETSC_FALSE;
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_Diagonal(Mat C)
815: {
816: Mat P = C->product->B;
817: Mat_Diagonal *p = (Mat_Diagonal *)P->data;
818: Vec cdiag;
820: PetscFunctionBegin;
821: MatCheckProduct(C, 1);
822: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
823: PetscCall(MatSetSizes(C, P->cmap->n, P->cmap->n, P->cmap->N, P->cmap->N));
824: PetscCall(MatSetType(C, MATDIAGONAL));
825: /* Duplicate P's diagonal vec so C inherits the correct vec type (e.g., Kokkos, CUDA, HIP) */
826: PetscCall(MatDiagonalSetUpDiagonal(P));
827: PetscCall(VecDuplicate(p->diag, &cdiag));
828: PetscCall(MatDiagonalSetDiagonal(C, cdiag));
829: PetscCall(VecDestroy(&cdiag));
830: C->assembled = PETSC_TRUE;
831: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_Diagonal;
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /* PtAP for any * diagonal: C = D * A * D (bilateral diagonal scaling) */
836: static PetscErrorCode MatProductNumeric_PtAP_Any_Diagonal(Mat C)
837: {
838: Mat A = C->product->A, P = C->product->B;
839: Mat_Diagonal *p = (Mat_Diagonal *)P->data;
841: PetscFunctionBegin;
842: MatCheckProduct(C, 1);
843: PetscCall(MatDiagonalSetUpDiagonal(P));
844: PetscCall(MatCopy(A, C, SAME_NONZERO_PATTERN));
845: PetscCall(MatDiagonalScale(C, p->diag, p->diag));
846: PetscFunctionReturn(PETSC_SUCCESS);
847: }
849: static PetscErrorCode MatProductSymbolic_PtAP_Any_Diagonal(Mat C)
850: {
851: Mat A = C->product->A;
852: Mat_Product *product = C->product;
853: Mat Cwork;
855: PetscFunctionBegin;
856: MatCheckProduct(C, 1);
857: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
858: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &Cwork));
859: C->product = NULL;
860: PetscCall(MatHeaderReplace(C, &Cwork));
861: C->product = product;
862: C->ops->productnumeric = MatProductNumeric_PtAP_Any_Diagonal;
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: /* PtAP for diagonal A and any non-diagonal P: C = P^T * A * P
867: Computed as P^T * AP where AP = A*P is a diagonal-scaled copy of P.
868: The Unsafe decomposition is not used because it requires the old-style
869: transposematmultnumeric pointer, which GPU types do not set. */
870: typedef struct {
871: Mat AP; /* AP = A*P (scaled copy of P) */
872: Mat PtAP; /* P^T * AP result via MatProduct AtB */
873: PetscObjectState pstate; /* P's state when inner product was last built */
874: } MatProductCtx_PtAP_DiagAny;
876: static PetscErrorCode MatProductCtxDestroy_PtAP_DiagAny(PetscCtxRt data)
877: {
878: MatProductCtx_PtAP_DiagAny *ctx = *(MatProductCtx_PtAP_DiagAny **)data;
880: PetscFunctionBegin;
881: PetscCall(MatDestroy(&ctx->AP));
882: PetscCall(MatDestroy(&ctx->PtAP));
883: PetscCall(PetscFree(ctx));
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: static PetscErrorCode MatProductNumeric_PtAP_Diagonal_Any(Mat C)
888: {
889: Mat_Product *product = C->product;
890: Mat A = product->A, P = product->B;
891: MatProductCtx_PtAP_DiagAny *ctx = (MatProductCtx_PtAP_DiagAny *)product->data;
892: Mat_Diagonal *a = (Mat_Diagonal *)A->data;
893: PetscObjectState pstate;
895: PetscFunctionBegin;
896: MatCheckProduct(C, 1);
897: /* Update AP = A * P */
898: PetscCall(MatDiagonalSetUpDiagonal(A));
899: PetscCall(MatCopy(P, ctx->AP, SAME_NONZERO_PATTERN));
900: PetscCall(MatDiagonalScale(ctx->AP, a->diag, NULL));
901: /* Rebuild inner product if P has changed (some backends, e.g. CUSPARSE,
902: do not properly recompute the transpose of the left operand on reuse) */
903: PetscCall(PetscObjectStateGet((PetscObject)P, &pstate));
904: if (pstate != ctx->pstate) {
905: PetscCall(MatDestroy(&ctx->PtAP));
906: PetscCall(MatProductCreate(P, ctx->AP, NULL, &ctx->PtAP));
907: PetscCall(MatProductSetType(ctx->PtAP, MATPRODUCT_AtB));
908: PetscCall(MatProductSetFill(ctx->PtAP, product->fill));
909: PetscCall(MatProductSetFromOptions(ctx->PtAP));
910: PetscCall(MatProductSymbolic(ctx->PtAP));
911: ctx->pstate = pstate;
912: }
913: /* Recompute P^T * AP */
914: PetscCall(MatProductNumeric(ctx->PtAP));
915: PetscCall(MatCopy(ctx->PtAP, C, SAME_NONZERO_PATTERN));
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: static PetscErrorCode MatProductSymbolic_PtAP_Diagonal_Any(Mat C)
920: {
921: Mat_Product *product = C->product;
922: Mat P = product->B;
923: MatProductCtx_PtAP_DiagAny *ctx;
924: Mat Cwork;
926: PetscFunctionBegin;
927: MatCheckProduct(C, 1);
928: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
929: PetscCall(PetscNew(&ctx));
931: /* AP = A * P (structure only; numeric phase fills values) */
932: PetscCall(MatDuplicate(P, MAT_DO_NOT_COPY_VALUES, &ctx->AP));
934: /* PtAP = P^T * AP (symbolic only: D*P has same sparsity as P) */
935: PetscCall(MatProductCreate(P, ctx->AP, NULL, &ctx->PtAP));
936: PetscCall(MatProductSetType(ctx->PtAP, MATPRODUCT_AtB));
937: PetscCall(MatProductSetFill(ctx->PtAP, product->fill));
938: PetscCall(MatProductSetFromOptions(ctx->PtAP));
939: PetscCall(MatProductSymbolic(ctx->PtAP));
941: /* Record P's state so numeric phase can detect changes */
942: PetscCall(PetscObjectStateGet((PetscObject)P, &ctx->pstate));
944: /* Set up C with the same structure as PtAP */
945: PetscCall(MatDuplicate(ctx->PtAP, MAT_DO_NOT_COPY_VALUES, &Cwork));
946: C->product = NULL;
947: PetscCall(MatHeaderReplace(C, &Cwork));
948: C->product = product;
949: C->assembled = PETSC_TRUE;
950: product->data = ctx;
951: product->destroy = MatProductCtxDestroy_PtAP_DiagAny;
952: C->ops->productnumeric = MatProductNumeric_PtAP_Diagonal_Any;
953: PetscFunctionReturn(PETSC_SUCCESS);
954: }
956: static PetscErrorCode MatProductSetFromOptions_Diagonal_Diagonal(Mat C)
957: {
958: Mat_Product *product = C->product;
960: PetscFunctionBegin;
961: if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_Diagonal;
962: PetscFunctionReturn(PETSC_SUCCESS);
963: }
965: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense_AB(Mat C)
966: {
967: PetscFunctionBegin;
968: C->ops->productsymbolic = MatProductSymbolic_Diagonal_Dense;
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense(Mat C)
973: {
974: Mat_Product *product = C->product;
976: PetscFunctionBegin;
977: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_Diagonal_Dense_AB(C));
978: else if (product->type == MATPRODUCT_PtAP) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_Any;
979: PetscFunctionReturn(PETSC_SUCCESS);
980: }
982: /*
983: Fallback dispatcher for matrix type combinations that have no specific registration
984: (e.g., GPU types like seqaijkokkos). Queried via the "anytype" mechanism in
985: MatProductSetFromOptions_Private() after specific type queries fail.
986: */
987: static PetscErrorCode MatProductSetFromOptions_Diagonal_Anytype(Mat C)
988: {
989: Mat_Product *product = C->product;
990: PetscBool Adiag, Bdiag;
992: PetscFunctionBegin;
993: PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATDIAGONAL, &Adiag));
994: PetscCall(PetscObjectTypeCompare((PetscObject)product->B, MATDIAGONAL, &Bdiag));
995: if (Adiag && !Bdiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Diagonal_Any;
996: else if (Bdiag && !Adiag && (product->type == MATPRODUCT_PtAP)) C->ops->productsymbolic = MatProductSymbolic_PtAP_Any_Diagonal;
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*MC
1001: MATDIAGONAL - MATDIAGONAL = "diagonal" - A diagonal matrix type with the diagonal implemented as a `Vec`. Useful for
1002: cases where `VecPointwiseMult()` or `VecPointwiseDivide()` should be thought of as the actions of a linear operator.
1004: Options Database Key:
1005: . -mat_vec_type type - set the `VecType` of the vector defining the diagonal
1007: Level: advanced
1009: Note:
1010: The first use case is simply to call `MatCreateDiagonal()' to provide the vector containing the diagonal matrix. The input vector will be
1011: referenced internally by the matrix: any changes to it will affect the matrix. Similar changes to the matrix will affect the vector.
1013: For the second use case call `MatSetType()` with a type of `MATDIAGONAL` followed by a call to `MatDiagonalSetDiagonal()`. The input vector will be
1014: referenced internally by the matrix: any changes to it will affect the matrix. Similar changes to the matrix will affect the vector.
1016: For the third use case call `MatSetType()` with a type of `MATDIAGONAL` followed by calls to `MatSetSizes()` and `MatDiagonalSet()`
1017: (in this case the diagonal vector will not be referenced internally by the matrix, its values will be copied) or some other
1018: operation to provide the matrix entries. One can control the `VecType` of the diagonal by calling `MatSetVecType()` or using `-mat_vec_type type`.
1020: .seealso: [](ch_matrices), `Mat`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`,
1021: `MatDiagonalGetInverseDiagonal()`, `MatDiagonalSet()`, `MatDiagonalSetDiagonal()`, `MatSetVecType()`
1022: M*/
1023: PETSC_INTERN PetscErrorCode MatCreate_Diagonal(Mat A)
1024: {
1025: Mat_Diagonal *ctx;
1027: PetscFunctionBegin;
1028: PetscCall(PetscNew(&ctx));
1029: A->data = (void *)ctx;
1031: A->structurally_symmetric = PETSC_BOOL3_TRUE;
1032: A->structural_symmetry_eternal = PETSC_TRUE;
1033: A->symmetry_eternal = PETSC_TRUE;
1034: A->symmetric = PETSC_BOOL3_TRUE;
1035: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
1037: A->ops->getrow = MatGetRow_Diagonal;
1038: A->ops->mult = MatMult_Diagonal;
1039: A->ops->multadd = MatMultAdd_Diagonal;
1040: A->ops->multtranspose = MatMult_Diagonal;
1041: A->ops->multtransposeadd = MatMultAdd_Diagonal;
1042: A->ops->norm = MatNorm_Diagonal;
1043: A->ops->duplicate = MatDuplicate_Diagonal;
1044: A->ops->solve = MatSolve_Diagonal;
1045: A->ops->solvetranspose = MatSolve_Diagonal;
1046: A->ops->shift = MatShift_Diagonal;
1047: A->ops->scale = MatScale_Diagonal;
1048: A->ops->diagonalscale = MatDiagonalScale_Diagonal;
1049: A->ops->getdiagonal = MatGetDiagonal_Diagonal;
1050: A->ops->diagonalset = MatDiagonalSet_Diagonal;
1051: A->ops->view = MatView_Diagonal;
1052: A->ops->zeroentries = MatZeroEntries_Diagonal;
1053: A->ops->destroy = MatDestroy_Diagonal;
1054: A->ops->getinfo = MatGetInfo_Diagonal;
1055: A->ops->axpy = MatAXPY_Diagonal;
1056: A->ops->setup = MatSetUp_Diagonal;
1057: A->ops->permute = MatPermute_Diagonal;
1058: A->ops->setrandom = MatSetRandom_Diagonal;
1059: A->ops->conjugate = MatConjugate_Diagonal;
1060: A->ops->transpose = MatTranspose_Diagonal;
1062: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetDiagonal_C", MatDiagonalGetDiagonal_Diagonal));
1063: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreDiagonal_C", MatDiagonalRestoreDiagonal_Diagonal));
1064: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetInverseDiagonal_C", MatDiagonalGetInverseDiagonal_Diagonal));
1065: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreInverseDiagonal_C", MatDiagonalRestoreInverseDiagonal_Diagonal));
1066: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_seqdense_C", MatProductSetFromOptions_Diagonal_Dense));
1067: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_mpidense_C", MatProductSetFromOptions_Diagonal_Dense));
1068: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_diagonal_C", MatProductSetFromOptions_Diagonal_Diagonal));
1069: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Diagonal_Anytype));
1070: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATDIAGONAL));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }