Actual source code: mcomposite.c


  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:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
 38:   Mat_CompositeLink next  = shell->head, oldnext;
 39:   PetscInt          i;

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

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

 64:   PetscFree(shell->scalings);
 65:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL);
 66:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL);
 67:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL);
 68:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL);
 69:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL);
 70:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL);
 71:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL);
 72:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL);
 73:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL);
 74:   PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL);
 75:   PetscFree(mat->data);
 76:   return 0;
 77: }

 79: PetscErrorCode MatMult_Composite_Multiplicative(Mat A, Vec x, Vec y)
 80: {
 81:   Mat_Composite    *shell = (Mat_Composite *)A->data;
 82:   Mat_CompositeLink next  = shell->head;
 83:   Vec               in, out;
 84:   PetscScalar       scale;
 85:   PetscInt          i;

 88:   in = x;
 89:   if (shell->right) {
 90:     if (!shell->rightwork) VecDuplicate(shell->right, &shell->rightwork);
 91:     VecPointwiseMult(shell->rightwork, shell->right, in);
 92:     in = shell->rightwork;
 93:   }
 94:   while (next->next) {
 95:     if (!next->work) { /* should reuse previous work if the same size */
 96:       MatCreateVecs(next->mat, NULL, &next->work);
 97:     }
 98:     out = next->work;
 99:     MatMult(next->mat, in, out);
100:     in   = out;
101:     next = next->next;
102:   }
103:   MatMult(next->mat, in, y);
104:   if (shell->left) VecPointwiseMult(y, shell->left, y);
105:   scale = shell->scale;
106:   if (shell->scalings) {
107:     for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
108:   }
109:   VecScale(y, scale);
110:   return 0;
111: }

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

122:   in = x;
123:   if (shell->left) {
124:     if (!shell->leftwork) VecDuplicate(shell->left, &shell->leftwork);
125:     VecPointwiseMult(shell->leftwork, shell->left, in);
126:     in = shell->leftwork;
127:   }
128:   while (tail->prev) {
129:     if (!tail->prev->work) { /* should reuse previous work if the same size */
130:       MatCreateVecs(tail->mat, NULL, &tail->prev->work);
131:     }
132:     out = tail->prev->work;
133:     MatMultTranspose(tail->mat, in, out);
134:     in   = out;
135:     tail = tail->prev;
136:   }
137:   MatMultTranspose(tail->mat, in, y);
138:   if (shell->right) VecPointwiseMult(y, shell->right, y);

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

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

161:   in = x;
162:   if (shell->right) {
163:     if (!shell->rightwork) VecDuplicate(shell->right, &shell->rightwork);
164:     VecPointwiseMult(shell->rightwork, shell->right, in);
165:     in = shell->rightwork;
166:   }

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

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

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

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

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

256: skip_merge_mvctx:
257:   VecSet(y, 0);
258:   if (!shell->leftwork2) VecDuplicate(y, &shell->leftwork2);
259:   y2 = shell->leftwork2;

261:   if (shell->Mvctx) { /* Have a merged Mvctx */
262:     /* 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
263:        in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an oppertunity to
264:        overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
265:      */
266:     VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD);
267:     VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD);

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

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

294:   if (shell->left) VecPointwiseMult(y, shell->left, y);
295:   VecScale(y, shell->scale);
296:   return 0;
297: }

299: PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y)
300: {
301:   Mat_Composite    *shell = (Mat_Composite *)A->data;
302:   Mat_CompositeLink next  = shell->head;
303:   Vec               in, y2 = NULL;
304:   PetscInt          i;

307:   in = x;
308:   if (shell->left) {
309:     if (!shell->leftwork) VecDuplicate(shell->left, &shell->leftwork);
310:     VecPointwiseMult(shell->leftwork, shell->left, in);
311:     in = shell->leftwork;
312:   }

314:   MatMultTranspose(next->mat, in, y);
315:   if (shell->scalings) {
316:     VecScale(y, shell->scalings[0]);
317:     if (!shell->rightwork2) VecDuplicate(y, &shell->rightwork2);
318:     y2 = shell->rightwork2;
319:   }
320:   i = 1;
321:   while ((next = next->next)) {
322:     if (!shell->scalings) MatMultTransposeAdd(next->mat, in, y, y);
323:     else {
324:       MatMultTranspose(next->mat, in, y2);
325:       VecAXPY(y, shell->scalings[i++], y2);
326:     }
327:   }
328:   if (shell->right) VecPointwiseMult(y, shell->right, y);
329:   VecScale(y, shell->scale);
330:   return 0;
331: }

