Actual source code: mcomposite.c

petsc-3.9.4 2018-09-11
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: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);

387:   Level: advanced

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

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

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

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

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

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

416:   Collective on MPI_Comm

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

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

426:    Level: advanced

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

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

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

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

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

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

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

469:    Collective on Mat

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

475:    Level: advanced

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

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

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

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

510:   Collective on MPI_Comm

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


516:    Level: advanced

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

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

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

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


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

553:   Collective on MPI_Comm

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


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

562:    Level: advanced

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

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

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

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

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

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

602:   MatHeaderReplace(mat,&tmat);

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