Actual source code: mcomposite.c


  2: #include <petsc/private/matimpl.h>

  4: const char *const MatCompositeMergeTypes[] = {"left", "right", "MatCompositeMergeType", "MAT_COMPOSITE_", NULL};

  6: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
  7: struct _Mat_CompositeLink {
  8:   Mat               mat;
  9:   Vec               work;
 10:   Mat_CompositeLink next, prev;
 11: };

 13: typedef struct {
 14:   MatCompositeType      type;
 15:   Mat_CompositeLink     head, tail;
 16:   Vec                   work;
 17:   PetscScalar           scale;                                      /* scale factor supplied with MatScale() */
 18:   Vec                   left, right;                                /* left and right diagonal scaling provided with MatDiagonalScale() */
 19:   Vec                   leftwork, rightwork, leftwork2, rightwork2; /* Two pairs of working vectors */
 20:   PetscInt              nmat;
 21:   PetscBool             merge;
 22:   MatCompositeMergeType mergetype;
 23:   MatStructure          structure;

 25:   PetscScalar *scalings;
 26:   PetscBool    merge_mvctx; /* Whether need to merge mvctx of component matrices */
 27:   Vec         *lvecs;       /* [nmat] Basically, they are Mvctx->lvec of each component matrix */
 28:   PetscScalar *larray;      /* [len] Data arrays of lvecs[] are stored consecutively in larray */
 29:   PetscInt     len;         /* Length of larray[] */
 30:   Vec          gvec;        /* Union of lvecs[] without duplicated entries */
 31:   PetscInt    *location;    /* A map that maps entries in garray[] to larray[] */
 32:   VecScatter   Mvctx;
 33: } Mat_Composite;

 35: PetscErrorCode MatDestroy_Composite(Mat mat)
 36: {
 37:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
 38:   Mat_CompositeLink next  = shell->head, oldnext;
 39:   PetscInt          i;

 41:   PetscFunctionBegin;
 42:   while (next) {
 43:     PetscCall(MatDestroy(&next->mat));
 44:     if (next->work && (!next->next || next->work != next->next->work)) PetscCall(VecDestroy(&next->work));
 45:     oldnext = next;
 46:     next    = next->next;
 47:     PetscCall(PetscFree(oldnext));
 48:   }
 49:   PetscCall(VecDestroy(&shell->work));
 50:   PetscCall(VecDestroy(&shell->left));
 51:   PetscCall(VecDestroy(&shell->right));
 52:   PetscCall(VecDestroy(&shell->leftwork));
 53:   PetscCall(VecDestroy(&shell->rightwork));
 54:   PetscCall(VecDestroy(&shell->leftwork2));
 55:   PetscCall(VecDestroy(&shell->rightwork2));

 57:   if (shell->Mvctx) {
 58:     for (i = 0; i < shell->nmat; i++) PetscCall(VecDestroy(&shell->lvecs[i]));
 59:     PetscCall(PetscFree3(shell->location, shell->larray, shell->lvecs));
 60:     PetscCall(PetscFree(shell->larray));
 61:     PetscCall(VecDestroy(&shell->gvec));
 62:     PetscCall(VecScatterDestroy(&shell->Mvctx));
 63:   }

 65:   PetscCall(PetscFree(shell->scalings));
 66:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL));
 67:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL));
 68:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL));
 69:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL));
 70:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL));
 71:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL));
 72:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL));
 73:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL));
 74:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL));
 75:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL));
 76:   PetscCall(PetscFree(mat->data));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: PetscErrorCode MatMult_Composite_Multiplicative(Mat A, Vec x, Vec y)
 81: {
 82:   Mat_Composite    *shell = (Mat_Composite *)A->data;
 83:   Mat_CompositeLink next  = shell->head;
 84:   Vec               in, out;
 85:   PetscScalar       scale;
 86:   PetscInt          i;

 88:   PetscFunctionBegin;
 89:   PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
 90:   in = x;
 91:   if (shell->right) {
 92:     if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
 93:     PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
 94:     in = shell->rightwork;
 95:   }
 96:   while (next->next) {
 97:     if (!next->work) { /* should reuse previous work if the same size */
 98:       PetscCall(MatCreateVecs(next->mat, NULL, &next->work));
 99:     }
100:     out = next->work;
101:     PetscCall(MatMult(next->mat, in, out));
102:     in   = out;
103:     next = next->next;
104:   }
105:   PetscCall(MatMult(next->mat, in, y));
106:   if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
107:   scale = shell->scale;
108:   if (shell->scalings) {
109:     for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
110:   }
111:   PetscCall(VecScale(y, scale));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A, Vec x, Vec y)
116: {
117:   Mat_Composite    *shell = (Mat_Composite *)A->data;
118:   Mat_CompositeLink tail  = shell->tail;
119:   Vec               in, out;
120:   PetscScalar       scale;
121:   PetscInt          i;

123:   PetscFunctionBegin;
124:   PetscCheck(tail, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
125:   in = x;
126:   if (shell->left) {
127:     if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
128:     PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
129:     in = shell->leftwork;
130:   }
131:   while (tail->prev) {
132:     if (!tail->prev->work) { /* should reuse previous work if the same size */
133:       PetscCall(MatCreateVecs(tail->mat, NULL, &tail->prev->work));
134:     }
135:     out = tail->prev->work;
136:     PetscCall(MatMultTranspose(tail->mat, in, out));
137:     in   = out;
138:     tail = tail->prev;
139:   }
140:   PetscCall(MatMultTranspose(tail->mat, in, y));
141:   if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));

143:   scale = shell->scale;
144:   if (shell->scalings) {
145:     for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
146:   }
147:   PetscCall(VecScale(y, scale));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: PetscErrorCode MatMult_Composite(Mat mat, Vec x, Vec y)
152: {
153:   Mat_Composite     *shell = (Mat_Composite *)mat->data;
154:   Mat_CompositeLink  cur   = shell->head;
155:   Vec                in, y2, xin;
156:   Mat                A, B;
157:   PetscInt           i, j, k, n, nuniq, lo, hi, mid, *gindices, *buf, *tmp, tot;
158:   const PetscScalar *vals;
159:   const PetscInt    *garray;
160:   IS                 ix, iy;
161:   PetscBool          match;

163:   PetscFunctionBegin;
164:   PetscCheck(cur, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
165:   in = x;
166:   if (shell->right) {
167:     if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
168:     PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
169:     in = shell->rightwork;
170:   }

172:   /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time
173:      we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and
174:      it is legal to merge Mvctx, because all component matrices have the same size.
175:    */
176:   if (shell->merge_mvctx && !shell->Mvctx) {
177:     /* Currently only implemented for MATMPIAIJ */
178:     for (cur = shell->head; cur; cur = cur->next) {
179:       PetscCall(PetscObjectTypeCompare((PetscObject)cur->mat, MATMPIAIJ, &match));
180:       if (!match) {
181:         shell->merge_mvctx = PETSC_FALSE;
182:         goto skip_merge_mvctx;
183:       }
184:     }

186:     /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */
187:     tot = 0;
188:     for (cur = shell->head; cur; cur = cur->next) {
189:       PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, NULL));
190:       PetscCall(MatGetLocalSize(B, NULL, &n));
191:       tot += n;
192:     }
193:     PetscCall(PetscMalloc3(tot, &shell->location, tot, &shell->larray, shell->nmat, &shell->lvecs));
194:     shell->len = tot;

196:     /* Go through matrices second time to sort off-diag columns and remove dups */
197:     PetscCall(PetscMalloc1(tot, &gindices)); /* No Malloc2() since we will give one to petsc and free the other */
198:     PetscCall(PetscMalloc1(tot, &buf));
199:     nuniq = 0; /* Number of unique nonzero columns */
200:     for (cur = shell->head; cur; cur = cur->next) {
201:       PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
202:       PetscCall(MatGetLocalSize(B, NULL, &n));
203:       /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */
204:       i = j = k = 0;
205:       while (i < n && j < nuniq) {
206:         if (garray[i] < gindices[j]) buf[k++] = garray[i++];
207:         else if (garray[i] > gindices[j]) buf[k++] = gindices[j++];
208:         else {
209:           buf[k++] = garray[i++];
210:           j++;
211:         }
212:       }
213:       /* Copy leftover in garray[] or gindices[] */
214:       if (i < n) {
215:         PetscCall(PetscArraycpy(buf + k, garray + i, n - i));
216:         nuniq = k + n - i;
217:       } else if (j < nuniq) {
218:         PetscCall(PetscArraycpy(buf + k, gindices + j, nuniq - j));
219:         nuniq = k + nuniq - j;
220:       } else nuniq = k;
221:       /* Swap gindices and buf to merge garray of the next matrix */
222:       tmp      = gindices;
223:       gindices = buf;
224:       buf      = tmp;
225:     }
226:     PetscCall(PetscFree(buf));

228:     /* Go through matrices third time to build a map from gindices[] to garray[] */
229:     tot = 0;
230:     for (cur = shell->head, j = 0; cur; cur = cur->next, j++) { /* j-th matrix */
231:       PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
232:       PetscCall(MatGetLocalSize(B, NULL, &n));
233:       PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &shell->lvecs[j]));
234:       /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */
235:       lo = 0;
236:       for (i = 0; i < n; i++) {
237:         hi = nuniq;
238:         while (hi - lo > 1) {
239:           mid = lo + (hi - lo) / 2;
240:           if (garray[i] < gindices[mid]) hi = mid;
241:           else lo = mid;
242:         }
243:         shell->location[tot + i] = lo; /* gindices[lo] = garray[i] */
244:         lo++;                          /* Since garray[i+1] > garray[i], we can safely advance lo */
245:       }
246:       tot += n;
247:     }

