Actual source code: mcomposite.c

petsc-3.10.5 2019-03-28
Report Typos and Errors

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

  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;

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

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

 45: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
 46: {
 47:   Mat_Composite     *shell = (Mat_Composite*)A->data;
 48:   Mat_CompositeLink next   = shell->head;
 49:   PetscErrorCode    ierr;
 50:   Vec               in,out;

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

 79: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
 80: {
 81:   Mat_Composite     *shell = (Mat_Composite*)A->data;
 82:   Mat_CompositeLink tail   = shell->tail;
 83:   PetscErrorCode    ierr;
 84:   Vec               in,out;

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

113: PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
114: {
115:   Mat_Composite     *shell = (Mat_Composite*)A->data;
116:   Mat_CompositeLink next   = shell->head;
117:   PetscErrorCode    ierr;
118:   Vec               in;

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

141: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
142: {
143:   Mat_Composite     *shell = (Mat_Composite*)A->data;
144:   Mat_CompositeLink next   = shell->head;
145:   PetscErrorCode    ierr;
146:   Vec               in;

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

169: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
170: {
171:   Mat_Composite     *shell = (Mat_Composite*)A->data;
172:   Mat_CompositeLink next   = shell->head;
173:   PetscErrorCode    ierr;

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

179:   MatGetDiagonal(next->mat,v);
180:   if (next->next && !shell->work) {
181:     VecDuplicate(v,&shell->work);
182:   }
183:   while ((next = next->next)) {
184:     MatGetDiagonal(next->mat,shell->work);
185:     VecAXPY(v,1.0,shell->work);
186:   }
187:   VecScale(v,shell->scale);
188:   return(0);
189: }

191: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
192: {
194:   PetscBool      flg = PETSC_FALSE;

197:   PetscOptionsGetBool(((PetscObject)Y)->options,((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);
198:   if (flg) {
199:     MatCompositeMerge(Y);
200:   }
201:   return(0);
202: }

204: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
205: {
206:   Mat_Composite *a = (Mat_Composite*)inA->data;

209:   a->scale *= alpha;
210:   return(0);
211: }

213: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
214: {
215:   Mat_Composite  *a = (Mat_Composite*)inA->data;

219:   if (left) {
220:     if (!a->left) {
221:       VecDuplicate(left,&a->left);
222:       VecCopy(left,a->left);
223:     } else {
224:       VecPointwiseMult(a->left,left,a->left);
225:     }
226:   }
227:   if (right) {
228:     if (!a->right) {
229:       VecDuplicate(right,&a->right);
230:       VecCopy(right,a->right);
231:     } else {
232:       VecPointwiseMult(a->right,right,a->right);
233:     }
234:   }
235:   return(0);
236: }

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

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

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

388:   Level: advanced

390: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
391: M*/

393: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
394: {
395:   Mat_Composite  *b;

399:   PetscNewLog(A,&b);
400:   A->data = (void*)b;
401:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

403:   PetscLayoutSetUp(A->rmap);
404:   PetscLayoutSetUp(A->cmap);

406:   A->assembled    = PETSC_TRUE;
407:   A->preallocated = PETSC_TRUE;
408:   b->type         = MAT_COMPOSITE_ADDITIVE;
409:   b->scale        = 1.0;
410:   PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
411:   return(0);
412: }

414: /*@
415:    MatCreateComposite - Creates a matrix as the sum of zero or more matrices

417:   Collective on MPI_Comm

419:    Input Parameters:
420: +  comm - MPI communicator
421: .  nmat - number of matrices to put in
422: -  mats - the matrices

424:    Output Parameter:
425: .  mat - the matrix

427:    Level: advanced

429:    Notes:
430:      Alternative construction
431: $       MatCreate(comm,&mat);
432: $       MatSetSizes(mat,m,n,M,N);
433: $       MatSetType(mat,MATCOMPOSITE);
434: $       MatCompositeAddMat(mat,mats[0]);
435: $       ....
436: $       MatCompositeAddMat(mat,mats[nmat-1]);
437: $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
438: $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

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

444: @*/
445: PetscErrorCode  MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
446: {
448:   PetscInt       m,n,M,N,i;

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

454:   MatGetLocalSize(mats[0],&m,&n);
455:   MatGetSize(mats[0],&M,&N);
456:   MatCreate(comm,mat);
457:   MatSetSizes(*mat,m,n,M,N);
458:   MatSetType(*mat,MATCOMPOSITE);
459:   for (i=0; i<nmat; i++) {
460:     MatCompositeAddMat(*mat,mats[i]);
461:   }
462:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
463:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
464:   return(0);
465: }

467: /*@
468:     MatCompositeAddMat - add another matrix to a composite matrix

470:    Collective on Mat

472:     Input Parameters:
473: +   mat - the composite matrix
474: -   smat - the partial matrix

476:    Level: advanced

478: .seealso: MatCreateComposite()
479: @*/
480: PetscErrorCode  MatCompositeAddMat(Mat mat,Mat smat)
481: {
482:   Mat_Composite     *shell;
483:   PetscErrorCode    ierr;
484:   Mat_CompositeLink ilink,next;

489:   PetscNewLog(mat,&ilink);
490:   ilink->next = 0;
491:   PetscObjectReference((PetscObject)smat);
492:   ilink->mat  = smat;

494:   shell = (Mat_Composite*)mat->data;
495:   next  = shell->head;
496:   if (!next) shell->head = ilink;
497:   else {
498:     while (next->next) {
499:       next = next->next;
500:     }
501:     next->next  = ilink;
502:     ilink->prev = next;
503:   }
504:   shell->tail = ilink;
505:   return(0);
506: }

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

511:   Collective on MPI_Comm

513:    Input Parameters:
514: .  mat - the composite matrix


517:    Level: advanced

519:    Notes:
520:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
521:     matrix in the composite matrix.

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

525: @*/
526: PetscErrorCode  MatCompositeSetType(Mat mat,MatCompositeType type)
527: {
528:   Mat_Composite  *b = (Mat_Composite*)mat->data;
529:   PetscBool      flg;

533:   PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
534:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
535:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
536:     mat->ops->getdiagonal   = 0;
537:     mat->ops->mult          = MatMult_Composite_Multiplicative;
538:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
539:     b->type                 = MAT_COMPOSITE_MULTIPLICATIVE;
540:   } else {
541:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
542:     mat->ops->mult          = MatMult_Composite;
543:     mat->ops->multtranspose = MatMultTranspose_Composite;
544:     b->type                 = MAT_COMPOSITE_ADDITIVE;
545:   }
546:   return(0);
547: }


550: /*@
551:    MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
552:      by summing all the matrices inside the composite matrix.

554:   Collective on MPI_Comm

556:    Input Parameters:
557: .  mat - the composite matrix


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

563:    Level: advanced

565:    Notes:
566:       The MatType of the resulting matrix will be the same as the MatType of the FIRST
567:     matrix in the composite matrix.

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

571: @*/
572: PetscErrorCode  MatCompositeMerge(Mat mat)
573: {
574:   Mat_Composite     *shell = (Mat_Composite*)mat->data;
575:   Mat_CompositeLink next   = shell->head, prev = shell->tail;
576:   PetscErrorCode    ierr;
577:   Mat               tmat,newmat;
578:   Vec               left,right;
579:   PetscScalar       scale;

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

585:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
586:     MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
587:     while ((next = next->next)) {
588:       MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
589:     }
590:   } else {
591:     MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
592:     while ((prev = prev->prev)) {
593:       MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
594:       MatDestroy(&tmat);
595:       tmat = newmat;
596:     }
597:   }

599:   scale = shell->scale;
600:   if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
601:   if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}

603:   MatHeaderReplace(mat,&tmat);

605:   MatDiagonalScale(mat,left,right);
606:   MatScale(mat,scale);
607:   VecDestroy(&left);
608:   VecDestroy(&right);
609:   return(0);
610: }