xref: /petsc/src/mat/impls/composite/mcomposite.c (revision d7f81bf2c76b9fc3b3edf4bce4e5de7778a41349)
1 
2 #include <petsc/private/matimpl.h>        /*I "petscmat.h" I*/
3 
4 typedef struct _Mat_CompositeLink *Mat_CompositeLink;
5 struct _Mat_CompositeLink {
6   Mat               mat;
7   Vec               work;
8   Mat_CompositeLink next,prev;
9 };
10 
11 typedef struct {
12   MatCompositeType  type;
13   Mat_CompositeLink head,tail;
14   Vec               work;
15   PetscScalar       scale;        /* scale factor supplied with MatScale() */
16   Vec               left,right;   /* left and right diagonal scaling provided with MatDiagonalScale() */
17   Vec               leftwork,rightwork;
18   PetscInt          nmat;
19 } Mat_Composite;
20 
21 PetscErrorCode MatDestroy_Composite(Mat mat)
22 {
23   PetscErrorCode    ierr;
24   Mat_Composite     *shell = (Mat_Composite*)mat->data;
25   Mat_CompositeLink next   = shell->head,oldnext;
26 
27   PetscFunctionBegin;
28   while (next) {
29     ierr = MatDestroy(&next->mat);CHKERRQ(ierr);
30     if (next->work && (!next->next || next->work != next->next->work)) {
31       ierr = VecDestroy(&next->work);CHKERRQ(ierr);
32     }
33     oldnext = next;
34     next    = next->next;
35     ierr    = PetscFree(oldnext);CHKERRQ(ierr);
36   }
37   ierr = VecDestroy(&shell->work);CHKERRQ(ierr);
38   ierr = VecDestroy(&shell->left);CHKERRQ(ierr);
39   ierr = VecDestroy(&shell->right);CHKERRQ(ierr);
40   ierr = VecDestroy(&shell->leftwork);CHKERRQ(ierr);
41   ierr = VecDestroy(&shell->rightwork);CHKERRQ(ierr);
42   ierr = PetscFree(mat->data);CHKERRQ(ierr);
43   PetscFunctionReturn(0);
44 }
45 
46 PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
47 {
48   Mat_Composite     *shell = (Mat_Composite*)A->data;
49   Mat_CompositeLink next   = shell->head;
50   PetscErrorCode    ierr;
51   Vec               in,out;
52 
53   PetscFunctionBegin;
54   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
55   in = x;
56   if (shell->right) {
57     if (!shell->rightwork) {
58       ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr);
59     }
60     ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr);
61     in   = shell->rightwork;
62   }
63   while (next->next) {
64     if (!next->work) { /* should reuse previous work if the same size */
65       ierr = MatCreateVecs(next->mat,NULL,&next->work);CHKERRQ(ierr);
66     }
67     out  = next->work;
68     ierr = MatMult(next->mat,in,out);CHKERRQ(ierr);
69     in   = out;
70     next = next->next;
71   }
72   ierr = MatMult(next->mat,in,y);CHKERRQ(ierr);
73   if (shell->left) {
74     ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr);
75   }
76   ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
77   PetscFunctionReturn(0);
78 }
79 
80 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
81 {
82   Mat_Composite     *shell = (Mat_Composite*)A->data;
83   Mat_CompositeLink tail   = shell->tail;
84   PetscErrorCode    ierr;
85   Vec               in,out;
86 
87   PetscFunctionBegin;
88   if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
89   in = x;
90   if (shell->left) {
91     if (!shell->leftwork) {
92       ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr);
93     }
94     ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr);
95     in   = shell->leftwork;
96   }
97   while (tail->prev) {
98     if (!tail->prev->work) { /* should reuse previous work if the same size */
99       ierr = MatCreateVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr);
100     }
101     out  = tail->prev->work;
102     ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr);
103     in   = out;
104     tail = tail->prev;
105   }
106   ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr);
107   if (shell->right) {
108     ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr);
109   }
110   ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
111   PetscFunctionReturn(0);
112 }
113 
114 PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
115 {
116   Mat_Composite     *shell = (Mat_Composite*)A->data;
117   Mat_CompositeLink next   = shell->head;
118   PetscErrorCode    ierr;
119   Vec               in;
120 
121   PetscFunctionBegin;
122   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
123   in = x;
124   if (shell->right) {
125     if (!shell->rightwork) {
126       ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr);
127     }
128     ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr);
129     in   = shell->rightwork;
130   }
131   ierr = MatMult(next->mat,in,y);CHKERRQ(ierr);
132   while ((next = next->next)) {
133     ierr = MatMultAdd(next->mat,in,y,y);CHKERRQ(ierr);
134   }
135   if (shell->left) {
136     ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr);
137   }
138   ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
139   PetscFunctionReturn(0);
140 }
141 
142 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
143 {
144   Mat_Composite     *shell = (Mat_Composite*)A->data;
145   Mat_CompositeLink next   = shell->head;
146   PetscErrorCode    ierr;
147   Vec               in;
148 
149   PetscFunctionBegin;
150   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
151   in = x;
152   if (shell->left) {
153     if (!shell->leftwork) {
154       ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr);
155     }
156     ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr);
157     in   = shell->leftwork;
158   }
159   ierr = MatMultTranspose(next->mat,in,y);CHKERRQ(ierr);
160   while ((next = next->next)) {
161     ierr = MatMultTransposeAdd(next->mat,in,y,y);CHKERRQ(ierr);
162   }
163   if (shell->right) {
164     ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr);
165   }
166   ierr = VecScale(y,shell->scale);CHKERRQ(ierr);
167   PetscFunctionReturn(0);
168 }
169 
170 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
171 {
172   Mat_Composite     *shell = (Mat_Composite*)A->data;
173   Mat_CompositeLink next   = shell->head;
174   PetscErrorCode    ierr;
175 
176   PetscFunctionBegin;
177   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
178   if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling");
179 
180   ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr);
181   if (next->next && !shell->work) {
182     ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr);
183   }
184   while ((next = next->next)) {
185     ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr);
186     ierr = VecAXPY(v,1.0,shell->work);CHKERRQ(ierr);
187   }
188   ierr = VecScale(v,shell->scale);CHKERRQ(ierr);
189   PetscFunctionReturn(0);
190 }
191 
192 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
193 {
194   PetscErrorCode ierr;
195   PetscBool      flg = PETSC_FALSE;
196 
197   PetscFunctionBegin;
198   ierr = PetscOptionsGetBool(((PetscObject)Y)->options,((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);CHKERRQ(ierr);
199   if (flg) {
200     ierr = MatCompositeMerge(Y);CHKERRQ(ierr);
201   }
202   PetscFunctionReturn(0);
203 }
204 
205 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
206 {
207   Mat_Composite *a = (Mat_Composite*)inA->data;
208 
209   PetscFunctionBegin;
210   a->scale *= alpha;
211   PetscFunctionReturn(0);
212 }
213 
214 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
215 {
216   Mat_Composite  *a = (Mat_Composite*)inA->data;
217   PetscErrorCode ierr;
218 
219   PetscFunctionBegin;
220   if (left) {
221     if (!a->left) {
222       ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr);
223       ierr = VecCopy(left,a->left);CHKERRQ(ierr);
224     } else {
225       ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr);
226     }
227   }
228   if (right) {
229     if (!a->right) {
230       ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr);
231       ierr = VecCopy(right,a->right);CHKERRQ(ierr);
232     } else {
233       ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr);
234     }
235   }
236   PetscFunctionReturn(0);
237 }
238 
239 static struct _MatOps MatOps_Values = {0,
240                                        0,
241                                        0,
242                                        MatMult_Composite,
243                                        0,
244                                 /*  5*/ MatMultTranspose_Composite,
245                                        0,
246                                        0,
247                                        0,
248                                        0,
249                                 /* 10*/ 0,
250                                        0,
251                                        0,
252                                        0,
253                                        0,
254                                 /* 15*/ 0,
255                                        0,
256                                        MatGetDiagonal_Composite,
257                                        MatDiagonalScale_Composite,
258                                        0,
259                                 /* 20*/ 0,
260                                        MatAssemblyEnd_Composite,
261                                        0,
262                                        0,
263                                /* 24*/ 0,
264                                        0,
265                                        0,
266                                        0,
267                                        0,
268                                /* 29*/ 0,
269                                        0,
270                                        0,
271                                        0,
272                                        0,
273                                /* 34*/ 0,
274                                        0,
275                                        0,
276                                        0,
277                                        0,
278                                /* 39*/ 0,
279                                        0,
280                                        0,
281                                        0,
282                                        0,
283                                /* 44*/ 0,
284                                        MatScale_Composite,
285                                        MatShift_Basic,
286                                        0,
287                                        0,
288                                /* 49*/ 0,
289                                        0,
290                                        0,
291                                        0,
292                                        0,
293                                /* 54*/ 0,
294                                        0,
295                                        0,
296                                        0,
297                                        0,
298                                /* 59*/ 0,
299                                        MatDestroy_Composite,
300                                        0,
301                                        0,
302                                        0,
303                                /* 64*/ 0,
304                                        0,
305                                        0,
306                                        0,
307                                        0,
308                                /* 69*/ 0,
309                                        0,
310                                        0,
311                                        0,
312                                        0,
313                                /* 74*/ 0,
314                                        0,
315                                        0,
316                                        0,
317                                        0,
318                                /* 79*/ 0,
319                                        0,
320                                        0,
321                                        0,
322                                        0,
323                                /* 84*/ 0,
324                                        0,
325                                        0,
326                                        0,
327                                        0,
328                                /* 89*/ 0,
329                                        0,
330                                        0,
331                                        0,
332                                        0,
333                                /* 94*/ 0,
334                                        0,
335                                        0,
336                                        0,
337                                        0,
338                                 /*99*/ 0,
339                                        0,
340                                        0,
341                                        0,
342                                        0,
343                                /*104*/ 0,
344                                        0,
345                                        0,
346                                        0,
347                                        0,
348                                /*109*/ 0,
349                                        0,
350                                        0,
351                                        0,
352                                        0,
353                                /*114*/ 0,
354                                        0,
355                                        0,
356                                        0,
357                                        0,
358                                /*119*/ 0,
359                                        0,
360                                        0,
361                                        0,
362                                        0,
363                                /*124*/ 0,
364                                        0,
365                                        0,
366                                        0,
367                                        0,
368                                /*129*/ 0,
369                                        0,
370                                        0,
371                                        0,
372                                        0,
373                                /*134*/ 0,
374                                        0,
375                                        0,
376                                        0,
377                                        0,
378                                /*139*/ 0,
379                                        0,
380                                        0
381 };
382 
383 /*MC
384    MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
385     The matrices need to have a correct size and parallel layout for the sum or product to be valid.
386 
387    Notes:
388     to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
389 
390   Level: advanced
391 
392 .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
393 M*/
394 
395 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
396 {
397   Mat_Composite  *b;
398   PetscErrorCode ierr;
399 
400   PetscFunctionBegin;
401   ierr    = PetscNewLog(A,&b);CHKERRQ(ierr);
402   A->data = (void*)b;
403   ierr    = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr);
404 
405   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
406   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
407 
408   A->assembled    = PETSC_TRUE;
409   A->preallocated = PETSC_TRUE;
410   b->type         = MAT_COMPOSITE_ADDITIVE;
411   b->scale        = 1.0;
412   b->nmat         = 0;
413   ierr            = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr);
414 
415   ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr);
416   ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr);
417   ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr);
418   ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNMat_C",MatCompositeGetNMat_Composite);CHKERRQ(ierr);
419   ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr);
420 
421   ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr);
422   PetscFunctionReturn(0);
423 }
424 
425 /*@
426    MatCreateComposite - Creates a matrix as the sum of zero or more matrices
427 
428   Collective on MPI_Comm
429 
430    Input Parameters:
431 +  comm - MPI communicator
432 .  nmat - number of matrices to put in
433 -  mats - the matrices
434 
435    Output Parameter:
436 .  mat - the matrix
437 
438    Level: advanced
439 
440    Notes:
441      Alternative construction
442 $       MatCreate(comm,&mat);
443 $       MatSetSizes(mat,m,n,M,N);
444 $       MatSetType(mat,MATCOMPOSITE);
445 $       MatCompositeAddMat(mat,mats[0]);
446 $       ....
447 $       MatCompositeAddMat(mat,mats[nmat-1]);
448 $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
449 $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
450 
451      For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
452 
453 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
454 
455 @*/
456 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
457 {
458   PetscErrorCode ierr;
459   PetscInt       m,n,M,N,i;
460 
461   PetscFunctionBegin;
462   if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
463   PetscValidPointer(mat,3);
464 
465   ierr = MatGetLocalSize(mats[0],PETSC_IGNORE,&n);CHKERRQ(ierr);
466   ierr = MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);CHKERRQ(ierr);
467   ierr = MatGetSize(mats[0],PETSC_IGNORE,&N);CHKERRQ(ierr);
468   ierr = MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);CHKERRQ(ierr);
469   ierr = MatCreate(comm,mat);CHKERRQ(ierr);
470   ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr);
471   ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr);
472   for (i=0; i<nmat; i++) {
473     ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr);
474   }
475   ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
476   ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
477   PetscFunctionReturn(0);
478 }
479 
480 
481 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat)
482 {
483   Mat_Composite     *shell = (Mat_Composite*)mat->data;
484   Mat_CompositeLink ilink,next = shell->head;
485   PetscErrorCode    ierr;
486 
487   PetscFunctionBegin;
488   ierr        = PetscNewLog(mat,&ilink);CHKERRQ(ierr);
489   ilink->next = 0;
490   ierr        = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr);
491   ilink->mat  = smat;
492 
493   if (!next) shell->head = ilink;
494   else {
495     while (next->next) {
496       next = next->next;
497     }
498     next->next  = ilink;
499     ilink->prev = next;
500   }
501   shell->tail =  ilink;
502   shell->nmat += 1;
503   PetscFunctionReturn(0);
504 }
505 
506 /*@
507     MatCompositeAddMat - add another matrix to a composite matrix
508 
509    Collective on Mat
510 
511     Input Parameters:
512 +   mat - the composite matrix
513 -   smat - the partial matrix
514 
515    Level: advanced
516 
517 .seealso: MatCreateComposite()
518 @*/
519 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
520 {
521   PetscErrorCode    ierr;
522 
523   PetscFunctionBegin;
524   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
525   PetscValidHeaderSpecific(smat,MAT_CLASSID,2);
526   ierr = PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));CHKERRQ(ierr);
527   PetscFunctionReturn(0);
528 }
529 
530 static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type)
531 {
532   Mat_Composite  *b = (Mat_Composite*)mat->data;
533 
534   PetscFunctionBegin;
535   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
536     mat->ops->getdiagonal   = 0;
537     mat->ops->mult          = MatMult_Composite_Multiplicative;
538     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
539     b->type                 = MAT_COMPOSITE_MULTIPLICATIVE;
540   } else {
541     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
542     mat->ops->mult          = MatMult_Composite;
543     mat->ops->multtranspose = MatMultTranspose_Composite;
544     b->type                 = MAT_COMPOSITE_ADDITIVE;
545   }
546   PetscFunctionReturn(0);
547 }
548 
549 /*@
550    MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
551 
552   Collective on MPI_Comm
553 
554    Input Parameters:
555 .  mat - the composite matrix
556 
557 
558    Level: advanced
559 
560    Notes:
561       The MatType of the resulting matrix will be the same as the MatType of the FIRST
562     matrix in the composite matrix.
563 
564 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
565 
566 @*/
567 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
568 {
569   PetscErrorCode ierr;
570 
571   PetscFunctionBegin;
572   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
573   ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr);
574   PetscFunctionReturn(0);
575 }
576 
577 static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
578 {
579   Mat_Composite     *shell = (Mat_Composite*)mat->data;
580   Mat_CompositeLink next   = shell->head, prev = shell->tail;
581   PetscErrorCode    ierr;
582   Mat               tmat,newmat;
583   Vec               left,right;
584   PetscScalar       scale;
585 
586   PetscFunctionBegin;
587   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
588 
589   PetscFunctionBegin;
590   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
591     ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr);
592     while ((next = next->next)) {
593       ierr = MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
594     }
595   } else {
596     ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr);
597     while ((next = next->next)) {
598       ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr);
599       ierr = MatDestroy(&tmat);CHKERRQ(ierr);
600       tmat = newmat;
601     }
602   }
603 
604   scale = shell->scale;
605   if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);}
606   if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);}
607 
608   ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr);
609 
610   ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr);
611   ierr = MatScale(mat,scale);CHKERRQ(ierr);
612   ierr = VecDestroy(&left);CHKERRQ(ierr);
613   ierr = VecDestroy(&right);CHKERRQ(ierr);
614   PetscFunctionReturn(0);
615 }
616 
617 /*@
618    MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
619      by summing all the matrices inside the composite matrix.
620 
621   Collective on MPI_Comm
622 
623    Input Parameters:
624 .  mat - the composite matrix
625 
626 
627    Options Database:
628 .  -mat_composite_merge  (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
629 
630    Level: advanced
631 
632    Notes:
633       The MatType of the resulting matrix will be the same as the MatType of the FIRST
634     matrix in the composite matrix.
635 
636 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
637 
638 @*/
639 PetscErrorCode MatCompositeMerge(Mat mat)
640 {
641   PetscErrorCode ierr;
642 
643   PetscFunctionBegin;
644   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
645   ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr);
646   PetscFunctionReturn(0);
647 }
648 
649 static PetscErrorCode MatCompositeGetNMat_Composite(Mat mat,PetscInt *nmat)
650 {
651   Mat_Composite     *shell = (Mat_Composite*)mat->data;
652 
653   PetscFunctionBegin;
654   *nmat = shell->nmat;
655   PetscFunctionReturn(0);
656 }
657 
658 /*@
659    MatCompositeGetNMat - Returns the number of matrices in composite.
660 
661    Not Collective
662 
663    Input Parameter:
664 .  mat - the composite matrix
665 
666    Output Parameter:
667 .  size - the local size
668 .  nmat - number of matrices in composite
669 
670    Level: advanced
671 
672 .seealso: MatCreateComposite(), MatCompositeGetMat()
673 
674 @*/
675 PetscErrorCode MatCompositeGetNMat(Mat mat,PetscInt *nmat)
676 {
677   PetscErrorCode ierr;
678 
679   PetscFunctionBegin;
680   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
681   PetscValidPointer(nmat,2);
682   ierr = PetscUseMethod(mat,"MatCompositeGetNMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr);
683   PetscFunctionReturn(0);
684 }
685 
686 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai)
687 {
688   Mat_Composite     *shell = (Mat_Composite*)mat->data;
689   Mat_CompositeLink ilink;
690   PetscInt          k;
691 
692   PetscFunctionBegin;
693   if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat);
694   ilink = shell->head;
695   for (k=0; k<i; k++) {
696     ilink = ilink->next;
697   }
698   *Ai = ilink->mat;
699   PetscFunctionReturn(0);
700 }
701 
702 /*@
703    MatCompositeGetMat - Returns the ith matrix from composite.
704 
705    Logically Collective on Mat
706 
707    Input Parameter:
708 +  mat - the composite matrix
709 -  i - the number of requested matrix
710 
711    Output Parameter:
712 .  Ai - ith matrix in composite
713 
714    Level: advanced
715 
716 .seealso: MatCreateComposite(), MatCompositeGetNMat()
717 
718 @*/
719 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai)
720 {
721   PetscErrorCode ierr;
722 
723   PetscFunctionBegin;
724   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
725   PetscValidLogicalCollectiveInt(mat,i,2);
726   PetscValidPointer(Ai,3);
727   ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr);
728   PetscFunctionReturn(0);
729 }
730 
731