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 MatRestoreRow_Diagonal(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
78: {
79: PetscFunctionBegin;
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
83: static PetscErrorCode MatMult_Diagonal(Mat A, Vec x, Vec y)
84: {
85: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
87: PetscFunctionBegin;
88: PetscCall(MatDiagonalSetUpDiagonal(A));
89: PetscCall(VecPointwiseMult(y, ctx->diag, x));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode MatMultAdd_Diagonal(Mat mat, Vec v1, Vec v2, Vec v3)
94: {
95: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
97: PetscFunctionBegin;
98: PetscCall(MatDiagonalSetUpDiagonal(mat));
99: if (v2 != v3) {
100: PetscCall(VecPointwiseMult(v3, ctx->diag, v1));
101: PetscCall(VecAXPY(v3, 1.0, v2));
102: } else {
103: Vec w;
104: PetscCall(VecDuplicate(v3, &w));
105: PetscCall(VecPointwiseMult(w, ctx->diag, v1));
106: PetscCall(VecAXPY(v3, 1.0, w));
107: PetscCall(VecDestroy(&w));
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode MatNorm_Diagonal(Mat A, NormType type, PetscReal *nrm)
113: {
114: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
116: PetscFunctionBegin;
117: PetscCall(MatDiagonalSetUpDiagonal(A));
118: type = (type == NORM_FROBENIUS) ? NORM_2 : NORM_INFINITY;
119: PetscCall(VecNorm(ctx->diag, type, nrm));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: static PetscErrorCode MatDuplicate_Diagonal(Mat A, MatDuplicateOption op, Mat *B)
124: {
125: Mat_Diagonal *actx = (Mat_Diagonal *)A->data;
127: PetscFunctionBegin;
128: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
129: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
130: PetscCall(MatSetBlockSizesFromMats(*B, A, A));
131: PetscCall(MatSetType(*B, MATDIAGONAL));
132: PetscCall(PetscLayoutReference(A->rmap, &(*B)->rmap));
133: PetscCall(PetscLayoutReference(A->cmap, &(*B)->cmap));
134: if (op == MAT_COPY_VALUES) {
135: Mat_Diagonal *bctx = (Mat_Diagonal *)(*B)->data;
137: PetscCall(MatSetUp(A));
138: PetscCall(MatSetUp(*B));
139: PetscCall(VecCopy(actx->diag, bctx->diag));
140: PetscCall(VecCopy(actx->inv_diag, bctx->inv_diag));
141: bctx->diag_valid = actx->diag_valid;
142: bctx->inv_diag_valid = actx->inv_diag_valid;
143: }
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: /*@
148: MatDiagonalGetDiagonal - Get the diagonal of a `MATDIAGONAL`
150: Input Parameter:
151: . A - the `MATDIAGONAL`
153: Output Parameter:
154: . diag - the `Vec` that defines the diagonal
156: Level: developer
158: Note:
159: The user must call
160: `MatDiagonalRestoreDiagonal()` before using the matrix again.
162: For a copy of the diagonal values, rather than a reference, use `MatGetDiagonal()`
164: 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()`
166: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`, `MatGetDiagonal()`
167: @*/
168: PetscErrorCode MatDiagonalGetDiagonal(Mat A, Vec *diag)
169: {
170: PetscFunctionBegin;
172: PetscAssertPointer(diag, 2);
173: *diag = NULL;
174: PetscUseMethod((PetscObject)A, "MatDiagonalGetDiagonal_C", (Mat, Vec *), (A, diag));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode MatDiagonalGetDiagonal_Diagonal(Mat A, Vec *diag)
179: {
180: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
182: PetscFunctionBegin;
183: PetscCall(MatSetUp(A));
184: PetscCall(MatDiagonalSetUpDiagonal(A));
185: *diag = ctx->diag;
186: PetscCall(PetscObjectStateGet((PetscObject)*diag, &ctx->diag_state));
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: /*@
191: MatDiagonalRestoreDiagonal - Restore the diagonal of a `MATDIAGONAL`
193: Input Parameters:
194: + A - the `MATDIAGONAL`
195: - diag - the `Vec` obtained from `MatDiagonalGetDiagonal()`
197: Level: developer
199: Note:
200: Use `MatDiagonalSet()` to change the values by copy, rather than reference.
202: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetDiagonal()`
203: @*/
204: PetscErrorCode MatDiagonalRestoreDiagonal(Mat A, Vec *diag)
205: {
206: PetscFunctionBegin;
208: PetscAssertPointer(diag, 2);
209: PetscUseMethod((PetscObject)A, "MatDiagonalRestoreDiagonal_C", (Mat, Vec *), (A, diag));
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: static PetscErrorCode MatDiagonalRestoreDiagonal_Diagonal(Mat A, Vec *diag)
214: {
215: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
216: PetscObjectState diag_state;
218: PetscFunctionBegin;
219: PetscCheck(ctx->diag == *diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
220: ctx->diag_valid = PETSC_TRUE;
221: PetscCall(PetscObjectStateGet((PetscObject)*diag, &diag_state));
222: if (ctx->diag_state != diag_state) {
223: PetscCall(PetscObjectStateIncrease((PetscObject)A));
224: ctx->inv_diag_valid = PETSC_FALSE;
225: }
226: *diag = NULL;
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /*@
231: MatDiagonalGetInverseDiagonal - Get the inverse diagonal of a `MATDIAGONAL`
233: Input Parameter:
234: . A - the `MATDIAGONAL`
236: Output Parameter:
237: . inv_diag - the `Vec` that defines the inverse diagonal
239: Level: developer
241: Note:
242: The user must call
243: `MatDiagonalRestoreInverseDiagonal()` before using the matrix again.
245: If a matrix is created only to call `MatSolve()` (which happens for `MATLMVMDIAGBROYDEN`),
246: using `MatDiagonalGetInverseDiagonal()` and `MatDiagonalRestoreInverseDiagonal()` avoids copies
247: and avoids any call to `VecReciprocal()`.
249: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MATLMVMBROYDEN`, `MatSolve()`
250: @*/
251: PetscErrorCode MatDiagonalGetInverseDiagonal(Mat A, Vec *inv_diag)
252: {
253: PetscFunctionBegin;
255: PetscAssertPointer(inv_diag, 2);
256: *inv_diag = NULL;
257: PetscUseMethod((PetscObject)A, "MatDiagonalGetInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode MatDiagonalGetInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
262: {
263: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
265: PetscFunctionBegin;
266: PetscCall(MatSetUp(A));
267: PetscCall(MatDiagonalSetUpInverseDiagonal(A));
268: *inv_diag = ctx->inv_diag;
269: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &ctx->inv_diag_state));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: MatDiagonalRestoreInverseDiagonal - Restore the inverse diagonal of a `MATDIAGONAL`
276: Input Parameters:
277: + A - the `MATDIAGONAL`
278: - inv_diag - the `Vec` obtained from `MatDiagonalGetInverseDiagonal()`
280: Level: developer
282: .seealso: [](ch_matrices), `MATDIAGONAL`, `MatCreateDiagonal()`, `MatDiagonalGetInverseDiagonal()`
283: @*/
284: PetscErrorCode MatDiagonalRestoreInverseDiagonal(Mat A, Vec *inv_diag)
285: {
286: PetscFunctionBegin;
288: PetscAssertPointer(inv_diag, 2);
289: PetscUseMethod((PetscObject)A, "MatDiagonalRestoreInverseDiagonal_C", (Mat, Vec *), (A, inv_diag));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode MatDiagonalRestoreInverseDiagonal_Diagonal(Mat A, Vec *inv_diag)
294: {
295: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
296: PetscObjectState inv_diag_state;
298: PetscFunctionBegin;
299: PetscCheck(ctx->inv_diag == *inv_diag, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Restored a different diagonal vector");
300: ctx->inv_diag_valid = PETSC_TRUE;
301: PetscCall(PetscObjectStateGet((PetscObject)*inv_diag, &inv_diag_state));
302: if (ctx->inv_diag_state != inv_diag_state) {
303: PetscCall(PetscObjectStateIncrease((PetscObject)A));
304: ctx->diag_valid = PETSC_FALSE;
305: }
306: *inv_diag = NULL;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode MatPermute_Diagonal(Mat A, IS rowp, IS colp, Mat *B)
311: {
312: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
313: Vec v;
315: PetscFunctionBegin;
316: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
317: PetscCall(VecDuplicate(ctx->diag, &v));
318: PetscCall(VecCopy(ctx->diag, v));
319: PetscCall(VecPermute(v, rowp, PETSC_FALSE));
320: PetscCall(MatCreateDiagonal(v, B));
321: PetscCall(VecDestroy(&v));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode MatDestroy_Diagonal(Mat mat)
326: {
327: Mat_Diagonal *ctx = (Mat_Diagonal *)mat->data;
329: PetscFunctionBegin;
330: PetscCall(VecDestroy(&ctx->diag));
331: PetscCall(VecDestroy(&ctx->inv_diag));
332: PetscCall(PetscFree(ctx->col));
333: PetscCall(PetscFree(ctx->val));
334: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetDiagonal_C", NULL));
335: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreDiagonal_C", NULL));
336: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalGetInverseDiagonal_C", NULL));
337: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDiagonalRestoreInverseDiagonal_C", NULL));
338: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_seqdense_C", NULL));
339: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_diagonal_mpidense_C", NULL));
340: PetscCall(PetscFree(mat->data));
341: mat->structural_symmetry_eternal = PETSC_FALSE;
342: mat->symmetry_eternal = PETSC_FALSE;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: static PetscErrorCode MatView_Diagonal(Mat J, PetscViewer viewer)
347: {
348: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
349: PetscBool iascii;
351: PetscFunctionBegin;
352: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
353: if (iascii) {
354: PetscViewerFormat format;
356: PetscCall(PetscViewerGetFormat(viewer, &format));
357: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
358: PetscCall(MatDiagonalSetUpDiagonal(J));
359: PetscCall(VecView(ctx->diag, viewer));
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode MatGetDiagonal_Diagonal(Mat J, Vec x)
365: {
366: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
368: PetscFunctionBegin;
369: PetscCall(MatDiagonalSetUpDiagonal(J));
370: PetscCall(VecCopy(ctx->diag, x));
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode MatDiagonalSet_Diagonal(Mat J, Vec D, InsertMode is)
375: {
376: Mat_Diagonal *ctx = (Mat_Diagonal *)J->data;
378: PetscFunctionBegin;
379: switch (is) {
380: case ADD_VALUES:
381: case ADD_ALL_VALUES:
382: case ADD_BC_VALUES:
383: PetscCall(MatDiagonalSetUpDiagonal(J));
384: PetscCall(VecAXPY(ctx->diag, 1.0, D));
385: ctx->inv_diag_valid = PETSC_FALSE;
386: break;
387: case INSERT_VALUES:
388: case INSERT_BC_VALUES:
389: case INSERT_ALL_VALUES:
390: PetscCall(MatSetUp(J));
391: PetscCall(VecCopy(D, ctx->diag));
392: ctx->diag_valid = PETSC_TRUE;
393: ctx->inv_diag_valid = PETSC_FALSE;
394: break;
395: case MAX_VALUES:
396: PetscCall(MatDiagonalSetUpDiagonal(J));
397: PetscCall(VecPointwiseMax(ctx->diag, D, ctx->diag));
398: ctx->inv_diag_valid = PETSC_FALSE;
399: break;
400: case MIN_VALUES:
401: PetscCall(MatDiagonalSetUpDiagonal(J));
402: PetscCall(VecPointwiseMin(ctx->diag, D, ctx->diag));
403: ctx->inv_diag_valid = PETSC_FALSE;
404: break;
405: case NOT_SET_VALUES:
406: break;
407: }
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: static PetscErrorCode MatShift_Diagonal(Mat Y, PetscScalar a)
412: {
413: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
415: PetscFunctionBegin;
416: PetscCall(MatDiagonalSetUpDiagonal(Y));
417: PetscCall(VecShift(ctx->diag, a));
418: ctx->inv_diag_valid = PETSC_FALSE;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode MatScale_Diagonal(Mat Y, PetscScalar a)
423: {
424: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
426: PetscFunctionBegin;
427: PetscCall(MatSetUp(Y));
428: PetscCall(MatDiagonalSetUpDiagonal(Y));
429: PetscCall(VecScale(ctx->diag, a));
430: ctx->inv_diag_valid = PETSC_FALSE;
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: static PetscErrorCode MatDiagonalScale_Diagonal(Mat Y, Vec l, Vec r)
435: {
436: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
438: PetscFunctionBegin;
439: PetscCall(MatDiagonalSetUpDiagonal(Y));
440: if (l) {
441: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, l));
442: ctx->inv_diag_valid = PETSC_FALSE;
443: }
444: if (r) {
445: PetscCall(VecPointwiseMult(ctx->diag, ctx->diag, r));
446: ctx->inv_diag_valid = PETSC_FALSE;
447: }
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode MatSetUp_Diagonal(Mat A)
452: {
453: Mat_Diagonal *ctx = (Mat_Diagonal *)A->data;
455: PetscFunctionBegin;
456: if (!ctx->diag) {
457: PetscCall(PetscLayoutSetUp(A->cmap));
458: PetscCall(PetscLayoutSetUp(A->rmap));
459: PetscCall(MatCreateVecs(A, &ctx->diag, NULL));
460: PetscCall(VecDuplicate(ctx->diag, &ctx->inv_diag));
461: PetscCall(VecZeroEntries(ctx->diag));
462: ctx->diag_valid = PETSC_TRUE;
463: }
464: A->assembled = PETSC_TRUE;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: static PetscErrorCode MatZeroEntries_Diagonal(Mat Y)
469: {
470: Mat_Diagonal *ctx = (Mat_Diagonal *)Y->data;
472: PetscFunctionBegin;
473: PetscCall(MatSetUp(Y));
474: PetscCall(VecZeroEntries(ctx->diag));
475: ctx->diag_valid = PETSC_TRUE;
476: ctx->inv_diag_valid = PETSC_FALSE;
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode MatSolve_Diagonal(Mat matin, Vec b, Vec x)
481: {
482: Mat_Diagonal *ctx = (Mat_Diagonal *)matin->data;
484: PetscFunctionBegin;
485: PetscCall(MatDiagonalSetUpInverseDiagonal(matin));
486: PetscCall(VecPointwiseMult(x, b, ctx->inv_diag));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode MatGetInfo_Diagonal(Mat A, MatInfoType flag, MatInfo *info)
491: {
492: PetscFunctionBegin;
493: info->block_size = 1.0;
494: info->nz_allocated = A->cmap->N;
495: info->nz_used = A->cmap->N;
496: info->nz_unneeded = 0.0;
497: info->assemblies = A->num_ass;
498: info->mallocs = 0.0;
499: info->memory = 0;
500: info->fill_ratio_given = 0;
501: info->fill_ratio_needed = 0;
502: info->factor_mallocs = 0;
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: /*@
507: MatCreateDiagonal - Creates a matrix defined by a given vector along its diagonal.
509: Collective
511: Input Parameter:
512: . diag - vector for the diagonal
514: Output Parameter:
515: . J - the diagonal matrix
517: Level: advanced
519: Notes:
520: Only supports square matrices with the same number of local rows and columns.
522: The input vector `diag` will be referenced internally: any changes to `diag`
523: will affect the matrix `J`.
525: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MATCONSTANTDIAGONAL`, `MatScale()`, `MatShift()`, `MatMult()`, `MatGetDiagonal()`, `MatSolve()`
526: `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`
527: @*/
528: PetscErrorCode MatCreateDiagonal(Vec diag, Mat *J)
529: {
530: PetscFunctionBegin;
532: PetscCall(MatCreate(PetscObjectComm((PetscObject)diag), J));
533: PetscInt m, M;
534: PetscCall(VecGetLocalSize(diag, &m));
535: PetscCall(VecGetSize(diag, &M));
536: PetscCall(MatSetSizes(*J, m, m, M, M));
537: PetscCall(MatSetType(*J, MATDIAGONAL));
539: PetscLayout map;
540: PetscCall(VecGetLayout(diag, &map));
541: PetscCall(MatSetLayouts(*J, map, map));
542: Mat_Diagonal *ctx = (Mat_Diagonal *)(*J)->data;
543: PetscCall(PetscObjectReference((PetscObject)diag));
544: PetscCall(VecDestroy(&ctx->diag));
545: PetscCall(VecDestroy(&ctx->inv_diag));
546: ctx->diag = diag;
547: ctx->diag_valid = PETSC_TRUE;
548: ctx->inv_diag_valid = PETSC_FALSE;
549: VecType type;
550: PetscCall(VecDuplicate(ctx->diag, &ctx->inv_diag));
551: PetscCall(VecGetType(diag, &type));
552: PetscCall(PetscFree((*J)->defaultvectype));
553: PetscCall(PetscStrallocpy(type, &(*J)->defaultvectype));
554: PetscCall(MatSetUp(*J));
555: ctx->col = NULL;
556: ctx->val = NULL;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode MatProductNumeric_Diagonal_Dense(Mat C)
561: {
562: Mat A, B;
563: Mat_Diagonal *a;
564: PetscScalar *c;
565: const PetscScalar *b, *alpha;
566: PetscInt ldb, ldc;
568: PetscFunctionBegin;
569: MatCheckProduct(C, 1);
570: A = C->product->A;
571: B = C->product->B;
572: a = (Mat_Diagonal *)A->data;
573: PetscCall(VecGetArrayRead(a->diag, &alpha));
574: PetscCall(MatDenseGetLDA(B, &ldb));
575: PetscCall(MatDenseGetLDA(C, &ldc));
576: PetscCall(MatDenseGetArrayRead(B, &b));
577: PetscCall(MatDenseGetArrayWrite(C, &c));
578: for (PetscInt j = 0; j < B->cmap->N; j++)
579: for (PetscInt i = 0; i < B->rmap->n; i++) c[i + j * ldc] = alpha[i] * b[i + j * ldb];
580: PetscCall(MatDenseRestoreArrayWrite(C, &c));
581: PetscCall(MatDenseRestoreArrayRead(B, &b));
582: PetscCall(VecRestoreArrayRead(a->diag, &alpha));
583: PetscCall(PetscLogFlops(B->cmap->N * B->rmap->n));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: static PetscErrorCode MatProductSymbolic_Diagonal_Dense(Mat C)
588: {
589: Mat A, B;
590: PetscInt n, N, m, M;
592: PetscFunctionBegin;
593: MatCheckProduct(C, 1);
594: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
595: A = C->product->A;
596: B = C->product->B;
597: PetscCall(MatDiagonalSetUpDiagonal(A));
598: PetscCall(MatGetLocalSize(C, &m, &n));
599: PetscCall(MatGetSize(C, &M, &N));
600: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
601: PetscCall(MatGetLocalSize(B, NULL, &n));
602: PetscCall(MatGetSize(B, NULL, &N));
603: PetscCall(MatGetLocalSize(A, &m, NULL));
604: PetscCall(MatGetSize(A, &M, NULL));
605: PetscCall(MatSetSizes(C, m, n, M, N));
606: }
607: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
608: PetscCall(MatSetUp(C));
609: C->ops->productnumeric = MatProductNumeric_Diagonal_Dense;
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense_AB(Mat C)
614: {
615: PetscFunctionBegin;
616: C->ops->productsymbolic = MatProductSymbolic_Diagonal_Dense;
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: static PetscErrorCode MatProductSetFromOptions_Diagonal_Dense(Mat C)
621: {
622: Mat_Product *product = C->product;
624: PetscFunctionBegin;
625: if (product->type == MATPRODUCT_AB || product->type == MATPRODUCT_AtB) PetscCall(MatProductSetFromOptions_Diagonal_Dense_AB(C));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /*MC
630: MATDIAGONAL - MATDIAGONAL = "diagonal" - A diagonal matrix type with the diagonal implemented as a `Vec`. Useful for
631: cases where `VecPointwiseMult()` or `VecPointwiseDivide()` should be thought of as the actions of a linear operator.
633: Level: advanced
635: .seealso: [](ch_matrices), `Mat`, `MatCreateDiagonal()`, `MatDiagonalRestoreInverseDiagonal()`, `MatDiagonalGetDiagonal()`, `MatDiagonalRestoreDiagonal()`, `MatDiagonalGetInverseDiagonal()`
636: M*/
637: PETSC_INTERN PetscErrorCode MatCreate_Diagonal(Mat A)
638: {
639: Mat_Diagonal *ctx;
641: PetscFunctionBegin;
642: PetscCall(PetscNew(&ctx));
643: A->data = (void *)ctx;
645: A->structurally_symmetric = PETSC_BOOL3_TRUE;
646: A->structural_symmetry_eternal = PETSC_TRUE;
647: A->symmetry_eternal = PETSC_TRUE;
648: A->symmetric = PETSC_BOOL3_TRUE;
649: if (!PetscDefined(USE_COMPLEX)) A->hermitian = PETSC_BOOL3_TRUE;
651: A->ops->getrow = MatGetRow_Diagonal;
652: A->ops->restorerow = MatRestoreRow_Diagonal;
653: A->ops->mult = MatMult_Diagonal;
654: A->ops->multadd = MatMultAdd_Diagonal;
655: A->ops->multtranspose = MatMult_Diagonal;
656: A->ops->multtransposeadd = MatMultAdd_Diagonal;
657: A->ops->norm = MatNorm_Diagonal;
658: A->ops->duplicate = MatDuplicate_Diagonal;
659: A->ops->solve = MatSolve_Diagonal;
660: A->ops->solvetranspose = MatSolve_Diagonal;
661: A->ops->shift = MatShift_Diagonal;
662: A->ops->scale = MatScale_Diagonal;
663: A->ops->diagonalscale = MatDiagonalScale_Diagonal;
664: A->ops->getdiagonal = MatGetDiagonal_Diagonal;
665: A->ops->diagonalset = MatDiagonalSet_Diagonal;
666: A->ops->view = MatView_Diagonal;
667: A->ops->zeroentries = MatZeroEntries_Diagonal;
668: A->ops->destroy = MatDestroy_Diagonal;
669: A->ops->getinfo = MatGetInfo_Diagonal;
670: A->ops->axpy = MatAXPY_Diagonal;
671: A->ops->setup = MatSetUp_Diagonal;
672: A->ops->permute = MatPermute_Diagonal;
674: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetDiagonal_C", MatDiagonalGetDiagonal_Diagonal));
675: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreDiagonal_C", MatDiagonalRestoreDiagonal_Diagonal));
676: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalGetInverseDiagonal_C", MatDiagonalGetInverseDiagonal_Diagonal));
677: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatDiagonalRestoreInverseDiagonal_C", MatDiagonalRestoreInverseDiagonal_Diagonal));
678: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_seqdense_C", MatProductSetFromOptions_Diagonal_Dense));
679: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_diagonal_mpidense_C", MatProductSetFromOptions_Diagonal_Dense));
680: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATDIAGONAL));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }