Actual source code: mcomposite.c

petsc-master 2019-06-22
Report Typos and Errors

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

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

  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;
 20:   PetscInt              nmat;
 21:   PetscBool             merge;
 22:   MatCompositeMergeType mergetype;
 23:   MatStructure          structure;
 24: } Mat_Composite;

 26: PetscErrorCode MatDestroy_Composite(Mat mat)
 27: {
 28:   PetscErrorCode    ierr;
 29:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
 30:   Mat_CompositeLink next   = shell->head,oldnext;

 33:   while (next) {
 34:     MatDestroy(&next->mat);
 35:     if (next->work && (!next->next || next->work != next->next->work)) {
 36:       VecDestroy(&next->work);
 37:     }
 38:     oldnext = next;
 39:     next    = next->next;
 40:     PetscFree(oldnext);
 41:   }
 42:   VecDestroy(&shell->work);
 43:   VecDestroy(&shell->left);
 44:   VecDestroy(&shell->right);
 45:   VecDestroy(&shell->leftwork);
 46:   VecDestroy(&shell->rightwork);
 47:   PetscFree(mat->data);
 48:   return(0);
 49: }

 51: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
 52: {
 53:   Mat_Composite     *shell = (Mat_Composite*)A->data;
 54:   Mat_CompositeLink next   = shell->head;
 55:   PetscErrorCode    ierr;
 56:   Vec               in,out;

 59:   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
 60:   in = x;
 61:   if (shell->right) {
 62:     if (!shell->rightwork) {
 63:       VecDuplicate(shell->right,&shell->rightwork);
 64:     }
 65:     VecPointwiseMult(shell->rightwork,shell->right,in);
 66:     in   = shell->rightwork;
 67:   }
 68:   while (next->next) {
 69:     if (!next->work) { /* should reuse previous work if the same size */
 70:       MatCreateVecs(next->mat,NULL,&next->work);
 71:     }
 72:     out  = next->work;
 73:     MatMult(next->mat,in,out);
 74:     in   = out;
 75:     next = next->next;
 76:   }
 77:   MatMult(next->mat,in,y);
 78:   if (shell->left) {
 79:     VecPointwiseMult(y,shell->left,y);
 80:   }
 81:   VecScale(y,shell->scale);
 82:   return(0);
 83: }

 85: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
 86: {
 87:   Mat_Composite     *shell = (Mat_Composite*)A->data;
 88:   Mat_CompositeLink tail   = shell->tail;
 89:   PetscErrorCode    ierr;
 90:   Vec               in,out;

 93:   if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
 94:   in = x;
 95:   if (shell->left) {
 96:     if (!shell->leftwork) {
 97:       VecDuplicate(shell->left,&shell->leftwork);
 98:     }
 99:     VecPointwiseMult(shell->leftwork,shell->left,in);
100:     in   = shell->leftwork;
101:   }
102:   while (tail->prev) {
103:     if (!tail->prev->work) { /* should reuse previous work if the same size */
104:       MatCreateVecs(tail->mat,NULL,&tail->prev->work);
105:     }
106:     out  = tail->prev->work;
107:     MatMultTranspose(tail->mat,in,out);
108:     in   = out;
109:     tail = tail->prev;
110:   }
111:   MatMultTranspose(tail->mat,in,y);
112:   if (shell->right) {
113:     VecPointwiseMult(y,shell->right,y);
114:   }
115:   VecScale(y,shell->scale);
116:   return(0);
117: }

119: PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
120: {
121:   Mat_Composite     *shell = (Mat_Composite*)A->data;
122:   Mat_CompositeLink next   = shell->head;
123:   PetscErrorCode    ierr;
124:   Vec               in;

127:   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
128:   in = x;
129:   if (shell->right) {
130:     if (!shell->rightwork) {
131:       VecDuplicate(shell->right,&shell->rightwork);
132:     }
133:     VecPointwiseMult(shell->rightwork,shell->right,in);
134:     in   = shell->rightwork;
135:   }
136:   MatMult(next->mat,in,y);
137:   while ((next = next->next)) {
138:     MatMultAdd(next->mat,in,y,y);
139:   }
140:   if (shell->left) {
141:     VecPointwiseMult(y,shell->left,y);
142:   }
143:   VecScale(y,shell->scale);
144:   return(0);
145: }

