Actual source code: mcomposite.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

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

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

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

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

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

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

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

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

 68:   PetscFree(shell->scalings);
 69:   PetscFree(mat->data);
 70:   return(0);
 71: }

 73: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
 74: {
 75:   Mat_Composite     *shell = (Mat_Composite*)A->data;
 76:   Mat_CompositeLink next   = shell->head;
 77:   PetscErrorCode    ierr;
 78:   Vec               in,out;
 79:   PetscScalar       scale;
 80:   PetscInt          i;

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

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

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

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

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

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

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

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

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

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

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

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

264:   if (shell->Mvctx) { /* Have a merged Mvctx */
265:     /* Suppose we want to compute y = sMx, where s is the scaling factor and A, B are matrix M's diagonal/off-diagonal part. We could do
266:        in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an oppertunity to
267:        overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
268:      */
269:     VecScatterBegin(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);
270:     VecScatterEnd(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);

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

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

297:   if (shell->left) {VecPointwiseMult(y,shell->left,y);}
298:   VecScale(y,shell->scale);
299:   return(0);
300: }

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

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

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

342: PetscErrorCode MatMultAdd_Composite(Mat A,Vec x,Vec y,Vec z)
343: {
344:   Mat_Composite     *shell = (Mat_Composite*)A->data;
345:   PetscErrorCode    ierr;

348:   if (y != z) {
349:     MatMult(A,x,z);
350:     VecAXPY(z,1.0,y);
351:   } else {
352:     if (!shell->leftwork) {
353:       VecDuplicate(z,&shell->leftwork);
354:     }
355:     MatMult(A,x,shell->leftwork);
356:     VecCopy(y,z);
357:     VecAXPY(z,1.0,shell->leftwork);
358:   }
359:   return(0);
360: }

362: PetscErrorCode MatMultTransposeAdd_Composite(Mat A,Vec x,Vec y, Vec z)
363: {
364:   Mat_Composite     *shell = (Mat_Composite*)A->data;
365:   PetscErrorCode    ierr;

368:   if (y != z) {
369:     MatMultTranspose(A,x,z);
370:     VecAXPY(z,1.0,y);
371:   } else {
372:     if (!shell->rightwork) {
373:       VecDuplicate(z,&shell->rightwork);
374:     }
375:     MatMultTranspose(A,x,shell->rightwork);
376:     VecCopy(y,z);
377:     VecAXPY(z,1.0,shell->rightwork);
378:   }
379:   return(0);
380: }

382: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
383: {
384:   Mat_Composite     *shell = (Mat_Composite*)A->data;
385:   Mat_CompositeLink next   = shell->head;
386:   PetscErrorCode    ierr;
387:   PetscInt          i;

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

393:   MatGetDiagonal(next->mat,v);
394:   if (shell->scalings) {VecScale(v,shell->scalings[0]);}

396:   if (next->next && !shell->work) {
397:     VecDuplicate(v,&shell->work);
398:   }
399:   i = 1;
400:   while ((next = next->next)) {
401:     MatGetDiagonal(next->mat,shell->work);
402:     VecAXPY(v,(shell->scalings ? shell->scalings[i++] : 1.0),shell->work);
403:   }
404:   VecScale(v,shell->scale);
405:   return(0);
406: }

408: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
409: {
410:   Mat_Composite     *shell = (Mat_Composite*)Y->data;
411:   PetscErrorCode    ierr;

414:   if (shell->merge) {
415:     MatCompositeMerge(Y);
416:   }

418:   return(0);
419: }

421: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
422: {
423:   Mat_Composite *a = (Mat_Composite*)inA->data;

426:   a->scale *= alpha;
427:   return(0);
428: }

430: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
431: {
432:   Mat_Composite  *a = (Mat_Composite*)inA->data;

436:   if (left) {
437:     if (!a->left) {
438:       VecDuplicate(left,&a->left);
439:       VecCopy(left,a->left);
440:     } else {
441:       VecPointwiseMult(a->left,left,a->left);
442:     }
443:   }
444:   if (right) {
445:     if (!a->right) {
446:       VecDuplicate(right,&a->right);
447:       VecCopy(right,a->right);
448:     } else {
449:       VecPointwiseMult(a->right,right,a->right);
450:     }
451:   }
452:   return(0);
453: }

455: PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A)
456: {
457:   Mat_Composite  *a = (Mat_Composite*)A->data;

461:   PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");
462:   PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);
463:   PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);
464:   PetscOptionsBool("-mat_composite_merge_mvctx","Merge MatMult() vecscat contexts","MatCreateComposite",a->merge_mvctx,&a->merge_mvctx,NULL);
465:   PetscOptionsTail();
466:   return(0);
467: }

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

472:   Collective

474:    Input Parameters:
475: +  comm - MPI communicator
476: .  nmat - number of matrices to put in
477: -  mats - the matrices

479:    Output Parameter:
480: .  mat - the matrix

482:    Options Database Keys:
483: +  -mat_composite_merge         - merge in MatAssemblyEnd()
484: .  -mat_composite_merge_mvctx   - merge Mvctx of component matrices to optimize communication in MatMult() for ADDITIVE matrices
485: -  -mat_composite_merge_type    - set merge direction

487:    Level: advanced

489:    Notes:
490:      Alternative construction
491: $       MatCreate(comm,&mat);
492: $       MatSetSizes(mat,m,n,M,N);
493: $       MatSetType(mat,MATCOMPOSITE);
494: $       MatCompositeAddMat(mat,mats[0]);
495: $       ....
496: $       MatCompositeAddMat(mat,mats[nmat-1]);
497: $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
498: $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

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

504: @*/
505: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
506: {
508:   PetscInt       m,n,M,N,i;

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

514:   MatGetLocalSize(mats[0],PETSC_IGNORE,&n);
515:   MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);
516:   MatGetSize(mats[0],PETSC_IGNORE,&N);
517:   MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);
518:   MatCreate(comm,mat);
519:   MatSetSizes(*mat,m,n,M,N);
520:   MatSetType(*mat,MATCOMPOSITE);
521:   for (i=0; i<nmat; i++) {
522:     MatCompositeAddMat(*mat,mats[i]);
523:   }
524:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
525:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
526:   return(0);
527: }


530: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat)
531: {
532:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
533:   Mat_CompositeLink ilink,next = shell->head;
534:   PetscErrorCode    ierr;

537:   PetscNewLog(mat,&ilink);
538:   ilink->next = NULL;
539:   PetscObjectReference((PetscObject)smat);
540:   ilink->mat  = smat;

542:   if (!next) shell->head = ilink;
543:   else {
544:     while (next->next) {
545:       next = next->next;
546:     }
547:     next->next  = ilink;
548:     ilink->prev = next;
549:   }
550:   shell->tail =  ilink;
551:   shell->nmat += 1;

553:   /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
554:   if (shell->scalings) {
555:     PetscRealloc(sizeof(PetscScalar)*shell->nmat,&shell->scalings);
556:     shell->scalings[shell->nmat-1] = 1.0;
557:   }
558:   return(0);
559: }

561: /*@
562:     MatCompositeAddMat - Add another matrix to a composite matrix.

564:    Collective on Mat

566:     Input Parameters:
567: +   mat - the composite matrix
568: -   smat - the partial matrix

570:    Level: advanced

572: .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE
573: @*/
574: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
575: {
576:   PetscErrorCode    ierr;

581:   PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));
582:   return(0);
583: }

585: static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type)
586: {
587:   Mat_Composite  *b = (Mat_Composite*)mat->data;

590:   b->type = type;
591:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
592:     mat->ops->getdiagonal   = NULL;
593:     mat->ops->mult          = MatMult_Composite_Multiplicative;
594:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
595:     b->merge_mvctx          = PETSC_FALSE;
596:   } else {
597:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
598:     mat->ops->mult          = MatMult_Composite;
599:     mat->ops->multtranspose = MatMultTranspose_Composite;
600:   }
601:   return(0);
602: }

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

607:    Logically Collective on Mat

609:    Input Parameters:
610: .  mat - the composite matrix

