Actual source code: mcomposite.c
1: #include <petsc/private/matimpl.h>
3: const char *const MatCompositeMergeTypes[] = {"left", "right", "MatCompositeMergeType", "MAT_COMPOSITE_", NULL};
5: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
6: struct _Mat_CompositeLink {
7: Mat mat;
8: Vec work;
9: Mat_CompositeLink next, prev;
10: };
12: typedef struct {
13: MatCompositeType type;
14: Mat_CompositeLink head, tail;
15: Vec work;
16: PetscScalar scale; /* scale factor supplied with MatScale() */
17: Vec left, right; /* left and right diagonal scaling provided with MatDiagonalScale() */
18: Vec leftwork, rightwork, leftwork2, rightwork2; /* Two pairs of working vectors */
19: PetscInt nmat;
20: PetscBool merge;
21: MatCompositeMergeType mergetype;
22: MatStructure structure;
24: PetscScalar *scalings;
25: PetscBool merge_mvctx; /* Whether need to merge mvctx of component matrices */
26: Vec *lvecs; /* [nmat] Basically, they are Mvctx->lvec of each component matrix */
27: PetscScalar *larray; /* [len] Data arrays of lvecs[] are stored consecutively in larray */
28: PetscInt len; /* Length of larray[] */
29: Vec gvec; /* Union of lvecs[] without duplicated entries */
30: PetscInt *location; /* A map that maps entries in garray[] to larray[] */
31: VecScatter Mvctx;
32: } Mat_Composite;
34: static PetscErrorCode MatDestroy_Composite(Mat mat)
35: {
36: Mat_Composite *shell = (Mat_Composite *)mat->data;
37: Mat_CompositeLink next = shell->head, oldnext;
38: PetscInt i;
40: PetscFunctionBegin;
41: while (next) {
42: PetscCall(MatDestroy(&next->mat));
43: if (next->work && (!next->next || next->work != next->next->work)) PetscCall(VecDestroy(&next->work));
44: oldnext = next;
45: next = next->next;
46: PetscCall(PetscFree(oldnext));
47: }
48: PetscCall(VecDestroy(&shell->work));
49: PetscCall(VecDestroy(&shell->left));
50: PetscCall(VecDestroy(&shell->right));
51: PetscCall(VecDestroy(&shell->leftwork));
52: PetscCall(VecDestroy(&shell->rightwork));
53: PetscCall(VecDestroy(&shell->leftwork2));
54: PetscCall(VecDestroy(&shell->rightwork2));
56: if (shell->Mvctx) {
57: for (i = 0; i < shell->nmat; i++) PetscCall(VecDestroy(&shell->lvecs[i]));
58: PetscCall(PetscFree3(shell->location, shell->larray, shell->lvecs));
59: PetscCall(PetscFree(shell->larray));
60: PetscCall(VecDestroy(&shell->gvec));
61: PetscCall(VecScatterDestroy(&shell->Mvctx));
62: }
64: PetscCall(PetscFree(shell->scalings));
65: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL));
66: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL));
67: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL));
68: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL));
69: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL));
70: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL));
71: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL));
72: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL));
73: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL));
74: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL));
75: PetscCall(PetscFree(mat->data));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode MatMult_Composite_Multiplicative(Mat A, Vec x, Vec y)
80: {
81: Mat_Composite *shell = (Mat_Composite *)A->data;
82: Mat_CompositeLink next = shell->head;
83: Vec in, out;
84: PetscScalar scale;
85: PetscInt i;
87: PetscFunctionBegin;
88: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
89: in = x;
90: if (shell->right) {
91: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
92: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
93: in = shell->rightwork;
94: }
95: while (next->next) {
96: if (!next->work) { /* should reuse previous work if the same size */
97: PetscCall(MatCreateVecs(next->mat, NULL, &next->work));
98: }
99: out = next->work;
100: PetscCall(MatMult(next->mat, in, out));
101: in = out;
102: next = next->next;
103: }
104: PetscCall(MatMult(next->mat, in, y));
105: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
106: scale = shell->scale;
107: if (shell->scalings) {
108: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
109: }
110: PetscCall(VecScale(y, scale));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A, Vec x, Vec y)
115: {
116: Mat_Composite *shell = (Mat_Composite *)A->data;
117: Mat_CompositeLink tail = shell->tail;
118: Vec in, out;
119: PetscScalar scale;
120: PetscInt i;
122: PetscFunctionBegin;
123: PetscCheck(tail, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
124: in = x;
125: if (shell->left) {
126: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
127: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
128: in = shell->leftwork;
129: }
130: while (tail->prev) {
131: if (!tail->prev->work) { /* should reuse previous work if the same size */
132: PetscCall(MatCreateVecs(tail->mat, NULL, &tail->prev->work));
133: }
134: out = tail->prev->work;
135: PetscCall(MatMultTranspose(tail->mat, in, out));
136: in = out;
137: tail = tail->prev;
138: }
139: PetscCall(MatMultTranspose(tail->mat, in, y));
140: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
142: scale = shell->scale;
143: if (shell->scalings) {
144: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
145: }
146: PetscCall(VecScale(y, scale));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatMult_Composite(Mat mat, Vec x, Vec y)
151: {
152: Mat_Composite *shell = (Mat_Composite *)mat->data;
153: Mat_CompositeLink cur = shell->head;
154: Vec in, y2, xin;
155: Mat A, B;
156: PetscInt i, j, k, n, nuniq, lo, hi, mid, *gindices, *buf, *tmp, tot;
157: const PetscScalar *vals;
158: const PetscInt *garray;
159: IS ix, iy;
160: PetscBool match;
162: PetscFunctionBegin;
163: PetscCheck(cur, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
164: in = x;
165: if (shell->right) {
166: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
167: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
168: in = shell->rightwork;
169: }
171: /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time
172: we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and
173: it is legal to merge Mvctx, because all component matrices have the same size.
174: */
175: if (shell->merge_mvctx && !shell->Mvctx) {
176: /* Currently only implemented for MATMPIAIJ */
177: for (cur = shell->head; cur; cur = cur->next) {
178: PetscCall(PetscObjectTypeCompare((PetscObject)cur->mat, MATMPIAIJ, &match));
179: if (!match) {
180: shell->merge_mvctx = PETSC_FALSE;
181: goto skip_merge_mvctx;
182: }
183: }
185: /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */
186: tot = 0;
187: for (cur = shell->head; cur; cur = cur->next) {
188: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, NULL));
189: PetscCall(MatGetLocalSize(B, NULL, &n));
190: tot += n;
191: }
192: PetscCall(PetscMalloc3(tot, &shell->location, tot, &shell->larray, shell->nmat, &shell->lvecs));
193: shell->len = tot;
195: /* Go through matrices second time to sort off-diag columns and remove dups */
196: PetscCall(PetscMalloc1(tot, &gindices)); /* No Malloc2() since we will give one to petsc and free the other */
197: PetscCall(PetscMalloc1(tot, &buf));
198: nuniq = 0; /* Number of unique nonzero columns */
199: for (cur = shell->head; cur; cur = cur->next) {
200: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
201: PetscCall(MatGetLocalSize(B, NULL, &n));
202: /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */
203: i = j = k = 0;
204: while (i < n && j < nuniq) {
205: if (garray[i] < gindices[j]) buf[k++] = garray[i++];
206: else if (garray[i] > gindices[j]) buf[k++] = gindices[j++];
207: else {
208: buf[k++] = garray[i++];
209: j++;
210: }
211: }
212: /* Copy leftover in garray[] or gindices[] */
213: if (i < n) {
214: PetscCall(PetscArraycpy(buf + k, garray + i, n - i));
215: nuniq = k + n - i;
216: } else if (j < nuniq) {
217: PetscCall(PetscArraycpy(buf + k, gindices + j, nuniq - j));
218: nuniq = k + nuniq - j;
219: } else nuniq = k;
220: /* Swap gindices and buf to merge garray of the next matrix */
221: tmp = gindices;
222: gindices = buf;
223: buf = tmp;
224: }
225: PetscCall(PetscFree(buf));
227: /* Go through matrices third time to build a map from gindices[] to garray[] */
228: tot = 0;
229: for (cur = shell->head, j = 0; cur; cur = cur->next, j++) { /* j-th matrix */
230: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
231: PetscCall(MatGetLocalSize(B, NULL, &n));
232: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &shell->lvecs[j]));
233: /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */
234: lo = 0;
235: for (i = 0; i < n; i++) {
236: hi = nuniq;
237: while (hi - lo > 1) {
238: mid = lo + (hi - lo) / 2;
239: if (garray[i] < gindices[mid]) hi = mid;
240: else lo = mid;
241: }
242: shell->location[tot + i] = lo; /* gindices[lo] = garray[i] */
243: lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */
244: }
245: tot += n;
246: }
248: /* Build merged Mvctx */
249: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nuniq, gindices, PETSC_OWN_POINTER, &ix));
250: PetscCall(ISCreateStride(PETSC_COMM_SELF, nuniq, 0, 1, &iy));
251: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &xin));
252: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nuniq, &shell->gvec));
253: PetscCall(VecScatterCreate(xin, ix, shell->gvec, iy, &shell->Mvctx));
254: PetscCall(VecDestroy(&xin));
255: PetscCall(ISDestroy(&ix));
256: PetscCall(ISDestroy(&iy));
257: }
259: skip_merge_mvctx:
260: PetscCall(VecSet(y, 0));
261: if (!shell->leftwork2) PetscCall(VecDuplicate(y, &shell->leftwork2));
262: y2 = shell->leftwork2;
264: if (shell->Mvctx) { /* Have a merged Mvctx */
265: /* Suppose we want to compute y = sMx, where s is the scaling factor and A, B are matrix M's diagonal/off-diagonal part. We could do
266: in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an opportunity to
267: overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
268: */
269: PetscCall(VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
270: PetscCall(VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
272: PetscCall(VecGetArrayRead(shell->gvec, &vals));
273: for (i = 0; i < shell->len; i++) shell->larray[i] = vals[shell->location[i]];
274: PetscCall(VecRestoreArrayRead(shell->gvec, &vals));
276: for (cur = shell->head, tot = i = 0; cur; cur = cur->next, i++) { /* i-th matrix */
277: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, &A, &B, NULL));
278: PetscUseTypeMethod(A, mult, in, y2);
279: PetscCall(MatGetLocalSize(B, NULL, &n));
280: PetscCall(VecPlaceArray(shell->lvecs[i], &shell->larray[tot]));
281: PetscCall((*B->ops->multadd)(B, shell->lvecs[i], y2, y2));
282: PetscCall(VecResetArray(shell->lvecs[i]));
283: PetscCall(VecAXPY(y, (shell->scalings ? shell->scalings[i] : 1.0), y2));
284: tot += n;
285: }
286: } else {
287: if (shell->scalings) {
288: for (cur = shell->head, i = 0; cur; cur = cur->next, i++) {
289: PetscCall(MatMult(cur->mat, in, y2));
290: PetscCall(VecAXPY(y, shell->scalings[i], y2));
291: }
292: } else {
293: for (cur = shell->head; cur; cur = cur->next) PetscCall(MatMultAdd(cur->mat, in, y, y));
294: }
295: }
297: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
298: PetscCall(VecScale(y, shell->scale));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: static PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y)
303: {
304: Mat_Composite *shell = (Mat_Composite *)A->data;
305: Mat_CompositeLink next = shell->head;
306: Vec in, y2 = NULL;
307: PetscInt i;
309: PetscFunctionBegin;
310: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
311: in = x;
312: if (shell->left) {
313: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
314: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
315: in = shell->leftwork;
316: }
318: PetscCall(MatMultTranspose(next->mat, in, y));
319: if (shell->scalings) {
320: PetscCall(VecScale(y, shell->scalings[0]));
321: if (!shell->rightwork2) PetscCall(VecDuplicate(y, &shell->rightwork2));
322: y2 = shell->rightwork2;
323: }
324: i = 1;
325: while ((next = next->next)) {
326: if (!shell->scalings) PetscCall(MatMultTransposeAdd(next->mat, in, y, y));
327: else {
328: PetscCall(MatMultTranspose(next->mat, in, y2));
329: PetscCall(VecAXPY(y, shell->scalings[i++], y2));
330: }
331: }
332: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
333: PetscCall(VecScale(y, shell->scale));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z)
338: {
339: Mat_Composite *shell = (Mat_Composite *)A->data;
341: PetscFunctionBegin;
342: if (y != z) {
343: PetscCall(MatMult(A, x, z));
344: PetscCall(VecAXPY(z, 1.0, y));
345: } else {
346: if (!shell->leftwork) PetscCall(VecDuplicate(z, &shell->leftwork));
347: PetscCall(MatMult(A, x, shell->leftwork));
348: PetscCall(VecCopy(y, z));
349: PetscCall(VecAXPY(z, 1.0, shell->leftwork));
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z)
355: {
356: Mat_Composite *shell = (Mat_Composite *)A->data;
358: PetscFunctionBegin;
359: if (y != z) {
360: PetscCall(MatMultTranspose(A, x, z));
361: PetscCall(VecAXPY(z, 1.0, y));
362: } else {
363: if (!shell->rightwork) PetscCall(VecDuplicate(z, &shell->rightwork));
364: PetscCall(MatMultTranspose(A, x, shell->rightwork));
365: PetscCall(VecCopy(y, z));
366: PetscCall(VecAXPY(z, 1.0, shell->rightwork));
367: }
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v)
372: {
373: Mat_Composite *shell = (Mat_Composite *)A->data;
374: Mat_CompositeLink next = shell->head;
375: PetscInt i;
377: PetscFunctionBegin;
378: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
379: PetscCheck(!shell->right && !shell->left, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get diagonal if left or right scaling");
381: PetscCall(MatGetDiagonal(next->mat, v));
382: if (shell->scalings) PetscCall(VecScale(v, shell->scalings[0]));
384: if (next->next && !shell->work) PetscCall(VecDuplicate(v, &shell->work));
385: i = 1;
386: while ((next = next->next)) {
387: PetscCall(MatGetDiagonal(next->mat, shell->work));
388: PetscCall(VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work));
389: }
390: PetscCall(VecScale(v, shell->scale));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t)
395: {
396: Mat_Composite *shell = (Mat_Composite *)Y->data;
398: PetscFunctionBegin;
399: if (shell->merge) PetscCall(MatCompositeMerge(Y));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha)
404: {
405: Mat_Composite *a = (Mat_Composite *)inA->data;
407: PetscFunctionBegin;
408: a->scale *= alpha;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: static PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right)
413: {
414: Mat_Composite *a = (Mat_Composite *)inA->data;
416: PetscFunctionBegin;
417: if (left) {
418: if (!a->left) {
419: PetscCall(VecDuplicate(left, &a->left));
420: PetscCall(VecCopy(left, a->left));
421: } else {
422: PetscCall(VecPointwiseMult(a->left, left, a->left));
423: }
424: }
425: if (right) {
426: if (!a->right) {
427: PetscCall(VecDuplicate(right, &a->right));
428: PetscCall(VecCopy(right, a->right));
429: } else {
430: PetscCall(VecPointwiseMult(a->right, right, a->right));
431: }
432: }
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject)
437: {
438: Mat_Composite *a = (Mat_Composite *)A->data;
440: PetscFunctionBegin;
441: PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options");
442: PetscCall(PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL));
443: PetscCall(PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL));
444: PetscCall(PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL));
445: PetscOptionsHeadEnd();
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*@
450: MatCreateComposite - Creates a matrix as the sum or product of one or more matrices
452: Collective
454: Input Parameters:
455: + comm - MPI communicator
456: . nmat - number of matrices to put in
457: - mats - the matrices
459: Output Parameter:
460: . mat - the matrix
462: Options Database Keys:
463: + -mat_composite_merge - merge in `MatAssemblyEnd()`
464: . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices
465: - -mat_composite_merge_type - set merge direction
467: Level: advanced
469: Note:
470: Alternative construction
471: .vb
472: MatCreate(comm,&mat);
473: MatSetSizes(mat,m,n,M,N);
474: MatSetType(mat,MATCOMPOSITE);
475: MatCompositeAddMat(mat,mats[0]);
476: ....
477: MatCompositeAddMat(mat,mats[nmat-1]);
478: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
479: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
480: .ve
482: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
484: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`,
485: `MATCOMPOSITE`, `MatCompositeType`
486: @*/
487: PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat)
488: {
489: PetscInt m, n, M, N, i;
491: PetscFunctionBegin;
492: PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix");
493: PetscAssertPointer(mat, 4);
495: PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n));
496: PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE));
497: PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N));
498: PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE));
499: PetscCall(MatCreate(comm, mat));
500: PetscCall(MatSetSizes(*mat, m, n, M, N));
501: PetscCall(MatSetType(*mat, MATCOMPOSITE));
502: for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i]));
503: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
504: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat)
509: {
510: Mat_Composite *shell = (Mat_Composite *)mat->data;
511: Mat_CompositeLink ilink, next = shell->head;
512: VecType vtype_mat, vtype_smat;
513: PetscBool match;
515: PetscFunctionBegin;
516: PetscCall(PetscNew(&ilink));
517: ilink->next = NULL;
518: PetscCall(PetscObjectReference((PetscObject)smat));
519: ilink->mat = smat;
521: if (!next) shell->head = ilink;
522: else {
523: while (next->next) next = next->next;
524: next->next = ilink;
525: ilink->prev = next;
526: }
527: shell->tail = ilink;
528: shell->nmat += 1;
530: /* If all of the partial matrices have the same default vector type, then the composite matrix should also have this default type.
531: Otherwise, the default type should be "standard". */
532: PetscCall(MatGetVecType(smat, &vtype_smat));
533: if (shell->nmat == 1) PetscCall(MatSetVecType(mat, vtype_smat));
534: else {
535: PetscCall(MatGetVecType(mat, &vtype_mat));
536: PetscCall(PetscStrcmp(vtype_smat, vtype_mat, &match));
537: if (!match) PetscCall(MatSetVecType(mat, VECSTANDARD));
538: }
540: /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
541: if (shell->scalings) {
542: PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings));
543: shell->scalings[shell->nmat - 1] = 1.0;
544: }
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: /*@
549: MatCompositeAddMat - Add another matrix to a composite matrix.
551: Collective
553: Input Parameters:
554: + mat - the composite matrix
555: - smat - the partial matrix
557: Level: advanced
559: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
560: @*/
561: PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat)
562: {
563: PetscFunctionBegin;
566: PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type)
571: {
572: Mat_Composite *b = (Mat_Composite *)mat->data;
574: PetscFunctionBegin;
575: b->type = type;
576: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
577: mat->ops->getdiagonal = NULL;
578: mat->ops->mult = MatMult_Composite_Multiplicative;
579: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
580: b->merge_mvctx = PETSC_FALSE;
581: } else {
582: mat->ops->getdiagonal = MatGetDiagonal_Composite;
583: mat->ops->mult = MatMult_Composite;
584: mat->ops->multtranspose = MatMultTranspose_Composite;
585: }
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@
590: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product.
592: Logically Collective
594: Input Parameters:
595: + mat - the composite matrix
596: - type - the `MatCompositeType` to use for the matrix
598: Level: advanced
600: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE`,
601: `MatCompositeType`
602: @*/
603: PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type)
604: {
605: PetscFunctionBegin;
608: PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type));
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type)
613: {
614: Mat_Composite *b = (Mat_Composite *)mat->data;
616: PetscFunctionBegin;
617: *type = b->type;
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /*@
622: MatCompositeGetType - Returns type of composite.
624: Not Collective
626: Input Parameter:
627: . mat - the composite matrix
629: Output Parameter:
630: . type - type of composite
632: Level: advanced
634: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType`
635: @*/
636: PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type)
637: {
638: PetscFunctionBegin;
640: PetscAssertPointer(type, 2);
641: PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str)
646: {
647: Mat_Composite *b = (Mat_Composite *)mat->data;
649: PetscFunctionBegin;
650: b->structure = str;
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: /*@
655: MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.
657: Not Collective
659: Input Parameters:
660: + mat - the composite matrix
661: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN`
663: Level: advanced
665: Note:
666: Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix.
668: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE`
669: @*/
670: PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str)
671: {
672: PetscFunctionBegin;
674: PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str));
675: PetscFunctionReturn(PETSC_SUCCESS);
676: }
678: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str)
679: {
680: Mat_Composite *b = (Mat_Composite *)mat->data;
682: PetscFunctionBegin;
683: *str = b->structure;
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@
688: MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.
690: Not Collective
692: Input Parameter:
693: . mat - the composite matrix
695: Output Parameter:
696: . str - structure of the matrices
698: Level: advanced
700: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE`
701: @*/
702: PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str)
703: {
704: PetscFunctionBegin;
706: PetscAssertPointer(str, 2);
707: PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str));
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type)
712: {
713: Mat_Composite *shell = (Mat_Composite *)mat->data;
715: PetscFunctionBegin;
716: shell->mergetype = type;
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`.
723: Logically Collective
725: Input Parameters:
726: + mat - the composite matrix
727: - type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]),
728: `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1])
730: Level: advanced
732: Note:
733: The resulting matrix is the same regardless of the `MatCompositeMergeType`. Only the order of operation is changed.
734: If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
735: otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].
737: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE`
738: @*/
739: PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type)
740: {
741: PetscFunctionBegin;
744: PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type));
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
749: {
750: Mat_Composite *shell = (Mat_Composite *)mat->data;
751: Mat_CompositeLink next = shell->head, prev = shell->tail;
752: Mat tmat, newmat;
753: Vec left, right;
754: PetscScalar scale;
755: PetscInt i;
757: PetscFunctionBegin;
758: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
759: scale = shell->scale;
760: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
761: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
762: i = 0;
763: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
764: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++]));
765: while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure));
766: } else {
767: i = shell->nmat - 1;
768: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
769: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--]));
770: while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure));
771: }
772: } else {
773: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
774: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
775: while ((next = next->next)) {
776: PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
777: PetscCall(MatDestroy(&tmat));
778: tmat = newmat;
779: }
780: } else {
781: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
782: while ((prev = prev->prev)) {
783: PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
784: PetscCall(MatDestroy(&tmat));
785: tmat = newmat;
786: }
787: }
788: if (shell->scalings) {
789: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
790: }
791: }
793: if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left));
794: if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right));
796: PetscCall(MatHeaderReplace(mat, &tmat));
798: PetscCall(MatDiagonalScale(mat, left, right));
799: PetscCall(MatScale(mat, scale));
800: PetscCall(VecDestroy(&left));
801: PetscCall(VecDestroy(&right));
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: /*@
806: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
807: by summing or computing the product of all the matrices inside the composite matrix.
809: Collective
811: Input Parameter:
812: . mat - the composite matrix
814: Options Database Keys:
815: + -mat_composite_merge - merge in `MatAssemblyEnd()`
816: - -mat_composite_merge_type - set merge direction
818: Level: advanced
820: Note:
821: The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix.
823: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE`
824: @*/
825: PetscErrorCode MatCompositeMerge(Mat mat)
826: {
827: PetscFunctionBegin;
829: PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat)
834: {
835: Mat_Composite *shell = (Mat_Composite *)mat->data;
837: PetscFunctionBegin;
838: *nmat = shell->nmat;
839: PetscFunctionReturn(PETSC_SUCCESS);
840: }
842: /*@
843: MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.
845: Not Collective
847: Input Parameter:
848: . mat - the composite matrix
850: Output Parameter:
851: . nmat - number of matrices in the composite matrix
853: Level: advanced
855: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
856: @*/
857: PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat)
858: {
859: PetscFunctionBegin;
861: PetscAssertPointer(nmat, 2);
862: PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat));
863: PetscFunctionReturn(PETSC_SUCCESS);
864: }
866: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai)
867: {
868: Mat_Composite *shell = (Mat_Composite *)mat->data;
869: Mat_CompositeLink ilink;
870: PetscInt k;
872: PetscFunctionBegin;
873: PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat);
874: ilink = shell->head;
875: for (k = 0; k < i; k++) ilink = ilink->next;
876: *Ai = ilink->mat;
877: PetscFunctionReturn(PETSC_SUCCESS);
878: }
880: /*@
881: MatCompositeGetMat - Returns the ith matrix from the composite matrix.
883: Logically Collective
885: Input Parameters:
886: + mat - the composite matrix
887: - i - the number of requested matrix
889: Output Parameter:
890: . Ai - ith matrix in composite
892: Level: advanced
894: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE`
895: @*/
896: PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai)
897: {
898: PetscFunctionBegin;
901: PetscAssertPointer(Ai, 3);
902: PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai));
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: static PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings)
907: {
908: Mat_Composite *shell = (Mat_Composite *)mat->data;
909: PetscInt nmat;
911: PetscFunctionBegin;
912: PetscCall(MatCompositeGetNumberMat(mat, &nmat));
913: if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings));
914: PetscCall(PetscArraycpy(shell->scalings, scalings, nmat));
915: PetscFunctionReturn(PETSC_SUCCESS);
916: }
918: /*@
919: MatCompositeSetScalings - Sets separate scaling factors for component matrices.
921: Logically Collective
923: Input Parameters:
924: + mat - the composite matrix
925: - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)
927: Level: advanced
929: .seealso: [](ch_matrices), `Mat`, `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE`
930: @*/
931: PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings)
932: {
933: PetscFunctionBegin;
935: PetscAssertPointer(scalings, 2);
937: PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings));
938: PetscFunctionReturn(PETSC_SUCCESS);
939: }
941: static struct _MatOps MatOps_Values = {NULL,
942: NULL,
943: NULL,
944: MatMult_Composite,
945: MatMultAdd_Composite,
946: /* 5*/ MatMultTranspose_Composite,
947: MatMultTransposeAdd_Composite,
948: NULL,
949: NULL,
950: NULL,
951: /* 10*/ NULL,
952: NULL,
953: NULL,
954: NULL,
955: NULL,
956: /* 15*/ NULL,
957: NULL,
958: MatGetDiagonal_Composite,
959: MatDiagonalScale_Composite,
960: NULL,
961: /* 20*/ NULL,
962: MatAssemblyEnd_Composite,
963: NULL,
964: NULL,
965: /* 24*/ NULL,
966: NULL,
967: NULL,
968: NULL,
969: NULL,
970: /* 29*/ NULL,
971: NULL,
972: NULL,
973: NULL,
974: NULL,
975: /* 34*/ NULL,
976: NULL,
977: NULL,
978: NULL,
979: NULL,
980: /* 39*/ NULL,
981: NULL,
982: NULL,
983: NULL,
984: NULL,
985: /* 44*/ NULL,
986: MatScale_Composite,
987: MatShift_Basic,
988: NULL,
989: NULL,
990: /* 49*/ NULL,
991: NULL,
992: NULL,
993: NULL,
994: NULL,
995: /* 54*/ NULL,
996: NULL,
997: NULL,
998: NULL,
999: NULL,
1000: /* 59*/ NULL,
1001: MatDestroy_Composite,
1002: NULL,
1003: NULL,
1004: NULL,
1005: /* 64*/ NULL,
1006: NULL,
1007: NULL,
1008: NULL,
1009: NULL,
1010: /* 69*/ NULL,
1011: NULL,
1012: NULL,
1013: NULL,
1014: NULL,
1015: /* 74*/ NULL,
1016: NULL,
1017: MatSetFromOptions_Composite,
1018: NULL,
1019: NULL,
1020: /* 79*/ NULL,
1021: NULL,
1022: NULL,
1023: NULL,
1024: NULL,
1025: /* 84*/ NULL,
1026: NULL,
1027: NULL,
1028: NULL,
1029: NULL,
1030: /* 89*/ NULL,
1031: NULL,
1032: NULL,
1033: NULL,
1034: NULL,
1035: /* 94*/ NULL,
1036: NULL,
1037: NULL,
1038: NULL,
1039: NULL,
1040: /*99*/ NULL,
1041: NULL,
1042: NULL,
1043: NULL,
1044: NULL,
1045: /*104*/ NULL,
1046: NULL,
1047: NULL,
1048: NULL,
1049: NULL,
1050: /*109*/ NULL,
1051: NULL,
1052: NULL,
1053: NULL,
1054: NULL,
1055: /*114*/ NULL,
1056: NULL,
1057: NULL,
1058: NULL,
1059: NULL,
1060: /*119*/ NULL,
1061: NULL,
1062: NULL,
1063: NULL,
1064: NULL,
1065: /*124*/ NULL,
1066: NULL,
1067: NULL,
1068: NULL,
1069: NULL,
1070: /*129*/ NULL,
1071: NULL,
1072: NULL,
1073: NULL,
1074: NULL,
1075: /*134*/ NULL,
1076: NULL,
1077: NULL,
1078: NULL,
1079: NULL,
1080: /*139*/ NULL,
1081: NULL,
1082: NULL,
1083: NULL,
1084: NULL,
1085: /*144*/ NULL,
1086: NULL,
1087: NULL,
1088: NULL,
1089: NULL,
1090: NULL,
1091: /*150*/ NULL,
1092: NULL};
1094: /*MC
1095: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
1096: The matrices need to have a correct size and parallel layout for the sum or product to be valid.
1098: Level: advanced
1100: Note:
1101: To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`);
1103: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`,
1104: `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()`
1105: M*/
1107: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1108: {
1109: Mat_Composite *b;
1111: PetscFunctionBegin;
1112: PetscCall(PetscNew(&b));
1113: A->data = (void *)b;
1114: A->ops[0] = MatOps_Values;
1116: PetscCall(PetscLayoutSetUp(A->rmap));
1117: PetscCall(PetscLayoutSetUp(A->cmap));
1119: A->assembled = PETSC_TRUE;
1120: A->preallocated = PETSC_TRUE;
1121: b->type = MAT_COMPOSITE_ADDITIVE;
1122: b->scale = 1.0;
1123: b->nmat = 0;
1124: b->merge = PETSC_FALSE;
1125: b->mergetype = MAT_COMPOSITE_MERGE_RIGHT;
1126: b->structure = DIFFERENT_NONZERO_PATTERN;
1127: b->merge_mvctx = PETSC_TRUE;
1129: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite));
1130: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite));
1131: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite));
1132: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite));
1133: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite));
1134: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite));
1135: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite));
1136: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite));
1137: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite));
1138: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite));
1140: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE));
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }