Actual source code: mcomposite.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  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:                                        0
398: };

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

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

405:   Level: advanced

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

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

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

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

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

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

438:   Collective on MPI_Comm

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

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

448:    Level: advanced

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

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

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

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

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

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

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

493:    Collective on Mat

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

499:    Level: advanced

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

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

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

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

536:   Collective on MPI_Comm

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


542:    Level: advanced

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

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

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

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


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

581:   Collective on MPI_Comm

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


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

590:    Level: advanced

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

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

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

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

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

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

630:   MatHeaderReplace(mat,tmat);

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