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;

513:   PetscFunctionBegin;
514:   PetscCall(PetscNew(&ilink));
515:   ilink->next = NULL;
516:   PetscCall(PetscObjectReference((PetscObject)smat));
517:   ilink->mat = smat;

519:   if (!next) shell->head = ilink;
520:   else {
521:     while (next->next) next = next->next;
522:     next->next  = ilink;
523:     ilink->prev = next;
524:   }
525:   shell->tail = ilink;
526:   shell->nmat += 1;

528:   /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
529:   if (shell->scalings) {
530:     PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings));
531:     shell->scalings[shell->nmat - 1] = 1.0;
532:   }
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: /*@
537:   MatCompositeAddMat - Add another matrix to a composite matrix.

539:   Collective

541:   Input Parameters:
542: + mat  - the composite matrix
543: - smat - the partial matrix

545:   Level: advanced

547: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
548: @*/
549: PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat)
550: {
551:   PetscFunctionBegin;
554:   PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat));
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }

558: static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type)
559: {
560:   Mat_Composite *b = (Mat_Composite *)mat->data;

562:   PetscFunctionBegin;
563:   b->type = type;
564:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
565:     mat->ops->getdiagonal   = NULL;
566:     mat->ops->mult          = MatMult_Composite_Multiplicative;
567:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
568:     b->merge_mvctx          = PETSC_FALSE;
569:   } else {
570:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
571:     mat->ops->mult          = MatMult_Composite;
572:     mat->ops->multtranspose = MatMultTranspose_Composite;
573:   }
574:   PetscFunctionReturn(PETSC_SUCCESS);
575: }

577: /*@
578:   MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product.

580:   Logically Collective

582:   Input Parameters:
583: + mat  - the composite matrix
584: - type - the `MatCompositeType` to use for the matrix

586:   Level: advanced

588: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE`,
589:           `MatCompositeType`
590: @*/
591: PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type)
592: {
593:   PetscFunctionBegin;
596:   PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type));
597:   PetscFunctionReturn(PETSC_SUCCESS);
598: }

600: static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type)
601: {
602:   Mat_Composite *b = (Mat_Composite *)mat->data;

604:   PetscFunctionBegin;
605:   *type = b->type;
606:   PetscFunctionReturn(PETSC_SUCCESS);
607: }

609: /*@
610:   MatCompositeGetType - Returns type of composite.

612:   Not Collective

614:   Input Parameter:
615: . mat - the composite matrix

617:   Output Parameter:
618: . type - type of composite

620:   Level: advanced

622: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType`
623: @*/
624: PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type)
625: {
626:   PetscFunctionBegin;
628:   PetscAssertPointer(type, 2);
629:   PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str)
634: {
635:   Mat_Composite *b = (Mat_Composite *)mat->data;

637:   PetscFunctionBegin;
638:   b->structure = str;
639:   PetscFunctionReturn(PETSC_SUCCESS);
640: }

642: /*@
643:   MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.

645:   Not Collective

647:   Input Parameters:
648: + mat - the composite matrix
649: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN`

651:   Level: advanced

653:   Note:
654:   Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix.

656: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE`
657: @*/
658: PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str)
659: {
660:   PetscFunctionBegin;
662:   PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str));
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str)
667: {
668:   Mat_Composite *b = (Mat_Composite *)mat->data;

670:   PetscFunctionBegin;
671:   *str = b->structure;
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }

675: /*@
676:   MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.

678:   Not Collective

680:   Input Parameter:
681: . mat - the composite matrix

683:   Output Parameter:
684: . str - structure of the matrices

686:   Level: advanced

688: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE`
689: @*/
690: PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str)
691: {
692:   PetscFunctionBegin;
694:   PetscAssertPointer(str, 2);
695:   PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str));
696:   PetscFunctionReturn(PETSC_SUCCESS);
697: }

699: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type)
700: {
701:   Mat_Composite *shell = (Mat_Composite *)mat->data;

703:   PetscFunctionBegin;
704:   shell->mergetype = type;
705:   PetscFunctionReturn(PETSC_SUCCESS);
706: }

708: /*@
709:   MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`.

711:   Logically Collective

713:   Input Parameters:
714: + mat  - the composite matrix
715: - type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]),
716:           `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1])

718:   Level: advanced

720:   Note:
721:   The resulting matrix is the same regardless of the `MatCompositeMergeType`. Only the order of operation is changed.
722:   If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
723:   otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].

725: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE`
726: @*/
727: PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type)
728: {
729:   PetscFunctionBegin;
732:   PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type));
733:   PetscFunctionReturn(PETSC_SUCCESS);
734: }

736: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
737: {
738:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
739:   Mat_CompositeLink next = shell->head, prev = shell->tail;
740:   Mat               tmat, newmat;
741:   Vec               left, right;
742:   PetscScalar       scale;
743:   PetscInt          i;

745:   PetscFunctionBegin;
746:   PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
747:   scale = shell->scale;
748:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
749:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
750:       i = 0;
751:       PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
752:       if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++]));
753:       while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure));
754:     } else {
755:       i = shell->nmat - 1;
756:       PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
757:       if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--]));
758:       while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure));
759:     }
760:   } else {
761:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
762:       PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
763:       while ((next = next->next)) {
764:         PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
765:         PetscCall(MatDestroy(&tmat));
766:         tmat = newmat;
767:       }
768:     } else {
769:       PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
770:       while ((prev = prev->prev)) {
771:         PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
772:         PetscCall(MatDestroy(&tmat));
773:         tmat = newmat;
774:       }
775:     }
776:     if (shell->scalings) {
777:       for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
778:     }
779:   }

781:   if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left));
782:   if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right));

784:   PetscCall(MatHeaderReplace(mat, &tmat));

786:   PetscCall(MatDiagonalScale(mat, left, right));
787:   PetscCall(MatScale(mat, scale));
788:   PetscCall(VecDestroy(&left));
789:   PetscCall(VecDestroy(&right));
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

793: /*@
794:   MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
795:   by summing or computing the product of all the matrices inside the composite matrix.

797:   Collective

799:   Input Parameter:
800: . mat - the composite matrix

802:   Options Database Keys:
803: + -mat_composite_merge      - merge in `MatAssemblyEnd()`
804: - -mat_composite_merge_type - set merge direction

806:   Level: advanced

808:   Note:
809:   The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix.

811: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE`
812: @*/
813: PetscErrorCode MatCompositeMerge(Mat mat)
814: {
815:   PetscFunctionBegin;
817:   PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat));
818:   PetscFunctionReturn(PETSC_SUCCESS);
819: }

821: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat)
822: {
823:   Mat_Composite *shell = (Mat_Composite *)mat->data;

825:   PetscFunctionBegin;
826:   *nmat = shell->nmat;
827:   PetscFunctionReturn(PETSC_SUCCESS);
828: }

830: /*@
831:   MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.

833:   Not Collective

835:   Input Parameter:
836: . mat - the composite matrix

838:   Output Parameter:
839: . nmat - number of matrices in the composite matrix

841:   Level: advanced

843: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
844: @*/
845: PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat)
846: {
847:   PetscFunctionBegin;
849:   PetscAssertPointer(nmat, 2);
850:   PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat));
851:   PetscFunctionReturn(PETSC_SUCCESS);
852: }

854: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai)
855: {
856:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
857:   Mat_CompositeLink ilink;
858:   PetscInt          k;

860:   PetscFunctionBegin;
861:   PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat);
862:   ilink = shell->head;
863:   for (k = 0; k < i; k++) ilink = ilink->next;
864:   *Ai = ilink->mat;
865:   PetscFunctionReturn(PETSC_SUCCESS);
866: }

868: /*@
869:   MatCompositeGetMat - Returns the ith matrix from the composite matrix.

871:   Logically Collective

873:   Input Parameters:
874: + mat - the composite matrix
875: - i   - the number of requested matrix

877:   Output Parameter:
878: . Ai - ith matrix in composite

880:   Level: advanced

882: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE`
883: @*/
884: PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai)
885: {
886:   PetscFunctionBegin;
889:   PetscAssertPointer(Ai, 3);
890:   PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai));
891:   PetscFunctionReturn(PETSC_SUCCESS);
892: }

894: static PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings)
895: {
896:   Mat_Composite *shell = (Mat_Composite *)mat->data;
897:   PetscInt       nmat;

899:   PetscFunctionBegin;
900:   PetscCall(MatCompositeGetNumberMat(mat, &nmat));
901:   if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings));
902:   PetscCall(PetscArraycpy(shell->scalings, scalings, nmat));
903:   PetscFunctionReturn(PETSC_SUCCESS);
904: }

906: /*@
907:   MatCompositeSetScalings - Sets separate scaling factors for component matrices.

909:   Logically Collective

911:   Input Parameters:
912: + mat      - the composite matrix
913: - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)

915:   Level: advanced

917: .seealso: [](ch_matrices), `Mat`, `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE`
918: @*/
919: PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings)
920: {
921:   PetscFunctionBegin;
923:   PetscAssertPointer(scalings, 2);
925:   PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings));
926:   PetscFunctionReturn(PETSC_SUCCESS);
927: }