249:     /* Build merged Mvctx */
250:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nuniq, gindices, PETSC_OWN_POINTER, &ix));
251:     PetscCall(ISCreateStride(PETSC_COMM_SELF, nuniq, 0, 1, &iy));
252:     PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &xin));
253:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, nuniq, &shell->gvec));
254:     PetscCall(VecScatterCreate(xin, ix, shell->gvec, iy, &shell->Mvctx));
255:     PetscCall(VecDestroy(&xin));
256:     PetscCall(ISDestroy(&ix));
257:     PetscCall(ISDestroy(&iy));
258:   }

260: skip_merge_mvctx:
261:   PetscCall(VecSet(y, 0));
262:   if (!shell->leftwork2) PetscCall(VecDuplicate(y, &shell->leftwork2));
263:   y2 = shell->leftwork2;

265:   if (shell->Mvctx) { /* Have a merged Mvctx */
266:     /* 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
267:        in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an opportunity to
268:        overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
269:      */
270:     PetscCall(VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
271:     PetscCall(VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));

273:     PetscCall(VecGetArrayRead(shell->gvec, &vals));
274:     for (i = 0; i < shell->len; i++) shell->larray[i] = vals[shell->location[i]];
275:     PetscCall(VecRestoreArrayRead(shell->gvec, &vals));

277:     for (cur = shell->head, tot = i = 0; cur; cur = cur->next, i++) { /* i-th matrix */
278:       PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, &A, &B, NULL));
279:       PetscUseTypeMethod(A, mult, in, y2);
280:       PetscCall(MatGetLocalSize(B, NULL, &n));
281:       PetscCall(VecPlaceArray(shell->lvecs[i], &shell->larray[tot]));
282:       PetscCall((*B->ops->multadd)(B, shell->lvecs[i], y2, y2));
283:       PetscCall(VecResetArray(shell->lvecs[i]));
284:       PetscCall(VecAXPY(y, (shell->scalings ? shell->scalings[i] : 1.0), y2));
285:       tot += n;
286:     }
287:   } else {
288:     if (shell->scalings) {
289:       for (cur = shell->head, i = 0; cur; cur = cur->next, i++) {
290:         PetscCall(MatMult(cur->mat, in, y2));
291:         PetscCall(VecAXPY(y, shell->scalings[i], y2));
292:       }
293:     } else {
294:       for (cur = shell->head; cur; cur = cur->next) PetscCall(MatMultAdd(cur->mat, in, y, y));
295:     }
296:   }

