Actual source code: mcomposite.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/matimpl.h>        /*I "petscmat.h" I*/

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

 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: } Mat_Composite;

 22: PetscErrorCode MatDestroy_Composite(Mat mat)
 23: {
 24:   PetscErrorCode    ierr;
 25:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
 26:   Mat_CompositeLink next   = shell->head,oldnext;

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

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

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

 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:       MatGetVecs(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: }

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

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

151: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
152: {
153:   Mat_Composite     *shell = (Mat_Composite*)A->data;
154:   Mat_CompositeLink next   = shell->head;
155:   PetscErrorCode    ierr;
156:   Vec               in;

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

181: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
182: {
183:   Mat_Composite     *shell = (Mat_Composite*)A->data;
184:   Mat_CompositeLink next   = shell->head;
185:   PetscErrorCode    ierr;

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

191:   MatGetDiagonal(next->mat,v);
192:   if (next->next && !shell->work) {
193:     VecDuplicate(v,&shell->work);
194:   }
195:   while ((next = next->next)) {
196:     MatGetDiagonal(next->mat,shell->work);
197:     VecAXPY(v,1.0,shell->work);
198:   }
199:   VecScale(v,shell->scale);
200:   return(0);
201: }

205: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
206: {
208:   PetscBool      flg = PETSC_FALSE;

211:   PetscOptionsGetBool(((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);
212:   if (flg) {
213:     MatCompositeMerge(Y);
214:   }
215:   return(0);
216: }

220: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
221: {
222:   Mat_Composite *a = (Mat_Composite*)inA->data;

225:   a->scale *= alpha;
226:   return(0);
227: }

231: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
232: {
233:   Mat_Composite  *a = (Mat_Composite*)inA->data;

237:   if (left) {
238:     if (!a->left) {
239:       VecDuplicate(left,&a->left);
240:       VecCopy(left,a->left);
241:     } else {
242:       VecPointwiseMult(a->left,left,a->left);
243:     }
244:   }
245:   if (right) {
246:     if (!a->right) {
247:       VecDuplicate(right,&a->right);
248:       VecCopy(right,a->right);
249:     } else {
250:       VecPointwiseMult(a->right,right,a->right);
251:     }
252:   }
253:   return(0);
254: }

256: static struct _MatOps MatOps_Values = {0,
257:                                        0,
258:                                        0,
259:                                        MatMult_Composite,
260:                                        0,
261:                                 /*  5*/ MatMultTranspose_Composite,
262:                                        0,
263:                                        0,
264:                                        0,
265:                                        0,
266:                                 /* 10*/ 0,
267:                                        0,
268:                                        0,
269:                                        0,
270:                                        0,
271:                                 /* 15*/ 0,
272:                                        0,
273:                                        MatGetDiagonal_Composite,
274:                                        MatDiagonalScale_Composite,
275:                                        0,
276:                                 /* 20*/ 0,
277:                                        MatAssemblyEnd_Composite,
278:                                        0,
279:                                        0,
280:                                /* 24*/ 0,
281:                                        0,
282:                                        0,
283:                                        0,
284:                                        0,
285:                                /* 29*/ 0,
286:                                        0,
287:                                        0,
288:                                        0,
289:                                        0,
290:                                /* 34*/ 0,
291:                                        0,
292:                                        0,
293:                                        0,
294:                                        0,
295:                                /* 39*/ 0,
296:                                        0,
297:                                        0,
298:                                        0,
299:                                        0,
300:                                /* 44*/ 0,
301:                                        MatScale_Composite,
302:                                        0,
303:                                        0,
304:                                        0,
305:                                /* 49*/ 0,
306:                                        0,
307:                                        0,
308:                                        0,
309:                                        0,
310:                                /* 54*/ 0,
311:                                        0,
312:                                        0,
313:                                        0,
314:                                        0,
315:                                /* 59*/ 0,
316:                                        MatDestroy_Composite,
317:                                        0,
318:                                        0,
319:                                        0,
320:                                /* 64*/ 0,
321:                                        0,
322:                                        0,
323:                                        0,
324:                                        0,
325:                                /* 69*/ 0,
326:                                        0,
327:                                        0,
328:                                        0,
329:                                        0,
330:                                /* 74*/ 0,
331:                                        0,
332:                                        0,
333:                                        0,
334:                                        0,
335:                                /* 79*/ 0,
336:                                        0,
337:                                        0,
338:                                        0,
339:                                        0,
340:                                /* 84*/ 0,
341:                                        0,
342:                                        0,
343:                                        0,
344:                                        0,
345:                                /* 89*/ 0,
346:                                        0,
347:                                        0,
348:                                        0,
349:                                        0,
350:                                /* 94*/ 0,
351:                                        0,
352:                                        0,
353:                                        0,
354:                                        0,
355:                                 /*99*/ 0,
356:                                        0,
357:                                        0,
358:                                        0,
359:                                        0,
360:                                /*104*/ 0,
361:                                        0,
362:                                        0,
363:                                        0,
364:                                        0,
365:                                /*109*/ 0,
366:                                        0,
367:                                        0,
368:                                        0,
369:                                        0,
370:                                /*114*/ 0,
371:                                        0,
372:                                        0,
373:                                        0,
374:                                        0,
375:                                /*119*/ 0,
376:                                        0,
377:                                        0,
378:                                        0,
379:                                        0,
380:                                /*124*/ 0,
381:                                        0,
382:                                        0,
383:                                        0,
384:                                        0,
385:                                /*129*/ 0,
386:                                        0,
387:                                        0,
388:                                        0,
389:                                        0,
390:                                /*134*/ 0,
391:                                        0,
392:                                        0,
393:                                        0,
394:                                        0,
395:                                /*139*/ 0,
396:                                        0
397: };

399: /*MC
400:    MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout).

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

404:   Level: advanced

406: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
407: M*/

411: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
412: {
413:   Mat_Composite  *b;

417:   PetscNewLog(A,Mat_Composite,&b);
418:   A->data = (void*)b;
419:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

421:   PetscLayoutSetUp(A->rmap);
422:   PetscLayoutSetUp(A->cmap);

424:   A->assembled    = PETSC_TRUE;
425:   A->preallocated = PETSC_TRUE;
426:   b->type         = MAT_COMPOSITE_ADDITIVE;
427:   b->scale        = 1.0;
428:   PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
429:   return(0);
430: }

434: /*@C
435:    MatCreateComposite - Creates a matrix as the sum of zero or more matrices

437:   Collective on MPI_Comm

439:    Input Parameters:
440: +  comm - MPI communicator
441: .  nmat - number of matrices to put in
442: -  mats - the matrices

444:    Output Parameter:
445: .  A - the matrix

447:    Level: advanced

449:    Notes:
450:      Alternative construction
451: $       MatCreate(comm,&mat);
452: $       MatSetSizes(mat,m,n,M,N);
453: $       MatSetType(mat,MATCOMPOSITE);
454: $       MatCompositeAddMat(mat,mats[0]);
455: $       ....
456: $       MatCompositeAddMat(mat,mats[nmat-1]);
457: $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
458: $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

462: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType

464: @*/
465: PetscErrorCode  MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
466: {
468:   PetscInt       m,n,M,N,i;

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

474:   MatGetLocalSize(mats[0],&m,&n);
475:   MatGetSize(mats[0],&M,&N);
476:   MatCreate(comm,mat);
477:   MatSetSizes(*mat,m,n,M,N);
478:   MatSetType(*mat,MATCOMPOSITE);
479:   for (i=0; i<nmat; i++) {
480:     MatCompositeAddMat(*mat,mats[i]);
481:   }
482:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
483:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
484:   return(0);
485: }

489: /*@
490:     MatCompositeAddMat - add another matrix to a composite matrix

492:    Collective on Mat

494:     Input Parameters:
495: +   mat - the composite matrix
496: -   smat - the partial matrix

498:    Level: advanced

500: .seealso: MatCreateComposite()
501: @*/
502: PetscErrorCode  MatCompositeAddMat(Mat mat,Mat smat)
503: {
504:   Mat_Composite     *shell;
505:   PetscErrorCode    ierr;
506:   Mat_CompositeLink ilink,next;

511:   PetscNewLog(mat,struct _Mat_CompositeLink,&ilink);
512:   ilink->next = 0;
513:   PetscObjectReference((PetscObject)smat);
514:   ilink->mat  = smat;

516:   shell = (Mat_Composite*)mat->data;
517:   next  = shell->head;
518:   if (!next) shell->head = ilink;
519:   else {
520:     while (next->next) {
521:       next = next->next;
522:     }
523:     next->next  = ilink;
524:     ilink->prev = next;
525:   }
526:   shell->tail = ilink;
527:   return(0);
528: }

532: /*@C
533:    MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product

535:   Collective on MPI_Comm

537:    Input Parameters:
538: .  mat - the composite matrix


541:    Level: advanced

543:    Notes:
544:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
545:     matrix in the composite matrix.

547: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE

549: @*/
550: PetscErrorCode  MatCompositeSetType(Mat mat,MatCompositeType type)
551: {
552:   Mat_Composite  *b = (Mat_Composite*)mat->data;
553:   PetscBool      flg;

557:   PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
558:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
559:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
560:     mat->ops->getdiagonal   = 0;
561:     mat->ops->mult          = MatMult_Composite_Multiplicative;
562:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
563:     b->type                 = MAT_COMPOSITE_MULTIPLICATIVE;
564:   } else {
565:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
566:     mat->ops->mult          = MatMult_Composite;
567:     mat->ops->multtranspose = MatMultTranspose_Composite;
568:     b->type                 = MAT_COMPOSITE_ADDITIVE;
569:   }
570:   return(0);
571: }


576: /*@C
577:    MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
578:      by summing all the matrices inside the composite matrix.

580:   Collective on MPI_Comm

582:    Input Parameters:
583: .  mat - the composite matrix


586:    Options Database:
587: .  -mat_composite_merge  (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)

589:    Level: advanced

591:    Notes:
592:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
593:     matrix in the composite matrix.

595: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE

597: @*/
598: PetscErrorCode  MatCompositeMerge(Mat mat)
599: {
600:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
601:   Mat_CompositeLink next   = shell->head, prev = shell->tail;
602:   PetscErrorCode    ierr;
603:   Mat               tmat,newmat;
604:   Vec               left,right;
605:   PetscScalar       scale;

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

611:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
612:     MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
613:     while ((next = next->next)) {
614:       MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
615:     }
616:   } else {
617:     MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
618:     while ((prev = prev->prev)) {
619:       MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
620:       MatDestroy(&tmat);
621:       tmat = newmat;
622:     }
623:   }

625:   scale = shell->scale;
626:   if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
627:   if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}

629:   MatHeaderReplace(mat,tmat);

631:   MatDiagonalScale(mat,left,right);
632:   MatScale(mat,scale);
633:   VecDestroy(&left);
634:   VecDestroy(&right);
635:   return(0);
636: }