333: PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z)
334: {
335:   Mat_Composite *shell = (Mat_Composite *)A->data;

337:   if (y != z) {
338:     MatMult(A, x, z);
339:     VecAXPY(z, 1.0, y);
340:   } else {
341:     if (!shell->leftwork) VecDuplicate(z, &shell->leftwork);
342:     MatMult(A, x, shell->leftwork);
343:     VecCopy(y, z);
344:     VecAXPY(z, 1.0, shell->leftwork);
345:   }
346:   return 0;
347: }

349: PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z)
350: {
351:   Mat_Composite *shell = (Mat_Composite *)A->data;

353:   if (y != z) {
354:     MatMultTranspose(A, x, z);
355:     VecAXPY(z, 1.0, y);
356:   } else {
357:     if (!shell->rightwork) VecDuplicate(z, &shell->rightwork);
358:     MatMultTranspose(A, x, shell->rightwork);
359:     VecCopy(y, z);
360:     VecAXPY(z, 1.0, shell->rightwork);
361:   }
362:   return 0;
363: }

365: PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v)
366: {
367:   Mat_Composite    *shell = (Mat_Composite *)A->data;
368:   Mat_CompositeLink next  = shell->head;
369:   PetscInt          i;


374:   MatGetDiagonal(next->mat, v);
375:   if (shell->scalings) VecScale(v, shell->scalings[0]);

377:   if (next->next && !shell->work) VecDuplicate(v, &shell->work);
378:   i = 1;
379:   while ((next = next->next)) {
380:     MatGetDiagonal(next->mat, shell->work);
381:     VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work);
382:   }
383:   VecScale(v, shell->scale);
384:   return 0;
385: }

387: PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t)
388: {
389:   Mat_Composite *shell = (Mat_Composite *)Y->data;

391:   if (shell->merge) MatCompositeMerge(Y);
392:   return 0;
393: }

395: PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha)
396: {
397:   Mat_Composite *a = (Mat_Composite *)inA->data;

399:   a->scale *= alpha;
400:   return 0;
401: }

403: PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right)
404: {
405:   Mat_Composite *a = (Mat_Composite *)inA->data;

407:   if (left) {
408:     if (!a->left) {
409:       VecDuplicate(left, &a->left);
410:       VecCopy(left, a->left);
411:     } else {
412:       VecPointwiseMult(a->left, left, a->left);
413:     }
414:   }
415:   if (right) {
416:     if (!a->right) {
417:       VecDuplicate(right, &a->right);
418:       VecCopy(right, a->right);
419:     } else {
420:       VecPointwiseMult(a->right, right, a->right);
421:     }
422:   }
423:   return 0;
424: }

426: PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject)
427: {
428:   Mat_Composite *a = (Mat_Composite *)A->data;

430:   PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options");
431:   PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL);
432:   PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL);
433:   PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL);
434:   PetscOptionsHeadEnd();
435:   return 0;
436: }

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

441:   Collective

443:    Input Parameters:
444: +  comm - MPI communicator
445: .  nmat - number of matrices to put in
446: -  mats - the matrices

448:    Output Parameter:
449: .  mat - the matrix

451:    Options Database Keys:
452: +  -mat_composite_merge         - merge in `MatAssemblyEnd()`
453: .  -mat_composite_merge_mvctx   - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices
454: -  -mat_composite_merge_type    - set merge direction

456:    Level: advanced

458:    Note:
459:      Alternative construction
460: $       MatCreate(comm,&mat);
461: $       MatSetSizes(mat,m,n,M,N);
462: $       MatSetType(mat,MATCOMPOSITE);
463: $       MatCompositeAddMat(mat,mats[0]);
464: $       ....
465: $       MatCompositeAddMat(mat,mats[nmat-1]);
466: $       MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
467: $       MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);

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

