Actual source code: diagonal.c
1: #include <petsc/private/matimpl.h>
3: typedef struct {
4: Vec diag;
5: PetscBool diag_valid;
6: Vec inv_diag;
7: PetscBool inv_diag_valid;
8: PetscObjectState diag_state, inv_diag_state;
9: PetscInt *col;
10: PetscScalar *val;
11: } Mat_Diagonal;
13: static PetscErrorCode MatDiagonalSetUpDiagonal(Mat A)
14: {
15: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
17: PetscFunctionBegin;
18: if (!ctx->diag_valid) {
19: PetscAssert(ctx->inv_diag_valid, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Neither diagonal nor inverse diagonal is in a valid state");
20: PetscCall(VecCopy(ctx->inv_diag, ctx->diag));
21: PetscCall(VecReciprocal(ctx->diag));
22: ctx->diag_valid = PETSC_TRUE;
23: }
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode MatDiagonalSetUpInverseDiagonal(Mat A)
28: {
29: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
31: PetscFunctionBegin;
32: if (!ctx->inv_diag_valid) {
33: PetscAssert(ctx->diag_valid, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Neither diagonal nor inverse diagonal is in a valid state");
34: PetscCall(VecCopy(ctx->diag, ctx->inv_diag));
35: PetscCall(VecReciprocal(ctx->inv_diag));
36: ctx->inv_diag_valid = PETSC_TRUE;
37: }
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode MatAXPY_Diagonal(Mat Y, PetscScalar a, Mat X, MatStructure str)
42: {
43: Mat_Diagonal *yctx = (Mat_Diagonal *)Y->data;
44: Mat_Diagonal *xctx = (Mat_Diagonal *)X->data;
46: PetscFunctionBegin;
47: PetscCall(MatDiagonalSetUpDiagonal(Y));
48: PetscCall(MatDiagonalSetUpDiagonal(X));
49: PetscCall(VecAXPY(yctx->diag, a, xctx->diag));
50: yctx->inv_diag_valid = PETSC_FALSE;
51: PetscFunctionReturn(PETSC_SUCCESS);
52: }
54: static PetscErrorCode MatGetRow_Diagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
55: {
56: Mat_Diagonal *mat = (Mat_Diagonal *)A->data;
58: PetscFunctionBegin;
59: if (ncols) *ncols = 1;
60: if (cols) {
61: if (!mat->col) PetscCall(PetscMalloc1(1, &mat->col));
62: *mat->col = row;
63: *cols = mat->col;
64: }
65: if (vals) {
66: const PetscScalar *v;
68: if (!mat->val) PetscCall(PetscMalloc1(1, &mat->val));
69: PetscCall(VecGetArrayRead(mat->diag, &v));
70: *mat->val = v[row];
71: *vals = mat->val;
72: PetscCall(VecRestoreArrayRead(mat->diag, &v));
73: }
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode MatMult_Diagonal(Mat A, Vec x, Vec y)
78: {
79: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
81: PetscFunctionBegin;
82: PetscCall(MatDiagonalSetUpDiagonal(A));
83: PetscCall(VecPointwiseMult(y, ctx->diag, x));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode MatMultAdd_Diagonal(Mat mat, Vec v1, Vec v2, Vec v3)
88: {
89: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
91: PetscFunctionBegin;
92: PetscCall(MatDiagonalSetUpDiagonal(mat));
93: if (v2 != v3) {
94: PetscCall(VecPointwiseMult(v3, ctx->diag, v1));
95: PetscCall(VecAXPY(v3, 1.0, v2));
96: } else {
97: Vec w;
98: PetscCall(VecDuplicate(v3, &w));
99: PetscCall(VecPointwiseMult(w, ctx->diag, v1));
100: PetscCall(VecAXPY(v3, 1.0, w));
101: PetscCall(VecDestroy(&w));
102: }
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode MatNorm_Diagonal(Mat A, NormType type, PetscReal *nrm)
107: {
108: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
110: PetscFunctionBegin;
111: PetscCall(MatDiagonalSetUpDiagonal(A));
112: type = (type == NORM_FROBENIUS) ? NORM_2 : NORM_INFINITY;
113: PetscCall(VecNorm(ctx->diag, type, nrm));
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode MatDuplicate_Diagonal(Mat A, MatDuplicateOption op, Mat *B)
118: {
119: Mat_Diagonal *actx = (Mat_Diagonal *)A->data;
121: PetscFunctionBegin;
122: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
123: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
124: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
125: PetscCall(MatSetType(*B, MATDIAGONAL));
126: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
127: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
128: if (op == MAT_COPY_VALUES) {
129: Mat_Diagonal *bctx = (Mat_Diagonal *)(*B)->data;
131: PetscCall(MatSetUp(A));
132: PetscCall(MatSetUp(*B));
133: PetscCall(VecCopy(actx->diag, bctx->diag));
134: PetscCall(VecCopy(actx->inv_diag, bctx->inv_diag));
135: bctx->diag_valid = actx->diag_valid;
136: bctx->inv_diag_valid = actx->inv_diag_valid;
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@
142: MatDiagonalGetDiagonal - Get the diagonal of a `MATDIAGONAL`
144: Input Parameter:
145: . A - the `MATDIAGONAL`
147: Output Parameter:
148: . diag - the `Vec` that defines the diagonal
150: Level: developer
152: Note:
153: The user must call
154: `MatDiagonalRestoreDiagonal()` before using the matrix again.
156: For a copy of the diagonal values, rather than a reference, use `MatGetDiagonal()`
158: Any changes to the obtained vector immediately change the action of the `Mat`. The matrix can be changed more efficiently by accessing this vector and changing its values, instead of filling a work vector and using `MatDiagonalSet()`
160: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`, `MatGetDiagonal()`
161: @*/
162: PetscErrorCode MatDiagonalGetDiagonal(Mat A, Vec *diag)
163: {
164: PetscFunctionBegin;
166: PetscAssertPointer(diag, 2);
167: *diag = NULL;
168: PetscUseMethod(A, "MatDiagonalGetDiagonal_C", (Mat, Vec *), (A, diag));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode MatDiagonalGetDiagonal_Diagonal(Mat A, Vec *diag)
173: {
174: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
176: PetscFunctionBegin;
177: PetscCall(MatSetUp(A));
178: PetscCall(MatDiagonalSetUpDiagonal(A));
179: *diag = ctx->diag;
180: PetscCall(PetscObjectStateGet((PetscObject)*diag, &ctx->diag_state));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@
185: MatDiagonalRestoreDiagonal - Restore the diagonal of a `MATDIAGONAL`
187: Input Parameters:
188: + A - the `MATDIAGONAL`
189: - diag - the `Vec` obtained from `MatDiagonalGetDiagonal()`
191: Level: developer
193: Note:
194: Use `MatDiagonalSet()` to change the values by copy, rather than reference.
196: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetDiagonal()`
197: @*/
198: PetscErrorCode MatDiagonalRestoreDiagonal(Mat A, Vec *diag)
199: {
200: PetscFunctionBegin;
202: PetscAssertPointer(diag, 2);
203: PetscUseMethod(A, "MatDiagonalRestoreDiagonal_C", (Mat, Vec *), (A, diag));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: static PetscErrorCode MatDiagonalRestoreDiagonal_Diagonal(Mat A, Vec *diag)
208: {
209: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
210: PetscObjectState diag_state;
212: PetscFunctionBegin;
213: PetscCheck(ctx->diag == *diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
214: ctx->diag_valid = PETSC_TRUE;
215: PetscCall(PetscObjectStateGet((PetscObject)*diag, &diag_state));
216: if (ctx->diag_state != diag_state) {
217: PetscCall(PetscObjectStateIncrease((PetscObject)A));
218: ctx->inv_diag_valid = PETSC_FALSE;
219: }
220: *diag = NULL;
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*@
225: MatDiagonalGetInverseDiagonal - Get the inverse diagonal of a `MATDIAGONAL`
227: Input Parameter:
228: . A - the `MATDIAGONAL`
230: Output Parameter:
231: . inv_diag - the `Vec` that defines the inverse diagonal
233: Level: developer
235: Note:
236: The user must call
237: `MatDiagonalRestoreInverseDiagonal()` before using the matrix again.
239: If a matrix is created only to call `MatSolve()` (which happens for `MATLMVMDIAGBROYDEN`),
240: using `MatDiagonalGetInverseDiagonal()` and `MatDiagonalRestoreInverseDiagonal()` avoids copies
241: and avoids any call to `VecReciprocal()`.
243: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MATLMVMBROYDEN`, `MatSolve()`
244: @*/
245: PetscErrorCode MatDiagonalGetInverseDiagonal(Mat A, Vec *inv_diag)
246: {
247: PetscFunctionBegin;
249: PetscAssertPointer(inv_diag, 2);
250: *inv_diag = NULL;
251: PetscUseMethod(A, "MatDiagonalGetInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode MatDiagonalGetInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
256: {
257: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
259: PetscFunctionBegin;
260: PetscCall(MatSetUp(A));
261: PetscCall(MatDiagonalSetUpInverseDiagonal(A));
262: *inv_diag = ctx->inv_diag;
263: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &ctx->inv_diag_state));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: /*@
268: MatDiagonalRestoreInverseDiagonal - Restore the inverse diagonal of a `MATDIAGONAL`
270: Input Parameters:
271: + A - the `MATDIAGONAL`
272: - inv_diag - the `Vec` obtained from `MatDiagonalGetInverseDiagonal()`
274: Level: developer
276: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetInverseDiagonal()`
277: @*/
278: PetscErrorCode MatDiagonalRestoreInverseDiagonal(Mat A, Vec *inv_diag)
279: {
280: PetscFunctionBegin;
282: PetscAssertPointer(inv_diag, 2);
283: PetscUseMethod(A, "MatDiagonalRestoreInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode MatDiagonalRestoreInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
288: {
289: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
290: PetscObjectState inv_diag_state;
292: PetscFunctionBegin;
293: PetscCheck(ctx->inv_diag == *inv_diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
294: ctx->inv_diag_valid = PETSC_TRUE;
295: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &inv_diag_state));
296: if (ctx->inv_diag_state != inv_diag_state) {
297: PetscCall(PetscObjectStateIncrease((PetscObject)A));
298: ctx->diag_valid = PETSC_FALSE;
299: }
300: *inv_diag = NULL;
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode MatPermute_Diagonal(Mat A, IS rowp, IS colp, Mat *B)
305: {
306: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
307: Vec v;
309: PetscFunctionBegin;
310: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
311: PetscCall(VecDuplicate(ctx->diag, &v));
312: PetscCall(VecCopy(ctx->diag, v));
313: PetscCall(VecPermute(v, rowp, PETSC_FALSE));
314: PetscCall(MatCreateDiagonal(v, B));
315: PetscCall(VecDestroy(&v));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: static PetscErrorCode MatDestroy_Diagonal(Mat mat)
320: {
321: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
323: PetscFunctionBegin;
324: PetscCall(VecDestroy(&ctx->diag));
325: PetscCall(VecDestroy(&ctx->inv_diag));
326: PetscCall(PetscFree(ctx->col));
327: PetscCall(PetscFree(ctx->val));
328: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetDiagonal_C", NULL));
329: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreDiagonal_C", NULL));
330: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetInverseDiagonal_C", NULL));
331: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreInverseDiagonal_C", NULL));
332: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_seqdense_C", NULL));
333: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_mpidense_C", NULL));
334: PetscCall(PetscFree(mat->data));
335: mat->structural_symmetry_eternal = PETSC_FALSE;
336: mat->symmetry_eternal = PETSC_FALSE;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode MatView_Diagonal(Mat J, PetscViewer viewer)
341: {
342: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
343: PetscBool iascii;
345: PetscFunctionBegin;
346: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
347: if (iascii) {
348: PetscViewerFormat format;
350: PetscCall(PetscViewerGetFormat(viewer, &format));
351: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
352: PetscCall(MatDiagonalSetUpDiagonal(J));
353: PetscCall(VecView(ctx->diag, viewer));
354: }
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: static PetscErrorCode MatGetDiagonal_Diagonal(Mat J, Vec x)
359: {
360: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
362: PetscFunctionBegin;
363: PetscCall(MatDiagonalSetUpDiagonal(J));
364: PetscCall(VecCopy(ctx->diag, x));
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode MatDiagonalSet_Diagonal(Mat J, Vec D, InsertMode is)
369: {
370: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
372: PetscFunctionBegin;
373: switch (is) {
374: case ADD_VALUES:
375: case ADD_ALL_VALUES:
376: case ADD_BC_VALUES:
377: PetscCall(MatDiagonalSetUpDiagonal(J));
378: PetscCall(VecAXPY(ctx->diag, 1.0, D));
379: ctx->inv_diag_valid = PETSC_FALSE;
380: break;
381: case INSERT_VALUES:
382: case INSERT_BC_VALUES:
383: case INSERT_ALL_VALUES:
384: PetscCall(MatSetUp(J));
385: PetscCall(VecCopy(D, ctx->diag));
386: ctx->diag_valid = PETSC_TRUE;
387: ctx->inv_diag_valid = PETSC_FALSE;
388: break;
389: case MAX_VALUES:
390: PetscCall(MatDiagonalSetUpDiagonal(J));
391: PetscCall(VecPointwiseMax(ctx->diag, D, ctx->diag));
392: ctx->inv_diag_valid = PETSC_FALSE;
393: break;
394: case MIN_VALUES:
395: PetscCall(MatDiagonalSetUpDiagonal(J));
396: PetscCall(VecPointwiseMin(ctx->diag, D, ctx->diag));
397: ctx->inv_diag_valid = PETSC_FALSE;
398: break;
399: case NOT_SET_VALUES:
400: break;
401: }
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode MatShift_Diagonal(Mat Y, PetscScalar a)
406: {
407: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
409: PetscFunctionBegin;
410: PetscCall(MatDiagonalSetUpDiagonal(Y));
411: PetscCall(VecShift(ctx->diag, a));
412: ctx->inv_diag_valid = PETSC_FALSE;
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode MatScale_Diagonal(Mat Y, PetscScalar a)
417: {
418: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
420: PetscFunctionBegin;
421: PetscCall(MatSetUp(Y));
422: PetscCall(MatDiagonalSetUpDiagonal(Y));
423: PetscCall(VecScale(ctx->diag, a));
424: ctx->inv_diag_valid = PETSC_FALSE;
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode MatDiagonalScale_Diagonal(Mat Y, Vec l, Vec r)
429: {
430: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
432: PetscFunctionBegin;
433: PetscCall(MatDiagonalSetUpDiagonal(Y));
434: if (l) {
435: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, l));
436: ctx->inv_diag_valid = PETSC_FALSE;
437: }
438: if (r) {
439: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, r));
440: ctx->inv_diag_valid = PETSC_FALSE;
441: }
442: PetscFunctionReturn(PETSC_SUCCESS);
443: }
445: static PetscErrorCode MatSetUp_Diagonal(Mat A)
446: {
447: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
449: PetscFunctionBegin;
450: if (!ctx->diag) {
451: PetscCall(PetscLayoutSetUp(A->cmap));
452: PetscCall(PetscLayoutSetUp(A->rmap));
453: PetscCall(MatCreateVecs(A, &ctx->diag, NULL));
454: PetscCall(VecDuplicate(ctx->diag, &ctx->inv_diag));
455: PetscCall(VecZeroEntries(ctx->diag));
456: ctx->diag_valid = PETSC_TRUE;
457: }
458: A->assembled = PETSC_TRUE;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: static PetscErrorCode MatZeroEntries_Diagonal(Mat Y)
463: {
464: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
466: PetscFunctionBegin;
467: PetscCall(MatSetUp(Y));
468: PetscCall(VecZeroEntries(ctx->diag));
469: ctx->diag_valid = PETSC_TRUE;
470: ctx->inv_diag_valid = PETSC_FALSE;
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: static PetscErrorCode MatSolve_Diagonal(Mat matin, Vec b, Vec x)
475: {
476: Mat_Diagonal *ctx = (Mat_Diagonal *)matin->data;
478: PetscFunctionBegin;
479: PetscCall(MatDiagonalSetUpInverseDiagonal(matin));
480: PetscCall(VecPointwiseMult(x, b, ctx->inv_diag));
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: static PetscErrorCode MatGetInfo_Diagonal(Mat A, MatInfoType flag, MatInfo *info)
485: {
486: PetscFunctionBegin;
487: info->block_size = 1.0;
488: info->nz_allocated = A->cmap->N;
489: info->nz_used = A->cmap->N;
490: info->nz_unneeded = 0.0;
491: info->assemblies = A->num_ass;
492: info->mallocs = 0.0;
493: info->memory = 0;
494: info->fill_ratio_given = 0;
495: info->fill_ratio_needed = 0;
496: info->factor_mallocs = 0;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: MatCreateDiagonal - Creates a matrix defined by a given vector along its diagonal.
503: Collective
505: Input Parameter:
506: . diag - vector for the diagonal
508: Output Parameter:
509: . J - the diagonal matrix
511: Level: advanced
513: Notes:
514: Only supports square matrices with the same number of local rows and columns.
516: The input vector `diag` will be referenced internally: any changes to `diag`
517: will affect the matrix `J`.
519: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatSolve()`
520: `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`
521: @*/
522: PetscErrorCode MatCreateDiagonal(Vec diag, Mat *J)
523: {
524: PetscFunctionBegin;
526: PetscCall(MatCreate(PetscObjectComm((PetscObject)diag), J));
527: PetscInt m, M;
528: PetscCall(VecGetLocalSize(diag, &m));
529: PetscCall(VecGetSize(diag, &M));
530: PetscCall(MatSetSizes(*J, m, m, M, M));
531: PetscCall(MatSetType(*J, MATDIAGONAL));
533: PetscLayout map;
534: PetscCall(VecGetLayout(diag, &map));
535: PetscCall(MatSetLayouts(*J, map, map));
536: Mat_Diagonal *ctx = (Mat_Diagonal *)(*J)->data;
537: PetscCall(PetscObjectReference((PetscObject)diag));
538: PetscCall(VecDestroy(&ctx->diag));
539: PetscCall(VecDestroy(&ctx->inv_diag));
540: ctx->diag = diag;
541: ctx->diag_valid = PETSC_TRUE;
542: ctx->inv_diag_valid = PETSC_FALSE;
543: VecType type;
544: PetscCall(VecDuplicate(ctx->diag, &ctx->inv_diag));
545: PetscCall(VecGetType(diag, &type));
546: PetscCall(PetscFree((*J)->defaultvectype));
547: PetscCall(PetscStrallocpy(type, &(*J)->defaultvectype));
548: PetscCall(MatSetUp(*J));
549: ctx->col = NULL;
550: ctx->val = NULL;
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode MatProductNumeric_Diagonal_Dense(Mat C)
555: {
556: Mat A, B;
557: Mat_Diagonal *a;
558: PetscScalar *c;
559: const PetscScalar *b, *alpha;
560: PetscInt ldb, ldc;
562: PetscFunctionBegin;
563: MatCheckProduct(C, 1);
564: A = C->product->A;
565: B = C->product->B;
566: a = (Mat_Diagonal *)A->data;
567: PetscCall(VecGetArrayRead(a->diag, &alpha));
568: PetscCall(MatDenseGetLDA(B, &ldb));
569: PetscCall(MatDenseGetLDA(C, &ldc));
570: PetscCall(MatDenseGetArrayRead(B, &b));
571: PetscCall(MatDenseGetArrayWrite(C, &c));
572: for (PetscInt j = 0; j < B->cmap->N; j++)
573: for (PetscInt i = 0; i < B->rmap->n; i++) c[i + j * ldc] = alpha[i] * b[i + j * ldb];
574: PetscCall(MatDenseRestoreArrayWrite(C, &c));
575: PetscCall(MatDenseRestoreArrayRead(B, &b));
576: PetscCall(VecRestoreArrayRead(a->diag, &alpha));
577: PetscCall(PetscLogFlops(B->cmap->N * B->rmap->n));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: static PetscErrorCode MatProductSymbolic_Diagonal_Dense(Mat C)
582: {
583: Mat A, B;
584: PetscInt n, N, m, M;
586: PetscFunctionBegin;
587: MatCheckProduct(C, 1);
588: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
589: A = C->product->A;
590: B = C->product->B;
591: PetscCall(MatDiagonalSetUpDiagonal(A));
592: PetscCall(MatGetLocalSize(C, &m, &n));
593: PetscCall(MatGetSize(C, &M, &N));
594: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
595: PetscCall(MatGetLocalSize(B, NULL, &n));
596: PetscCall(MatGetSize(B, NULL, &N));
597: PetscCall(MatGetLocalSize(A, &m, NULL));
598: PetscCall(MatGetSize(A, &M, NULL));
599: PetscCall(MatSetSizes(C, m, n, M, N));
600: }
601: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
602: PetscCall(MatSetUp(C));
603: C->ops->productnumeric = MatProductNumeric_Diagonal_Dense;
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense_AB(Mat C)
608: {
609: PetscFunctionBegin;
610: C->ops->productsymbolic = MatProductSymbolic_Diagonal_Dense;
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense(Mat C)
615: {
616: Mat_Product *product = C->product;
618: PetscFunctionBegin;
619: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_Diagonal_Dense_AB(C));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: /*MC
624: MATDIAGONAL - MATDIAGONAL = "diagonal" - A diagonal matrix type with the diagonal implemented as a `Vec`. Useful for
625: cases where `VecPointwiseMult()` or `VecPointwiseDivide()` should be thought of as the actions of a linear operator.
627: Level: advanced
629: .seealso: [](ch_matrices), `Mat`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`
630: M*/
631: PETSC_INTERN PetscErrorCode MatCreate_Diagonal(Mat A)
632: {
633: Mat_Diagonal *ctx;
635: PetscFunctionBegin;
636: PetscCall(PetscNew(&ctx));
637: A->data = (void *)ctx;
639: A->structurally_symmetric = PETSC_BOOL3_TRUE;
640: A->structural_symmetry_eternal = PETSC_TRUE;
641: A->symmetry_eternal = PETSC_TRUE;
642: A->symmetric = PETSC_BOOL3_TRUE;
643: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
645: A->ops->getrow = MatGetRow_Diagonal;
646: A->ops->mult = MatMult_Diagonal;
647: A->ops->multadd = MatMultAdd_Diagonal;
648: A->ops->multtranspose = MatMult_Diagonal;
649: A->ops->multtransposeadd = MatMultAdd_Diagonal;
650: A->ops->norm = MatNorm_Diagonal;
651: A->ops->duplicate = MatDuplicate_Diagonal;
652: A->ops->solve = MatSolve_Diagonal;
653: A->ops->solvetranspose = MatSolve_Diagonal;
654: A->ops->shift = MatShift_Diagonal;
655: A->ops->scale = MatScale_Diagonal;
656: A->ops->diagonalscale = MatDiagonalScale_Diagonal;
657: A->ops->getdiagonal = MatGetDiagonal_Diagonal;
658: A->ops->diagonalset = MatDiagonalSet_Diagonal;
659: A->ops->view = MatView_Diagonal;
660: A->ops->zeroentries = MatZeroEntries_Diagonal;
661: A->ops->destroy = MatDestroy_Diagonal;
662: A->ops->getinfo = MatGetInfo_Diagonal;
663: A->ops->axpy = MatAXPY_Diagonal;
664: A->ops->setup = MatSetUp_Diagonal;
665: A->ops->permute = MatPermute_Diagonal;
667: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetDiagonal_C", MatDiagonalGetDiagonal_Diagonal));
668: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreDiagonal_C", MatDiagonalRestoreDiagonal_Diagonal));
669: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetInverseDiagonal_C", MatDiagonalGetInverseDiagonal_Diagonal));
670: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreInverseDiagonal_C", MatDiagonalRestoreInverseDiagonal_Diagonal));
671: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_seqdense_C", MatProductSetFromOptions_Diagonal_Dense));
672: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_mpidense_C", MatProductSetFromOptions_Diagonal_Dense));
673: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATDIAGONAL));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }