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