612:    Level: advanced

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

616: @*/
617: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
618: {

624:   PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));
625:   return(0);
626: }

628: static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type)
629: {
630:   Mat_Composite  *b = (Mat_Composite*)mat->data;

633:   *type = b->type;
634:   return(0);
635: }

637: /*@
638:    MatCompositeGetType - Returns type of composite.

640:    Not Collective

642:    Input Parameter:
643: .  mat - the composite matrix

645:    Output Parameter:
646: .  type - type of composite

648:    Level: advanced

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

652: @*/
653: PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type)
654: {

660:   PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));
661:   return(0);
662: }

664: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str)
665: {
666:   Mat_Composite  *b = (Mat_Composite*)mat->data;

669:   b->structure = str;
670:   return(0);
671: }

673: /*@
674:    MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.

676:    Not Collective

678:    Input Parameters:
679: +  mat - the composite matrix
680: -  str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN

682:    Level: advanced

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

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

689: @*/
690: PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str)
691: {

696:   PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));
697:   return(0);
698: }

700: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str)
701: {
702:   Mat_Composite  *b = (Mat_Composite*)mat->data;

705:   *str = b->structure;
706:   return(0);
707: }

709: /*@
710:    MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.

712:    Not Collective

714:    Input Parameter:
715: .  mat - the composite matrix

717:    Output Parameter:
718: .  str - structure of the matrices

720:    Level: advanced

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

724: @*/
725: PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str)
726: {

732:   PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));
733:   return(0);
734: }

736: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type)
737: {
738:   Mat_Composite     *shell = (Mat_Composite*)mat->data;

741:   shell->mergetype = type;
742:   return(0);
743: }

745: /*@
746:    MatCompositeSetMergeType - Sets order of MatCompositeMerge().

748:    Logically Collective on Mat

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

755:    Level: advanced

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

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

764: @*/
765: PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type)
766: {

772:   PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));
773:   return(0);
774: }

776: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
777: {
778:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
779:   Mat_CompositeLink next   = shell->head, prev = shell->tail;
780:   PetscErrorCode    ierr;
781:   Mat               tmat,newmat;
782:   Vec               left,right;
783:   PetscScalar       scale;
784:   PetscInt          i;

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

790:   scale = shell->scale;
791:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
792:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
793:       i = 0;
794:       MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
795:       if (shell->scalings) {MatScale(tmat,shell->scalings[i++]);}
796:       while ((next = next->next)) {
797:         MatAXPY(tmat,(shell->scalings ? shell->scalings[i++] : 1.0),next->mat,shell->structure);
798:       }
799:     } else {
800:       i = shell->nmat-1;
801:       MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
802:       if (shell->scalings) {MatScale(tmat,shell->scalings[i--]);}
803:       while ((prev = prev->prev)) {
804:         MatAXPY(tmat,(shell->scalings ? shell->scalings[i--] : 1.0),prev->mat,shell->structure);
805:       }
806:     }
807:   } else {
808:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
809:       MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
810:       while ((next = next->next)) {
811:         MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
812:         MatDestroy(&tmat);
813:         tmat = newmat;
814:       }
815:     } else {
816:       MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
817:       while ((prev = prev->prev)) {
818:         MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
819:         MatDestroy(&tmat);
820:         tmat = newmat;
821:       }
822:     }
823:     if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];}
824:   }

826:   if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
827:   if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}

829:   MatHeaderReplace(mat,&tmat);

831:   MatDiagonalScale(mat,left,right);
832:   MatScale(mat,scale);
833:   VecDestroy(&left);
834:   VecDestroy(&right);
835:   return(0);
836: }

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

842:   Collective

844:    Input Parameters:
845: .  mat - the composite matrix


848:    Options Database Keys:
849: +  -mat_composite_merge - merge in MatAssemblyEnd()
850: -  -mat_composite_merge_type - set merge direction

852:    Level: advanced

854:    Notes:
855:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
856:     matrix in the composite matrix.

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

860: @*/
861: PetscErrorCode MatCompositeMerge(Mat mat)
862: {

867:   PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));
868:   return(0);
869: }