471: .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`, `MATCOMPOSITE`
472: @*/
473: PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat)
474: {
475:   PetscInt m, n, M, N, i;


480:   MatGetLocalSize(mats[0], PETSC_IGNORE, &n);
481:   MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE);
482:   MatGetSize(mats[0], PETSC_IGNORE, &N);
483:   MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE);
484:   MatCreate(comm, mat);
485:   MatSetSizes(*mat, m, n, M, N);
486:   MatSetType(*mat, MATCOMPOSITE);
487:   for (i = 0; i < nmat; i++) MatCompositeAddMat(*mat, mats[i]);
488:   MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY);
489:   MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY);
490:   return 0;
491: }

493: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat)
494: {
495:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
496:   Mat_CompositeLink ilink, next = shell->head;

498:   PetscNew(&ilink);
499:   ilink->next = NULL;
500:   PetscObjectReference((PetscObject)smat);
501:   ilink->mat = smat;

503:   if (!next) shell->head = ilink;
504:   else {
505:     while (next->next) next = next->next;
506:     next->next  = ilink;
507:     ilink->prev = next;
508:   }
509:   shell->tail = ilink;
510:   shell->nmat += 1;

512:   /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
513:   if (shell->scalings) {
514:     PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings);
515:     shell->scalings[shell->nmat - 1] = 1.0;
516:   }
517:   return 0;
518: }

520: /*@
521:     MatCompositeAddMat - Add another matrix to a composite matrix.

523:    Collective on mat

525:     Input Parameters:
526: +   mat - the composite matrix
527: -   smat - the partial matrix

529:    Level: advanced

531: .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
532: @*/
533: PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat)
534: {
537:   PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat));
538:   return 0;
539: }

541: static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type)
542: {
543:   Mat_Composite *b = (Mat_Composite *)mat->data;

545:   b->type = type;
546:   if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
547:     mat->ops->getdiagonal   = NULL;
548:     mat->ops->mult          = MatMult_Composite_Multiplicative;
549:     mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
550:     b->merge_mvctx          = PETSC_FALSE;
551:   } else {
552:     mat->ops->getdiagonal   = MatGetDiagonal_Composite;
553:     mat->ops->mult          = MatMult_Composite;
554:     mat->ops->multtranspose = MatMultTranspose_Composite;
555:   }
556:   return 0;
557: }

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

562:    Logically Collective on mat

564:    Input Parameters:
565: .  mat - the composite matrix

567:    Level: advanced

569: .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE`
570: @*/
571: PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type)
572: {
575:   PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type));
576:   return 0;
577: }

579: static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type)
580: {
581:   Mat_Composite *b = (Mat_Composite *)mat->data;

583:   *type = b->type;
584:   return 0;
585: }

587: /*@
588:    MatCompositeGetType - Returns type of composite.

590:    Not Collective

592:    Input Parameter:
593: .  mat - the composite matrix

595:    Output Parameter:
596: .  type - type of composite

598:    Level: advanced

600: .seealso: `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType`
601: @*/
602: PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type)
603: {
606:   PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type));
607:   return 0;
608: }

610: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str)
611: {
612:   Mat_Composite *b = (Mat_Composite *)mat->data;

614:   b->structure = str;
615:   return 0;
616: }

618: /*@
619:    MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.

621:    Not Collective

623:    Input Parameters:
624: +  mat - the composite matrix
625: -  str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN`

627:    Level: advanced

629:    Note:
630:     Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix.

632: .seealso: `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE`
633: @*/
634: PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str)
635: {
637:   PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str));
638:   return 0;
639: }

641: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str)
642: {
643:   Mat_Composite *b = (Mat_Composite *)mat->data;

645:   *str = b->structure;
646:   return 0;
647: }

649: /*@
650:    MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.

652:    Not Collective

654:    Input Parameter:
655: .  mat - the composite matrix

657:    Output Parameter:
658: .  str - structure of the matrices

660:    Level: advanced

662: .seealso: `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE`
663: @*/
664: PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str)
665: {
668:   PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str));
669:   return 0;
670: }

672: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type)
673: {
674:   Mat_Composite *shell = (Mat_Composite *)mat->data;

676:   shell->mergetype = type;
677:   return 0;
678: }

680: /*@
681:    MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`.

683:    Logically Collective on mat

685:    Input Parameters:
686: +  mat - the composite matrix
687: -  type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]),
688:           `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1])

690:    Level: advanced

692:    Note:
693:     The resulting matrix is the same regardles of the `MatCompositeMergeType`. Only the order of operation is changed.
694:     If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
695:     otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].

697: .seealso: `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE`
698: @*/
699: PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type)
700: {
703:   PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type));
704:   return 0;
705: }

707: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
708: {
709:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
710:   Mat_CompositeLink next = shell->head, prev = shell->tail;
711:   Mat               tmat, newmat;
712:   Vec               left, right;
713:   PetscScalar       scale;
714:   PetscInt          i;

717:   scale = shell->scale;
718:   if (shell->type == MAT_COMPOSITE_ADDITIVE) {
719:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
720:       i = 0;
721:       MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat);
722:       if (shell->scalings) MatScale(tmat, shell->scalings[i++]);
723:       while ((next = next->next)) MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure);
724:     } else {
725:       i = shell->nmat - 1;
726:       MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat);
727:       if (shell->scalings) MatScale(tmat, shell->scalings[i--]);
728:       while ((prev = prev->prev)) MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure);
729:     }
730:   } else {
731:     if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
732:       MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat);
733:       while ((next = next->next)) {
734:         MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat);
735:         MatDestroy(&tmat);
736:         tmat = newmat;
737:       }
738:     } else {
739:       MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat);
740:       while ((prev = prev->prev)) {
741:         MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat);
742:         MatDestroy(&tmat);
743:         tmat = newmat;
744:       }
745:     }
746:     if (shell->scalings) {
747:       for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
748:     }
749:   }