147: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
148: {
149:   Mat_Composite     *shell = (Mat_Composite*)A->data;
150:   Mat_CompositeLink next   = shell->head;
151:   PetscErrorCode    ierr;
152:   Vec               in;

155:   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
156:   in = x;
157:   if (shell->left) {
158:     if (!shell->leftwork) {
159:       VecDuplicate(shell->left,&shell->leftwork);
160:     }
161:     VecPointwiseMult(shell->leftwork,shell->left,in);
162:     in   = shell->leftwork;
163:   }
164:   MatMultTranspose(next->mat,in,y);
165:   while ((next = next->next)) {
166:     MatMultTransposeAdd(next->mat,in,y,y);
167:   }
168:   if (shell->right) {
169:     VecPointwiseMult(y,shell->right,y);
170:   }
171:   VecScale(y,shell->scale);
172:   return(0);
173: }

175: PetscErrorCode MatMultAdd_Composite(Mat A,Vec x,Vec y,Vec z)
176: {
177:   Mat_Composite     *shell = (Mat_Composite*)A->data;
178:   PetscErrorCode    ierr;

181:   if (y != z) {
182:     MatMult(A,x,z);
183:     VecAXPY(z,1.0,y);
184:   } else {
185:     if (!shell->leftwork) {
186:       VecDuplicate(z,&shell->leftwork);
187:     }
188:     MatMult(A,x,shell->leftwork);
189:     VecCopy(y,z);
190:     VecAXPY(z,1.0,shell->leftwork);
191:   }
192:   return(0);
193: }

195: PetscErrorCode MatMultTransposeAdd_Composite(Mat A,Vec x,Vec y, Vec z)
196: {
197:   Mat_Composite     *shell = (Mat_Composite*)A->data;
198:   PetscErrorCode    ierr;

201:   if (y != z) {
202:     MatMultTranspose(A,x,z);
203:     VecAXPY(z,1.0,y);
204:   } else {
205:     if (!shell->rightwork) {
206:       VecDuplicate(z,&shell->rightwork);
207:     }
208:     MatMultTranspose(A,x,shell->rightwork);
209:     VecCopy(y,z);
210:     VecAXPY(z,1.0,shell->rightwork);
211:   }
212:   return(0);
213: }

215: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
216: {
217:   Mat_Composite     *shell = (Mat_Composite*)A->data;
218:   Mat_CompositeLink next   = shell->head;
219:   PetscErrorCode    ierr;

222:   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
223:   if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling");

225:   MatGetDiagonal(next->mat,v);
226:   if (next->next && !shell->work) {
227:     VecDuplicate(v,&shell->work);
228:   }
229:   while ((next = next->next)) {
230:     MatGetDiagonal(next->mat,shell->work);
231:     VecAXPY(v,1.0,shell->work);
232:   }
233:   VecScale(v,shell->scale);
234:   return(0);
235: }

237: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
238: {
239:   Mat_Composite     *shell = (Mat_Composite*)Y->data;
240:   PetscErrorCode    ierr;

243:   if (shell->merge) {
244:     MatCompositeMerge(Y);
245:   }
246:   return(0);
247: }

249: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
250: {
251:   Mat_Composite *a = (Mat_Composite*)inA->data;

254:   a->scale *= alpha;
255:   return(0);
256: }

258: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
259: {
260:   Mat_Composite  *a = (Mat_Composite*)inA->data;

264:   if (left) {
265:     if (!a->left) {
266:       VecDuplicate(left,&a->left);
267:       VecCopy(left,a->left);
268:     } else {
269:       VecPointwiseMult(a->left,left,a->left);
270:     }
271:   }
272:   if (right) {
273:     if (!a->right) {
274:       VecDuplicate(right,&a->right);
275:       VecCopy(right,a->right);
276:     } else {
277:       VecPointwiseMult(a->right,right,a->right);
278:     }
279:   }
280:   return(0);
281: }

283: PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A)
284: {
285:   Mat_Composite  *a = (Mat_Composite*)A->data;

289:   PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");
290:   PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);
291:   PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);
292:   PetscOptionsTail();
293:   return(0);
294: }

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

299:   Collective

301:    Input Parameters:
302: +  comm - MPI communicator
303: .  nmat - number of matrices to put in
304: -  mats - the matrices

306:    Output Parameter:
307: .  mat - the matrix

309:    Options Database Keys:
310: +  -mat_composite_merge - merge in MatAssemblyEnd()
311: -  -mat_composite_merge_type - set merge direction