298:   if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
299:   PetscCall(VecScale(y, shell->scale));
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y)
304: {
305:   Mat_Composite    *shell = (Mat_Composite *)A->data;
306:   Mat_CompositeLink next  = shell->head;
307:   Vec               in, y2 = NULL;
308:   PetscInt          i;

310:   PetscFunctionBegin;
311:   PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
312:   in = x;
313:   if (shell->left) {
314:     if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
315:     PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
316:     in = shell->leftwork;
317:   }

319:   PetscCall(MatMultTranspose(next->mat, in, y));
320:   if (shell->scalings) {
321:     PetscCall(VecScale(y, shell->scalings[0]));
322:     if (!shell->rightwork2) PetscCall(VecDuplicate(y, &shell->rightwork2));
323:     y2 = shell->rightwork2;
324:   }
325:   i = 1;
326:   while ((next = next->next)) {
327:     if (!shell->scalings) PetscCall(MatMultTransposeAdd(next->mat, in, y, y));
328:     else {
329:       PetscCall(MatMultTranspose(next->mat, in, y2));
330:       PetscCall(VecAXPY(y, shell->scalings[i++], y2));
331:     }
332:   }
333:   if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
334:   PetscCall(VecScale(y, shell->scale));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z)
339: {
340:   Mat_Composite *shell = (Mat_Composite *)A->data;

342:   PetscFunctionBegin;
343:   if (y != z) {
344:     PetscCall(MatMult(A, x, z));
345:     PetscCall(VecAXPY(z, 1.0, y));
346:   } else {
347:     if (!shell->leftwork) PetscCall(VecDuplicate(z, &shell->leftwork));
348:     PetscCall(MatMult(A, x, shell->leftwork));
349:     PetscCall(VecCopy(y, z));
350:     PetscCall(VecAXPY(z, 1.0, shell->leftwork));
351:   }
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z)
356: {
357:   Mat_Composite *shell = (Mat_Composite *)A->data;

359:   PetscFunctionBegin;
360:   if (y != z) {
361:     PetscCall(MatMultTranspose(A, x, z));
362:     PetscCall(VecAXPY(z, 1.0, y));
363:   } else {
364:     if (!shell->rightwork) PetscCall(VecDuplicate(z, &shell->rightwork));
365:     PetscCall(MatMultTranspose(A, x, shell->rightwork));
366:     PetscCall(VecCopy(y, z));
367:     PetscCall(VecAXPY(z, 1.0, shell->rightwork));
368:   }
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v)
373: {
374:   Mat_Composite    *shell = (Mat_Composite *)A->data;
375:   Mat_CompositeLink next  = shell->head;
376:   PetscInt          i;

378:   PetscFunctionBegin;
379:   PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
380:   PetscCheck(!shell->right && !shell->left, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get diagonal if left or right scaling");

382:   PetscCall(MatGetDiagonal(next->mat, v));
383:   if (shell->scalings) PetscCall(VecScale(v, shell->scalings[0]));

385:   if (next->next && !shell->work) PetscCall(VecDuplicate(v, &shell->work));
386:   i = 1;
387:   while ((next = next->next)) {
388:     PetscCall(MatGetDiagonal(next->mat, shell->work));
389:     PetscCall(VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work));
390:   }
391:   PetscCall(VecScale(v, shell->scale));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t)
396: {
397:   Mat_Composite *shell = (Mat_Composite *)Y->data;

399:   PetscFunctionBegin;
400:   if (shell->merge) PetscCall(MatCompositeMerge(Y));
401:   PetscFunctionReturn(PETSC_SUCCESS);
402: }