751:   if ((left = shell->left)) PetscObjectReference((PetscObject)left);
752:   if ((right = shell->right)) PetscObjectReference((PetscObject)right);

754:   MatHeaderReplace(mat, &tmat);

756:   MatDiagonalScale(mat, left, right);
757:   MatScale(mat, scale);
758:   VecDestroy(&left);
759:   VecDestroy(&right);
760:   return 0;
761: }

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

767:   Collective

769:    Input Parameter:
770: .  mat - the composite matrix

772:    Options Database Keys:
773: +  -mat_composite_merge - merge in `MatAssemblyEnd()`
774: -  -mat_composite_merge_type - set merge direction

776:    Level: advanced

778:    Note:
779:       The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix.

781: .seealso: `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE`
782: @*/
783: PetscErrorCode MatCompositeMerge(Mat mat)
784: {
786:   PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat));
787:   return 0;
788: }

790: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat)
791: {
792:   Mat_Composite *shell = (Mat_Composite *)mat->data;

794:   *nmat = shell->nmat;
795:   return 0;
796: }

798: /*@
799:    MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.

801:    Not Collective

803:    Input Parameter:
804: .  mat - the composite matrix

806:    Output Parameter:
807: .  nmat - number of matrices in the composite matrix

809:    Level: advanced

811: .seealso: `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
812: @*/
813: PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat)
814: {
817:   PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat));
818:   return 0;
819: }

821: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai)
822: {
823:   Mat_Composite    *shell = (Mat_Composite *)mat->data;
824:   Mat_CompositeLink ilink;
825:   PetscInt          k;

828:   ilink = shell->head;
829:   for (k = 0; k < i; k++) ilink = ilink->next;
830:   *Ai = ilink->mat;
831:   return 0;
832: }

834: /*@
835:    MatCompositeGetMat - Returns the ith matrix from the composite matrix.

837:    Logically Collective on mat

839:    Input Parameters:
840: +  mat - the composite matrix
841: -  i - the number of requested matrix

843:    Output Parameter:
844: .  Ai - ith matrix in composite

846:    Level: advanced

848: .seealso: `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE`
849: @*/
850: PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai)
851: {
855:   PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai));
856:   return 0;
857: }

859: PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings)
860: {
861:   Mat_Composite *shell = (Mat_Composite *)mat->data;
862:   PetscInt       nmat;

864:   MatCompositeGetNumberMat(mat, &nmat);
865:   if (!shell->scalings) PetscMalloc1(nmat, &shell->scalings);
866:   PetscArraycpy(shell->scalings, scalings, nmat);
867:   return 0;
868: }

870: /*@
871:    MatCompositeSetScalings - Sets separate scaling factors for component matrices.

873:    Logically Collective on mat

875:    Input Parameters:
876: +  mat      - the composite matrix
877: -  scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)

879:    Level: advanced

881: .seealso: `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE`
882: @*/
883: PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings)
884: {
888:   PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings));
889:   return 0;
890: }