929: static struct _MatOps MatOps_Values = {NULL,
930:                                        NULL,
931:                                        NULL,
932:                                        MatMult_Composite,
933:                                        MatMultAdd_Composite,
934:                                        /*  5*/ MatMultTranspose_Composite,
935:                                        MatMultTransposeAdd_Composite,
936:                                        NULL,
937:                                        NULL,
938:                                        NULL,
939:                                        /* 10*/ NULL,
940:                                        NULL,
941:                                        NULL,
942:                                        NULL,
943:                                        NULL,
944:                                        /* 15*/ NULL,
945:                                        NULL,
946:                                        MatGetDiagonal_Composite,
947:                                        MatDiagonalScale_Composite,
948:                                        NULL,
949:                                        /* 20*/ NULL,
950:                                        MatAssemblyEnd_Composite,
951:                                        NULL,
952:                                        NULL,
953:                                        /* 24*/ NULL,
954:                                        NULL,
955:                                        NULL,
956:                                        NULL,
957:                                        NULL,
958:                                        /* 29*/ NULL,
959:                                        NULL,
960:                                        NULL,
961:                                        NULL,
962:                                        NULL,
963:                                        /* 34*/ NULL,
964:                                        NULL,
965:                                        NULL,
966:                                        NULL,
967:                                        NULL,
968:                                        /* 39*/ NULL,
969:                                        NULL,
970:                                        NULL,
971:                                        NULL,
972:                                        NULL,
973:                                        /* 44*/ NULL,
974:                                        MatScale_Composite,
975:                                        MatShift_Basic,
976:                                        NULL,
977:                                        NULL,
978:                                        /* 49*/ NULL,
979:                                        NULL,
980:                                        NULL,
981:                                        NULL,
982:                                        NULL,
983:                                        /* 54*/ NULL,
984:                                        NULL,
985:                                        NULL,
986:                                        NULL,
987:                                        NULL,
988:                                        /* 59*/ NULL,
989:                                        MatDestroy_Composite,
990:                                        NULL,
991:                                        NULL,
992:                                        NULL,
993:                                        /* 64*/ NULL,
994:                                        NULL,
995:                                        NULL,
996:                                        NULL,
997:                                        NULL,
998:                                        /* 69*/ NULL,
999:                                        NULL,
1000:                                        NULL,
1001:                                        NULL,
1002:                                        NULL,
1003:                                        /* 74*/ NULL,
1004:                                        NULL,
1005:                                        MatSetFromOptions_Composite,
1006:                                        NULL,
1007:                                        NULL,
1008:                                        /* 79*/ NULL,
1009:                                        NULL,
1010:                                        NULL,
1011:                                        NULL,
1012:                                        NULL,
1013:                                        /* 84*/ NULL,
1014:                                        NULL,
1015:                                        NULL,
1016:                                        NULL,
1017:                                        NULL,
1018:                                        /* 89*/ NULL,
1019:                                        NULL,
1020:                                        NULL,
1021:                                        NULL,
1022:                                        NULL,
1023:                                        /* 94*/ NULL,
1024:                                        NULL,
1025:                                        NULL,
1026:                                        NULL,
1027:                                        NULL,
1028:                                        /*99*/ NULL,
1029:                                        NULL,
1030:                                        NULL,
1031:                                        NULL,
1032:                                        NULL,
1033:                                        /*104*/ NULL,
1034:                                        NULL,
1035:                                        NULL,
1036:                                        NULL,
1037:                                        NULL,
1038:                                        /*109*/ NULL,
1039:                                        NULL,
1040:                                        NULL,
1041:                                        NULL,
1042:                                        NULL,
1043:                                        /*114*/ NULL,
1044:                                        NULL,
1045:                                        NULL,
1046:                                        NULL,
1047:                                        NULL,
1048:                                        /*119*/ NULL,
1049:                                        NULL,
1050:                                        NULL,
1051:                                        NULL,
1052:                                        NULL,
1053:                                        /*124*/ NULL,
1054:                                        NULL,
1055:                                        NULL,
1056:                                        NULL,
1057:                                        NULL,
1058:                                        /*129*/ NULL,
1059:                                        NULL,
1060:                                        NULL,
1061:                                        NULL,
1062:                                        NULL,
1063:                                        /*134*/ NULL,
1064:                                        NULL,
1065:                                        NULL,
1066:                                        NULL,
1067:                                        NULL,
1068:                                        /*139*/ NULL,
1069:                                        NULL,
1070:                                        NULL,
1071:                                        NULL,
1072:                                        NULL,
1073:                                        /*144*/ NULL,
1074:                                        NULL,
1075:                                        NULL,
1076:                                        NULL,
1077:                                        NULL,
1078:                                        NULL,
1079:                                        /*150*/ NULL,
1080:                                        NULL};

1082: /*MC
1083:    MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
1084:     The matrices need to have a correct size and parallel layout for the sum or product to be valid.

1086:   Level: advanced

1088:    Note:
1089:    To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`);

1091: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`,
1092:           `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()`
1093: M*/

1095: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1096: {
1097:   Mat_Composite *b;

1099:   PetscFunctionBegin;
1100:   PetscCall(PetscNew(&b));
1101:   A->data   = (void *)b;
1102:   A->ops[0] = MatOps_Values;

1104:   PetscCall(PetscLayoutSetUp(A->rmap));
1105:   PetscCall(PetscLayoutSetUp(A->cmap));

1107:   A->assembled    = PETSC_TRUE;
1108:   A->preallocated = PETSC_TRUE;
1109:   b->type         = MAT_COMPOSITE_ADDITIVE;
1110:   b->scale        = 1.0;
1111:   b->nmat         = 0;
1112:   b->merge        = PETSC_FALSE;
1113:   b->mergetype    = MAT_COMPOSITE_MERGE_RIGHT;
1114:   b->structure    = DIFFERENT_NONZERO_PATTERN;
1115:   b->merge_mvctx  = PETSC_TRUE;

1117:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite));
1118:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite));
1119:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite));
1120:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite));
1121:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite));
1122:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite));
1123:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite));
1124:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite));
1125:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite));
1126:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite));

1128:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE));
1129:   PetscFunctionReturn(PETSC_SUCCESS);
1130: }