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