404: PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha)
405: {
406:   Mat_Composite *a = (Mat_Composite *)inA->data;

408:   PetscFunctionBegin;
409:   a->scale *= alpha;
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right)
414: {
415:   Mat_Composite *a = (Mat_Composite *)inA->data;

417:   PetscFunctionBegin;
418:   if (left) {
419:     if (!a->left) {
420:       PetscCall(VecDuplicate(left, &a->left));
421:       PetscCall(VecCopy(left, a->left));
422:     } else {
423:       PetscCall(VecPointwiseMult(a->left, left, a->left));
424:     }
425:   }
426:   if (right) {
427:     if (!a->right) {
428:       PetscCall(VecDuplicate(right, &a->right));
429:       PetscCall(VecCopy(right, a->right));
430:     } else {
431:       PetscCall(VecPointwiseMult(a->right, right, a->right));
432:     }
433:   }
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject)
438: {
439:   Mat_Composite *a = (Mat_Composite *)A->data;

441:   PetscFunctionBegin;
442:   PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options");
443:   PetscCall(PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL));
444:   PetscCall(PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL));
445:   PetscCall(PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL));
446:   PetscOptionsHeadEnd();
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: /*@
451:    MatCreateComposite - Creates a matrix as the sum or product of one or more matrices

453:   Collective

455:    Input Parameters:
456: +  comm - MPI communicator
457: .  nmat - number of matrices to put in
458: -  mats - the matrices

460:    Output Parameter:
461: .  mat - the matrix

463:    Options Database Keys:
464: +  -mat_composite_merge         - merge in `MatAssemblyEnd()`
465: .  -mat_composite_merge_mvctx   - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices
466: -  -mat_composite_merge_type    - set merge direction

468:    Level: advanced

470:    Note:
471:      Alternative construction
472: .vb
473:        MatCreate(comm,&mat);
474:        MatSetSizes(mat,m,n,M,N);
475:        MatSetType(mat,MATCOMPOSITE);
476:        MatCompositeAddMat(mat,mats[0]);
477:        ....
478:        MatCompositeAddMat(mat,mats[nmat-1]);
479:        MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
480:        MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
481: .ve

483:      For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]

485: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`,
486:           `MATCOMPOSITE`, `MatCompositeType`
487: @*/
488: PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat)
489: {
490:   PetscInt m, n, M, N, i;

492:   PetscFunctionBegin;
493:   PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix");

496:   PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n));
497:   PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE));
498:   PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N));
499:   PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE));
500:   PetscCall(MatCreate(comm, mat));
501:   PetscCall(MatSetSizes(*mat, m, n, M, N));
502:   PetscCall(MatSetType(*mat, MATCOMPOSITE));
503:   for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i]));
504:   PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
505:   PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
506:   PetscFunctionReturn(PETSC_SUCCESS);
507: }

509: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat)
510: {
511:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
512:   Mat_CompositeLink ilink, next = shell->head;

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

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

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

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

540:    Collective

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

546:    Level: advanced

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

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

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

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

581:    Logically Collective

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

587:    Level: advanced

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

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

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

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

613:    Not Collective

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

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

621:    Level: advanced

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

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

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

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

646:    Not Collective

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

652:    Level: advanced

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

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

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

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

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

679:    Not Collective

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

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

687:    Level: advanced

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

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

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

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

712:    Logically Collective

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

719:    Level: advanced

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

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

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

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

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

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

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

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

798:   Collective

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

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

807:    Level: advanced

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

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

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

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

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

834:    Not Collective

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

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

842:    Level: advanced

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

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

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

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

872:    Logically Collective

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

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

881:    Level: advanced

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

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

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

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

910:    Logically Collective

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

916:    Level: advanced

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

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

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

1087:   Level: advanced

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

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

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

1100:   PetscFunctionBegin;
1101:   PetscCall(PetscNew(&b));
1102:   A->data = (void *)b;
1103:   PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps)));

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

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

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

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