871: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat)
872: {
873:   Mat_Composite     *shell = (Mat_Composite*)mat->data;

876:   *nmat = shell->nmat;
877:   return(0);
878: }

880: /*@
881:    MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.

883:    Not Collective

885:    Input Parameter:
886: .  mat - the composite matrix

888:    Output Parameter:
889: .  nmat - number of matrices in the composite matrix

891:    Level: advanced

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

895: @*/
896: PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat)
897: {

903:   PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));
904:   return(0);
905: }

907: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai)
908: {
909:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
910:   Mat_CompositeLink ilink;
911:   PetscInt          k;

914:   if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat);
915:   ilink = shell->head;
916:   for (k=0; k<i; k++) {
917:     ilink = ilink->next;
918:   }
919:   *Ai = ilink->mat;
920:   return(0);
921: }

923: /*@
924:    MatCompositeGetMat - Returns the ith matrix from the composite matrix.

926:    Logically Collective on Mat

928:    Input Parameter:
929: +  mat - the composite matrix
930: -  i - the number of requested matrix

932:    Output Parameter:
933: .  Ai - ith matrix in composite

935:    Level: advanced

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

939: @*/
940: PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai)
941: {

948:   PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));
949:   return(0);
950: }

952: PetscErrorCode MatCompositeSetScalings_Composite(Mat mat,const PetscScalar *scalings)
953: {
955:   Mat_Composite  *shell = (Mat_Composite*)mat->data;
956:   PetscInt       nmat;

959:   MatCompositeGetNumberMat(mat,&nmat);
960:   if (!shell->scalings) {PetscMalloc1(nmat,&shell->scalings);}
961:   PetscArraycpy(shell->scalings,scalings,nmat);
962:   return(0);
963: }

965: /*@
966:    MatCompositeSetScalings - Sets separate scaling factors for component matrices.

968:    Logically Collective on Mat

970:    Input Parameter:
971: +  mat      - the composite matrix
972: -  scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)

974:    Level: advanced

976: .seealso: MatScale(), MatDiagonalScale(), MATCOMPOSITE

978: @*/
979: PetscErrorCode MatCompositeSetScalings(Mat mat,const PetscScalar *scalings)
980: {

987:   PetscUseMethod(mat,"MatCompositeSetScalings_C",(Mat,const PetscScalar*),(mat,scalings));
988:   return(0);
989: }

991: static struct _MatOps MatOps_Values = {NULL,
992:                                        NULL,
993:                                        NULL,
994:                                        MatMult_Composite,
995:                                        MatMultAdd_Composite,
996:                                /*  5*/ MatMultTranspose_Composite,
997:                                        MatMultTransposeAdd_Composite,
998:                                        NULL,
999:                                        NULL,
1000:                                        NULL,
1001:                                /* 10*/ NULL,
1002:                                        NULL,
1003:                                        NULL,
1004:                                        NULL,
1005:                                        NULL,
1006:                                /* 15*/ NULL,
1007:                                        NULL,
1008:                                        MatGetDiagonal_Composite,
1009:                                        MatDiagonalScale_Composite,
1010:                                        NULL,
1011:                                /* 20*/ NULL,
1012:                                        MatAssemblyEnd_Composite,
1013:                                        NULL,
1014:                                        NULL,
1015:                                /* 24*/ NULL,
1016:                                        NULL,
1017:                                        NULL,
1018:                                        NULL,
1019:                                        NULL,
1020:                                /* 29*/ NULL,
1021:                                        NULL,
1022:                                        NULL,
1023:                                        NULL,
1024:                                        NULL,
1025:                                /* 34*/ NULL,
1026:                                        NULL,
1027:                                        NULL,
1028:                                        NULL,
1029:                                        NULL,
1030:                                /* 39*/ NULL,
1031:                                        NULL,
1032:                                        NULL,
1033:                                        NULL,
1034:                                        NULL,
1035:                                /* 44*/ NULL,
1036:                                        MatScale_Composite,
1037:                                        MatShift_Basic,
1038:                                        NULL,
1039:                                        NULL,
1040:                                /* 49*/ NULL,
1041:                                        NULL,
1042:                                        NULL,
1043:                                        NULL,
1044:                                        NULL,
1045:                                /* 54*/ NULL,
1046:                                        NULL,
1047:                                        NULL,
1048:                                        NULL,
1049:                                        NULL,
1050:                                /* 59*/ NULL,
1051:                                        MatDestroy_Composite,
1052:                                        NULL,
1053:                                        NULL,
1054:                                        NULL,
1055:                                /* 64*/ NULL,
1056:                                        NULL,
1057:                                        NULL,
1058:                                        NULL,
1059:                                        NULL,
1060:                                /* 69*/ NULL,
1061:                                        NULL,
1062:                                        NULL,
1063:                                        NULL,
1064:                                        NULL,
1065:                                /* 74*/ NULL,
1066:                                        NULL,
1067:                                        MatSetFromOptions_Composite,
1068:                                        NULL,
1069:                                        NULL,
1070:                                /* 79*/ NULL,
1071:                                        NULL,
1072:                                        NULL,
1073:                                        NULL,
1074:                                        NULL,
1075:                                /* 84*/ NULL,
1076:                                        NULL,
1077:                                        NULL,
1078:                                        NULL,
1079:                                        NULL,
1080:                                /* 89*/ NULL,
1081:                                        NULL,
1082:                                        NULL,
1083:                                        NULL,
1084:                                        NULL,
1085:                                /* 94*/ NULL,
1086:                                        NULL,
1087:                                        NULL,
1088:                                        NULL,
1089:                                        NULL,
1090:                                 /*99*/ NULL,
1091:                                        NULL,
1092:                                        NULL,
1093:                                        NULL,
1094:                                        NULL,
1095:                                /*104*/ NULL,
1096:                                        NULL,
1097:                                        NULL,
1098:                                        NULL,
1099:                                        NULL,
1100:                                /*109*/ NULL,
1101:                                        NULL,
1102:                                        NULL,
1103:                                        NULL,
1104:                                        NULL,
1105:                                /*114*/ NULL,
1106:                                        NULL,
1107:                                        NULL,
1108:                                        NULL,
1109:                                        NULL,
1110:                                /*119*/ NULL,
1111:                                        NULL,
1112:                                        NULL,
1113:                                        NULL,
1114:                                        NULL,
1115:                                /*124*/ NULL,
1116:                                        NULL,
1117:                                        NULL,
1118:                                        NULL,
1119:                                        NULL,
1120:                                /*129*/ NULL,
1121:                                        NULL,
1122:                                        NULL,
1123:                                        NULL,
1124:                                        NULL,
1125:                                /*134*/ NULL,
1126:                                        NULL,
1127:                                        NULL,
1128:                                        NULL,
1129:                                        NULL,
1130:                                /*139*/ NULL,
1131:                                        NULL,
1132:                                        NULL
1133: };

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

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

1142:   Level: advanced

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

1147: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1148: {
1149:   Mat_Composite  *b;

1153:   PetscNewLog(A,&b);
1154:   A->data = (void*)b;
1155:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

1157:   PetscLayoutSetUp(A->rmap);
1158:   PetscLayoutSetUp(A->cmap);

1160:   A->assembled    = PETSC_TRUE;
1161:   A->preallocated = PETSC_TRUE;
1162:   b->type         = MAT_COMPOSITE_ADDITIVE;
1163:   b->scale        = 1.0;
1164:   b->nmat         = 0;
1165:   b->merge        = PETSC_FALSE;
1166:   b->mergetype    = MAT_COMPOSITE_MERGE_RIGHT;
1167:   b->structure    = DIFFERENT_NONZERO_PATTERN;
1168:   b->merge_mvctx  = PETSC_TRUE;



1172:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);
1173:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);
1174:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);
1175:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);
1176:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);
1177:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);
1178:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);
1179:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);
1180:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);
1181:   PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetScalings_C",MatCompositeSetScalings_Composite);

1183:   PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
1184:   return(0);
1185: }