313:    Level: advanced

315:    Notes:
316:      Alternative construction
317: $       MatCreate(comm,&mat);
318: $       MatSetSizes(mat,m,n,M,N);
319: $       MatSetType(mat,MATCOMPOSITE);
320: $       MatCompositeAddMat(mat,mats[0]);
321: $       ....
322: $       MatCompositeAddMat(mat,mats[nmat-1]);
323: $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
324: $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

328: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeGetMat(), MatCompositeMerge(), MatCompositeSetType(), MATCOMPOSITE

330: @*/
331: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
332: {
334:   PetscInt       m,n,M,N,i;

337:   if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");

340:   MatGetLocalSize(mats[0],PETSC_IGNORE,&n);
341:   MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);
342:   MatGetSize(mats[0],PETSC_IGNORE,&N);
343:   MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);
344:   MatCreate(comm,mat);
345:   MatSetSizes(*mat,m,n,M,N);
346:   MatSetType(*mat,MATCOMPOSITE);
347:   for (i=0; i<nmat; i++) {
348:     MatCompositeAddMat(*mat,mats[i]);
349:   }
350:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
351:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
352:   return(0);
353: }


356: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat)
357: {
358:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
359:   Mat_CompositeLink ilink,next = shell->head;
360:   PetscErrorCode    ierr;

363:   PetscNewLog(mat,&ilink);
364:   ilink->next = 0;
365:   PetscObjectReference((PetscObject)smat);
366:   ilink->mat  = smat;

368:   if (!next) shell->head = ilink;
369:   else {
370:     while (next->next) {
371:       next = next->next;
372:     }
373:     next->next  = ilink;
374:     ilink->prev = next;
375:   }
376:   shell->tail =  ilink;
377:   shell->nmat += 1;
378:   return(0);
379: }

381: /*@
382:     MatCompositeAddMat - Add another matrix to a composite matrix.

384:    Collective on Mat

386:     Input Parameters:
387: +   mat - the composite matrix
388: -   smat - the partial matrix

390:    Level: advanced

392: .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE
393: @*/
394: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
395: {
396:   PetscErrorCode    ierr;

401:   PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));
402:   return(0);
403: }

405: static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type)
406: {
407:   Mat_Composite  *b = (Mat_Composite*)mat->data;

410:   b->type = type;
411:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
412:     mat->ops->getdiagonal   = 0;
413:     mat->ops->mult          = MatMult_Composite_Multiplicative;
414:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
415:   } else {
416:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
417:     mat->ops->mult          = MatMult_Composite;
418:     mat->ops->multtranspose = MatMultTranspose_Composite;
419:   }
420:   return(0);
421: }

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

426:    Logically Collective on Mat

428:    Input Parameters:
429: .  mat - the composite matrix

431:    Level: advanced

433: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE

435: @*/
436: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
437: {

443:   PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));
444:   return(0);
445: }

447: static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type)
448: {
449:   Mat_Composite  *b = (Mat_Composite*)mat->data;

452:   *type = b->type;
453:   return(0);
454: }

456: /*@
457:    MatCompositeGetType - Returns type of composite.

459:    Not Collective

461:    Input Parameter:
462: .  mat - the composite matrix

464:    Output Parameter:
465: .  type - type of composite

467:    Level: advanced

469: .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE

471: @*/
472: PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type)
473: {

479:   PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));
480:   return(0);
481: }

483: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str)
484: {
485:   Mat_Composite  *b = (Mat_Composite*)mat->data;

488:   b->structure = str;
489:   return(0);
490: }

492: /*@
493:    MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.

495:    Not Collective

497:    Input Parameters:
498: +  mat - the composite matrix
499: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN

501:    Level: advanced

503:    Notes:
504:     Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix.

506: .seealso: MatAXPY(), MatCreateComposite(), MatCompositeMerge() MatCompositeGetMatStructure(), MATCOMPOSITE

508: @*/
509: PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str)
510: {

515:   PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));
516:   return(0);
517: }

519: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str)
520: {
521:   Mat_Composite  *b = (Mat_Composite*)mat->data;

524:   *str = b->structure;
525:   return(0);
526: }