892: static struct _MatOps MatOps_Values = {NULL,
893:                                        NULL,
894:                                        NULL,
895:                                        MatMult_Composite,
896:                                        MatMultAdd_Composite,
897:                                        /*  5*/ MatMultTranspose_Composite,
898:                                        MatMultTransposeAdd_Composite,
899:                                        NULL,
900:                                        NULL,
901:                                        NULL,
902:                                        /* 10*/ NULL,
903:                                        NULL,
904:                                        NULL,
905:                                        NULL,
906:                                        NULL,
907:                                        /* 15*/ NULL,
908:                                        NULL,
909:                                        MatGetDiagonal_Composite,
910:                                        MatDiagonalScale_Composite,
911:                                        NULL,
912:                                        /* 20*/ NULL,
913:                                        MatAssemblyEnd_Composite,
914:                                        NULL,
915:                                        NULL,
916:                                        /* 24*/ NULL,
917:                                        NULL,
918:                                        NULL,
919:                                        NULL,
920:                                        NULL,
921:                                        /* 29*/ NULL,
922:                                        NULL,
923:                                        NULL,
924:                                        NULL,
925:                                        NULL,
926:                                        /* 34*/ NULL,
927:                                        NULL,
928:                                        NULL,
929:                                        NULL,
930:                                        NULL,
931:                                        /* 39*/ NULL,
932:                                        NULL,
933:                                        NULL,
934:                                        NULL,
935:                                        NULL,
936:                                        /* 44*/ NULL,
937:                                        MatScale_Composite,
938:                                        MatShift_Basic,
939:                                        NULL,
940:                                        NULL,
941:                                        /* 49*/ NULL,
942:                                        NULL,
943:                                        NULL,
944:                                        NULL,
945:                                        NULL,
946:                                        /* 54*/ NULL,
947:                                        NULL,
948:                                        NULL,
949:                                        NULL,
950:                                        NULL,
951:                                        /* 59*/ NULL,
952:                                        MatDestroy_Composite,
953:                                        NULL,
954:                                        NULL,
955:                                        NULL,
956:                                        /* 64*/ NULL,
957:                                        NULL,
958:                                        NULL,
959:                                        NULL,
960:                                        NULL,
961:                                        /* 69*/ NULL,
962:                                        NULL,
963:                                        NULL,
964:                                        NULL,
965:                                        NULL,
966:                                        /* 74*/ NULL,
967:                                        NULL,
968:                                        MatSetFromOptions_Composite,
969:                                        NULL,
970:                                        NULL,
971:                                        /* 79*/ NULL,
972:                                        NULL,
973:                                        NULL,
974:                                        NULL,
975:                                        NULL,
976:                                        /* 84*/ NULL,
977:                                        NULL,
978:                                        NULL,
979:                                        NULL,
980:                                        NULL,
981:                                        /* 89*/ NULL,
982:                                        NULL,
983:                                        NULL,
984:                                        NULL,
985:                                        NULL,
986:                                        /* 94*/ NULL,
987:                                        NULL,
988:                                        NULL,
989:                                        NULL,
990:                                        NULL,
991:                                        /*99*/ NULL,
992:                                        NULL,
993:                                        NULL,
994:                                        NULL,
995:                                        NULL,
996:                                        /*104*/ NULL,
997:                                        NULL,
998:                                        NULL,
999:                                        NULL,
1000:                                        NULL,
1001:                                        /*109*/ NULL,
1002:                                        NULL,
1003:                                        NULL,
1004:                                        NULL,
1005:                                        NULL,
1006:                                        /*114*/ NULL,
1007:                                        NULL,
1008:                                        NULL,
1009:                                        NULL,
1010:                                        NULL,
1011:                                        /*119*/ NULL,
1012:                                        NULL,
1013:                                        NULL,
1014:                                        NULL,
1015:                                        NULL,
1016:                                        /*124*/ NULL,
1017:                                        NULL,
1018:                                        NULL,
1019:                                        NULL,
1020:                                        NULL,
1021:                                        /*129*/ NULL,
1022:                                        NULL,
1023:                                        NULL,
1024:                                        NULL,
1025:                                        NULL,
1026:                                        /*134*/ NULL,
1027:                                        NULL,
1028:                                        NULL,
1029:                                        NULL,
1030:                                        NULL,
1031:                                        /*139*/ NULL,
1032:                                        NULL,
1033:                                        NULL,
1034:                                        NULL,
1035:                                        NULL,
1036:                                        /*144*/ NULL,
1037:                                        NULL,
1038:                                        NULL,
1039:                                        NULL,
1040:                                        NULL,
1041:                                        NULL,
1042:                                        /*150*/ NULL};

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

1048:    Note:
1049:    To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`);

1051:   Level: advanced

1053: .seealso: `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`,
1054:           `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()`
1055: M*/

1057: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1058: {
1059:   Mat_Composite *b;

1061:   PetscNew(&b);
1062:   A->data = (void *)b;
1063:   PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps));

1065:   PetscLayoutSetUp(A->rmap);
1066:   PetscLayoutSetUp(A->cmap);

1068:   A->assembled    = PETSC_TRUE;
1069:   A->preallocated = PETSC_TRUE;
1070:   b->type         = MAT_COMPOSITE_ADDITIVE;
1071:   b->scale        = 1.0;
1072:   b->nmat         = 0;
1073:   b->merge        = PETSC_FALSE;
1074:   b->mergetype    = MAT_COMPOSITE_MERGE_RIGHT;
1075:   b->structure    = DIFFERENT_NONZERO_PATTERN;
1076:   b->merge_mvctx  = PETSC_TRUE;

1078:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite);
1079:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite);
1080:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite);
1081:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite);
1082:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite);
1083:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite);
1084:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite);
1085:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite);
1086:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite);
1087:   PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite);

1089:   PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE);
1090:   return 0;
1091: }