528: /*@
529:    MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.

531:    Not Collective

533:    Input Parameter:
534: .  mat - the composite matrix

536:    Output Parameter:
537: .  str - structure of the matrices

539:    Level: advanced

541: .seealso: MatCreateComposite(), MatCompositeSetMatStructure(), MATCOMPOSITE

543: @*/
544: PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str)
545: {

551:   PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));
552:   return(0);
553: }

555: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type)
556: {
557:   Mat_Composite     *shell = (Mat_Composite*)mat->data;

560:   shell->mergetype = type;
561:   return(0);
562: }

564: /*@
565:    MatCompositeSetMergeType - Sets order of MatCompositeMerge().

567:    Logically Collective on Mat

569:    Input Parameter:
570: +  mat - the composite matrix
571: -  type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]),
572:           MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1])

574:    Level: advanced

576:    Notes:
577:     The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed.
578:     If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
579:     otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].

581: .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE

583: @*/
584: PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type)
585: {

591:   PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));
592:   return(0);
593: }

595: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
596: {
597:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
598:   Mat_CompositeLink next   = shell->head, prev = shell->tail;
599:   PetscErrorCode    ierr;
600:   Mat               tmat,newmat;
601:   Vec               left,right;
602:   PetscScalar       scale;

605:   if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");

608:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
609:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
610:       MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
611:       while ((next = next->next)) {
612:         MatAXPY(tmat,1.0,next->mat,shell->structure);
613:       }
614:     } else {
615:       MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
616:       while ((prev = prev->prev)) {
617:         MatAXPY(tmat,1.0,prev->mat,shell->structure);
618:       }
619:     }
620:   } else {
621:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
622:       MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
623:       while ((next = next->next)) {
624:         MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
625:         MatDestroy(&tmat);
626:         tmat = newmat;
627:       }
628:     } else {
629:       MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
630:       while ((prev = prev->prev)) {
631:         MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
632:         MatDestroy(&tmat);
633:         tmat = newmat;
634:       }
635:     }
636:   }

638:   scale = shell->scale;
639:   if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
640:   if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}

642:   MatHeaderReplace(mat,&tmat);

644:   MatDiagonalScale(mat,left,right);
645:   MatScale(mat,scale);
646:   VecDestroy(&left);
647:   VecDestroy(&right);
648:   return(0);
649: }

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

655:   Collective

657:    Input Parameters:
658: .  mat - the composite matrix


661:    Options Database Keys:
662: +  -mat_composite_merge - merge in MatAssemblyEnd()
663: -  -mat_composite_merge_type - set merge direction

665:    Level: advanced

667:    Notes:
668:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
669:     matrix in the composite matrix.

671: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMatStructure(), MatCompositeSetMergeType(), MATCOMPOSITE

673: @*/
674: PetscErrorCode MatCompositeMerge(Mat mat)
675: {

680:   PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));
681:   return(0);
682: }

684: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat)
685: {
686:   Mat_Composite     *shell = (Mat_Composite*)mat->data;

689:   *nmat = shell->nmat;
690:   return(0);
691: }

693: /*@
694:    MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.

696:    Not Collective

698:    Input Parameter:
699: .  mat - the composite matrix

701:    Output Parameter:
702: .  nmat - number of matrices in the composite matrix

704:    Level: advanced

706: .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE

708: @*/
709: PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat)
710: {

716:   PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));
717:   return(0);
718: }

720: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai)
721: {
722:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
723:   Mat_CompositeLink ilink;
724:   PetscInt          k;

727:   if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat);
728:   ilink = shell->head;
729:   for (k=0; k<i; k++) {
730:     ilink = ilink->next;
731:   }
732:   *Ai = ilink->mat;
733:   return(0);
734: }

736: /*@
737:    MatCompositeGetMat - Returns the ith matrix from the composite matrix.

739:    Logically Collective on Mat

741:    Input Parameter:
742: +  mat - the composite matrix
743: -  i - the number of requested matrix

745:    Output Parameter:
746: .  Ai - ith matrix in composite

748:    Level: advanced

750: .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE

752: @*/
753: PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai)
754: {

761:   PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));
762:   return(0);
763: }

765: static struct _MatOps MatOps_Values = {0,
766:                                        0,
767:                                        0,
768:                                        MatMult_Composite,
769:                                        MatMultAdd_Composite,
770:                                 /*  5*/ MatMultTranspose_Composite,
771:                                        MatMultTransposeAdd_Composite,
772:                                        0,
773:                                        0,
774:                                        0,
775:                                 /* 10*/ 0,
776:                                        0,
777:                                        0,
778:                                        0,
779:                                        0,
780:                                 /* 15*/ 0,
781:                                        0,
782:                                        MatGetDiagonal_Composite,
783:                                        MatDiagonalScale_Composite,
784:                                        0,
785:                                 /* 20*/ 0,
786:                                        MatAssemblyEnd_Composite,
787:                                        0,
788:                                        0,
789:                                /* 24*/ 0,
790:                                        0,
791:                                        0,
792:                                        0,
793:                                        0,
794:                                /* 29*/ 0,
795:                                        0,
796:                                        0,
797:                                        0,
798:                                        0,
799:                                /* 34*/ 0,
800:                                        0,
801:                                        0,
802:                                        0,
803:                                        0,
804:                                /* 39*/ 0,
805:                                        0,
806:                                        0,
807:                                        0,
808:                                        0,
809:                                /* 44*/ 0,
810:                                        MatScale_Composite,
811:                                        MatShift_Basic,
812:                                        0,
813:                                        0,
814:                                /* 49*/ 0,
815:                                        0,
816:                                        0,
817:                                        0,
818:                                        0,
819:                                /* 54*/ 0,
820:                                        0,
821:                                        0,
822:                                        0,
823:                                        0,
824:                                /* 59*/ 0,
825:                                        MatDestroy_Composite,
826:                                        0,
827:                                        0,
828:                                        0,
829:                                /* 64*/ 0,
830:                                        0,
831:                                        0,
832:                                        0,
833:                                        0,
834:                                /* 69*/ 0,
835:                                        0,
836:                                        0,
837:                                        0,
838:                                        0,
839:                                /* 74*/ 0,
840:                                        0,
841:                                        MatSetFromOptions_Composite,
842:                                        0,
843:                                        0,
844:                                /* 79*/ 0,
845:                                        0,
846:                                        0,
847:                                        0,
848:                                        0,
849:                                /* 84*/ 0,
850:                                        0,
851:                                        0,
852:                                        0,
853:                                        0,
854:                                /* 89*/ 0,
855:                                        0,
856:                                        0,
857:                                        0,
858:                                        0,
859:                                /* 94*/ 0,
860:                                        0,
861:                                        0,
862:                                        0,
863:                                        0,
864:                                 /*99*/ 0,
865:                                        0,
866:                                        0,
867:                                        0,
868:                                        0,
869:                                /*104*/ 0,
870:                                        0,
871:                                        0,
872:                                        0,
873:                                        0,
874:                                /*109*/ 0,
875:                                        0,
876:                                        0,
877:                                        0,
878:                                        0,
879:                                /*114*/ 0,
880:                                        0,
881:                                        0,
882:                                        0,
883:                                        0,
884:                                /*119*/ 0,
885:                                        0,
886:                                        0,
887:                                        0,
888:                                        0,
889:                                /*124*/ 0,
890:                                        0,
891:                                        0,
892:                                        0,
893:                                        0,
894:                                /*129*/ 0,
895:                                        0,
896:                                        0,
897:                                        0,
898:                                        0,
899:                                /*134*/ 0,
900:                                        0,
901:                                        0,
902:                                        0,
903:                                        0,
904:                                /*139*/ 0,
905:                                        0,
906:                                        0
907: };

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

913:    Notes:
914:     to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);

916:   Level: advanced

918: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeSetMatStructure(), MatCompositeGetMatStructure(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat()
919: M*/

921: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
922: {
923:   Mat_Composite  *b;

927:   PetscNewLog(A,&b);
928:   A->data = (void*)b;
929:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

931:   PetscLayoutSetUp(A->rmap);
932:   PetscLayoutSetUp(A->cmap);

934:   A->assembled    = PETSC_TRUE;
935:   A->preallocated = PETSC_TRUE;
936:   b->type         = MAT_COMPOSITE_ADDITIVE;
937:   b->scale        = 1.0;
938:   b->nmat         = 0;
939:   b->merge        = PETSC_FALSE;
940:   b->mergetype    = MAT_COMPOSITE_MERGE_RIGHT;
941:   b->structure    = DIFFERENT_NONZERO_PATTERN;

943:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);
944:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);
945:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);
946:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);
947:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);
948:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);
949:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);
950:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);
951:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);

953:   PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
954:   return(0);
955: }