Actual source code: mcomposite.c
1: #include <petsc/private/matimpl.h>
3: const char *const MatCompositeMergeTypes[] = {"left", "right", "MatCompositeMergeType", "MAT_COMPOSITE_", NULL};
5: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
6: struct _Mat_CompositeLink {
7: Mat mat;
8: Vec work;
9: Mat_CompositeLink next, prev;
10: };
12: typedef struct {
13: MatCompositeType type;
14: Mat_CompositeLink head, tail;
15: Vec work;
16: PetscScalar scale; /* scale factor supplied with MatScale() */
17: Vec left, right; /* left and right diagonal scaling provided with MatDiagonalScale() */
18: Vec leftwork, rightwork, leftwork2, rightwork2; /* Two pairs of working vectors */
19: PetscInt nmat;
20: PetscBool merge;
21: MatCompositeMergeType mergetype;
22: MatStructure structure;
24: PetscScalar *scalings;
25: PetscBool merge_mvctx; /* Whether need to merge mvctx of component matrices */
26: Vec *lvecs; /* [nmat] Basically, they are Mvctx->lvec of each component matrix */
27: PetscScalar *larray; /* [len] Data arrays of lvecs[] are stored consecutively in larray */
28: PetscInt len; /* Length of larray[] */
29: Vec gvec; /* Union of lvecs[] without duplicated entries */
30: PetscInt *location; /* A map that maps entries in garray[] to larray[] */
31: VecScatter Mvctx;
32: } Mat_Composite;
34: static PetscErrorCode MatDestroy_Composite(Mat mat)
35: {
36: Mat_Composite *shell = (Mat_Composite *)mat->data;
37: Mat_CompositeLink next = shell->head, oldnext;
38: PetscInt i;
40: PetscFunctionBegin;
41: while (next) {
42: PetscCall(MatDestroy(&next->mat));
43: if (next->work && (!next->next || next->work != next->next->work)) PetscCall(VecDestroy(&next->work));
44: oldnext = next;
45: next = next->next;
46: PetscCall(PetscFree(oldnext));
47: }
48: PetscCall(VecDestroy(&shell->work));
49: PetscCall(VecDestroy(&shell->left));
50: PetscCall(VecDestroy(&shell->right));
51: PetscCall(VecDestroy(&shell->leftwork));
52: PetscCall(VecDestroy(&shell->rightwork));
53: PetscCall(VecDestroy(&shell->leftwork2));
54: PetscCall(VecDestroy(&shell->rightwork2));
56: if (shell->Mvctx) {
57: for (i = 0; i < shell->nmat; i++) PetscCall(VecDestroy(&shell->lvecs[i]));
58: PetscCall(PetscFree3(shell->location, shell->larray, shell->lvecs));
59: PetscCall(PetscFree(shell->larray));
60: PetscCall(VecDestroy(&shell->gvec));
61: PetscCall(VecScatterDestroy(&shell->Mvctx));
62: }
64: PetscCall(PetscFree(shell->scalings));
65: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL));
66: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL));
67: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL));
68: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL));
69: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL));
70: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL));
71: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL));
72: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL));
73: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL));
74: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL));
75: PetscCall(PetscFree(mat->data));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static 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;
87: PetscFunctionBegin;
88: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
89: in = x;
90: if (shell->right) {
91: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
92: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
93: in = shell->rightwork;
94: }
95: while (next->next) {
96: if (!next->work) { /* should reuse previous work if the same size */
97: PetscCall(MatCreateVecs(next->mat, NULL, &next->work));
98: }
99: out = next->work;
100: PetscCall(MatMult(next->mat, in, out));
101: in = out;
102: next = next->next;
103: }
104: PetscCall(MatMult(next->mat, in, y));
105: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
106: scale = shell->scale;
107: if (shell->scalings) {
108: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
109: }
110: PetscCall(VecScale(y, scale));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: static PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A, Vec x, Vec y)
115: {
116: Mat_Composite *shell = (Mat_Composite *)A->data;
117: Mat_CompositeLink tail = shell->tail;
118: Vec in, out;
119: PetscScalar scale;
120: PetscInt i;
122: PetscFunctionBegin;
123: PetscCheck(tail, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
124: in = x;
125: if (shell->left) {
126: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
127: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
128: in = shell->leftwork;
129: }
130: while (tail->prev) {
131: if (!tail->prev->work) { /* should reuse previous work if the same size */
132: PetscCall(MatCreateVecs(tail->mat, NULL, &tail->prev->work));
133: }
134: out = tail->prev->work;
135: PetscCall(MatMultTranspose(tail->mat, in, out));
136: in = out;
137: tail = tail->prev;
138: }
139: PetscCall(MatMultTranspose(tail->mat, in, y));
140: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
142: scale = shell->scale;
143: if (shell->scalings) {
144: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
145: }
146: PetscCall(VecScale(y, scale));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode MatMult_Composite(Mat mat, Vec x, Vec y)
151: {
152: Mat_Composite *shell = (Mat_Composite *)mat->data;
153: Mat_CompositeLink cur = shell->head;
154: Vec in, y2, xin;
155: Mat A, B;
156: PetscInt i, j, k, n, nuniq, lo, hi, mid, *gindices, *buf, *tmp, tot;
157: const PetscScalar *vals;
158: const PetscInt *garray;
159: IS ix, iy;
160: PetscBool match;
162: PetscFunctionBegin;
163: PetscCheck(cur, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
164: in = x;
165: if (shell->right) {
166: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
167: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
168: in = shell->rightwork;
169: }
171: /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time
172: we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and
173: it is legal to merge Mvctx, because all component matrices have the same size.
174: */
175: if (shell->merge_mvctx && !shell->Mvctx) {
176: /* Currently only implemented for MATMPIAIJ */
177: for (cur = shell->head; cur; cur = cur->next) {
178: PetscCall(PetscObjectTypeCompare((PetscObject)cur->mat, MATMPIAIJ, &match));
179: if (!match) {
180: shell->merge_mvctx = PETSC_FALSE;
181: goto skip_merge_mvctx;
182: }
183: }
185: /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */
186: tot = 0;
187: for (cur = shell->head; cur; cur = cur->next) {
188: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, NULL));
189: PetscCall(MatGetLocalSize(B, NULL, &n));
190: tot += n;
191: }
192: PetscCall(PetscMalloc3(tot, &shell->location, tot, &shell->larray, shell->nmat, &shell->lvecs));
193: shell->len = tot;
195: /* Go through matrices second time to sort off-diag columns and remove dups */
196: PetscCall(PetscMalloc1(tot, &gindices)); /* No Malloc2() since we will give one to petsc and free the other */
197: PetscCall(PetscMalloc1(tot, &buf));
198: nuniq = 0; /* Number of unique nonzero columns */
199: for (cur = shell->head; cur; cur = cur->next) {
200: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
201: PetscCall(MatGetLocalSize(B, NULL, &n));
202: /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */
203: i = j = k = 0;
204: while (i < n && j < nuniq) {
205: if (garray[i] < gindices[j]) buf[k++] = garray[i++];
206: else if (garray[i] > gindices[j]) buf[k++] = gindices[j++];
207: else {
208: buf[k++] = garray[i++];
209: j++;
210: }
211: }
212: /* Copy leftover in garray[] or gindices[] */
213: if (i < n) {
214: PetscCall(PetscArraycpy(buf + k, garray + i, n - i));
215: nuniq = k + n - i;
216: } else if (j < nuniq) {
217: PetscCall(PetscArraycpy(buf + k, gindices + j, nuniq - j));
218: nuniq = k + nuniq - j;
219: } else nuniq = k;
220: /* Swap gindices and buf to merge garray of the next matrix */
221: tmp = gindices;
222: gindices = buf;
223: buf = tmp;
224: }
225: PetscCall(PetscFree(buf));
227: /* Go through matrices third time to build a map from gindices[] to garray[] */
228: tot = 0;
229: for (cur = shell->head, j = 0; cur; cur = cur->next, j++) { /* j-th matrix */
230: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
231: PetscCall(MatGetLocalSize(B, NULL, &n));
232: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &shell->lvecs[j]));
233: /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */
234: lo = 0;
235: for (i = 0; i < n; i++) {
236: hi = nuniq;
237: while (hi - lo > 1) {
238: mid = lo + (hi - lo) / 2;
239: if (garray[i] < gindices[mid]) hi = mid;
240: else lo = mid;
241: }
242: shell->location[tot + i] = lo; /* gindices[lo] = garray[i] */
243: lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */
244: }
245: tot += n;
246: }
248: /* Build merged Mvctx */
249: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nuniq, gindices, PETSC_OWN_POINTER, &ix));
250: PetscCall(ISCreateStride(PETSC_COMM_SELF, nuniq, 0, 1, &iy));
251: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &xin));
252: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nuniq, &shell->gvec));
253: PetscCall(VecScatterCreate(xin, ix, shell->gvec, iy, &shell->Mvctx));
254: PetscCall(VecDestroy(&xin));
255: PetscCall(ISDestroy(&ix));
256: PetscCall(ISDestroy(&iy));
257: }
259: skip_merge_mvctx:
260: PetscCall(VecSet(y, 0));
261: if (!shell->leftwork2) PetscCall(VecDuplicate(y, &shell->leftwork2));
262: y2 = shell->leftwork2;
264: if (shell->Mvctx) { /* Have a merged Mvctx */
265: /* 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
266: in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an opportunity to
267: overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
268: */
269: PetscCall(VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
270: PetscCall(VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
272: PetscCall(VecGetArrayRead(shell->gvec, &vals));
273: for (i = 0; i < shell->len; i++) shell->larray[i] = vals[shell->location[i]];
274: PetscCall(VecRestoreArrayRead(shell->gvec, &vals));
276: for (cur = shell->head, tot = i = 0; cur; cur = cur->next, i++) { /* i-th matrix */
277: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, &A, &B, NULL));
278: PetscUseTypeMethod(A, mult, in, y2);
279: PetscCall(MatGetLocalSize(B, NULL, &n));
280: PetscCall(VecPlaceArray(shell->lvecs[i], &shell->larray[tot]));
281: PetscCall((*B->ops->multadd)(B, shell->lvecs[i], y2, y2));
282: PetscCall(VecResetArray(shell->lvecs[i]));
283: PetscCall(VecAXPY(y, (shell->scalings ? shell->scalings[i] : 1.0), y2));
284: tot += n;
285: }
286: } else {
287: if (shell->scalings) {
288: for (cur = shell->head, i = 0; cur; cur = cur->next, i++) {
289: PetscCall(MatMult(cur->mat, in, y2));
290: PetscCall(VecAXPY(y, shell->scalings[i], y2));
291: }
292: } else {
293: for (cur = shell->head; cur; cur = cur->next) PetscCall(MatMultAdd(cur->mat, in, y, y));
294: }
295: }
297: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
298: PetscCall(VecScale(y, shell->scale));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: static PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y)
303: {
304: Mat_Composite *shell = (Mat_Composite *)A->data;
305: Mat_CompositeLink next = shell->head;
306: Vec in, y2 = NULL;
307: PetscInt i;
309: PetscFunctionBegin;
310: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
311: in = x;
312: if (shell->left) {
313: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
314: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
315: in = shell->leftwork;
316: }
318: PetscCall(MatMultTranspose(next->mat, in, y));
319: if (shell->scalings) {
320: PetscCall(VecScale(y, shell->scalings[0]));
321: if (!shell->rightwork2) PetscCall(VecDuplicate(y, &shell->rightwork2));
322: y2 = shell->rightwork2;
323: }
324: i = 1;
325: while ((next = next->next)) {
326: if (!shell->scalings) PetscCall(MatMultTransposeAdd(next->mat, in, y, y));
327: else {
328: PetscCall(MatMultTranspose(next->mat, in, y2));
329: PetscCall(VecAXPY(y, shell->scalings[i++], y2));
330: }
331: }
332: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
333: PetscCall(VecScale(y, shell->scale));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z)
338: {
339: Mat_Composite *shell = (Mat_Composite *)A->data;
341: PetscFunctionBegin;
342: if (y != z) {
343: PetscCall(MatMult(A, x, z));
344: PetscCall(VecAXPY(z, 1.0, y));
345: } else {
346: if (!shell->leftwork) PetscCall(VecDuplicate(z, &shell->leftwork));
347: PetscCall(MatMult(A, x, shell->leftwork));
348: PetscCall(VecCopy(y, z));
349: PetscCall(VecAXPY(z, 1.0, shell->leftwork));
350: }
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z)
355: {
356: Mat_Composite *shell = (Mat_Composite *)A->data;
358: PetscFunctionBegin;
359: if (y != z) {
360: PetscCall(MatMultTranspose(A, x, z));
361: PetscCall(VecAXPY(z, 1.0, y));
362: } else {
363: if (!shell->rightwork) PetscCall(VecDuplicate(z, &shell->rightwork));
364: PetscCall(MatMultTranspose(A, x, shell->rightwork));
365: PetscCall(VecCopy(y, z));
366: PetscCall(VecAXPY(z, 1.0, shell->rightwork));
367: }
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v)
372: {
373: Mat_Composite *shell = (Mat_Composite *)A->data;
374: Mat_CompositeLink next = shell->head;
375: PetscInt i;
377: PetscFunctionBegin;
378: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
379: PetscCheck(!shell->right && !shell->left, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get diagonal if left or right scaling");
381: PetscCall(MatGetDiagonal(next->mat, v));
382: if (shell->scalings) PetscCall(VecScale(v, shell->scalings[0]));
384: if (next->next && !shell->work) PetscCall(VecDuplicate(v, &shell->work));
385: i = 1;
386: while ((next = next->next)) {
387: PetscCall(MatGetDiagonal(next->mat, shell->work));
388: PetscCall(VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work));
389: }
390: PetscCall(VecScale(v, shell->scale));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t)
395: {
396: Mat_Composite *shell = (Mat_Composite *)Y->data;
398: PetscFunctionBegin;
399: if (shell->merge) PetscCall(MatCompositeMerge(Y));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha)
404: {
405: Mat_Composite *a = (Mat_Composite *)inA->data;
407: PetscFunctionBegin;
408: a->scale *= alpha;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: static PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right)
413: {
414: Mat_Composite *a = (Mat_Composite *)inA->data;
416: PetscFunctionBegin;
417: if (left) {
418: if (!a->left) {
419: PetscCall(VecDuplicate(left, &a->left));
420: PetscCall(VecCopy(left, a->left));
421: } else {
422: PetscCall(VecPointwiseMult(a->left, left, a->left));
423: }
424: }
425: if (right) {
426: if (!a->right) {
427: PetscCall(VecDuplicate(right, &a->right));
428: PetscCall(VecCopy(right, a->right));
429: } else {
430: PetscCall(VecPointwiseMult(a->right, right, a->right));
431: }
432: }
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject)
437: {
438: Mat_Composite *a = (Mat_Composite *)A->data;
440: PetscFunctionBegin;
441: PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options");
442: PetscCall(PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL));
443: PetscCall(PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL));
444: PetscCall(PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL));
445: PetscOptionsHeadEnd();
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*@
450: MatCreateComposite - Creates a matrix as the sum or product of one or more matrices
452: Collective
454: Input Parameters:
455: + comm - MPI communicator
456: . nmat - number of matrices to put in
457: - mats - the matrices
459: Output Parameter:
460: . mat - the matrix
462: Options Database Keys:
463: + -mat_composite_merge - merge in `MatAssemblyEnd()`
464: . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices
465: - -mat_composite_merge_type - set merge direction
467: Level: advanced
469: Note:
470: Alternative construction
471: .vb
472: MatCreate(comm,&mat);
473: MatSetSizes(mat,m,n,M,N);
474: MatSetType(mat,MATCOMPOSITE);
475: MatCompositeAddMat(mat,mats[0]);
476: ....
477: MatCompositeAddMat(mat,mats[nmat-1]);
478: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
479: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
480: .ve
482: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
484: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`,
485: `MATCOMPOSITE`, `MatCompositeType`
486: @*/
487: PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat)
488: {
489: PetscInt m, n, M, N, i;
491: PetscFunctionBegin;
492: PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix");
493: PetscAssertPointer(mat, 4);
495: PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n));
496: PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE));
497: PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N));
498: PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE));
499: PetscCall(MatCreate(comm, mat));
500: PetscCall(MatSetSizes(*mat, m, n, M, N));
501: PetscCall(MatSetType(*mat, MATCOMPOSITE));
502: for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i]));
503: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
504: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat)
509: {
510: Mat_Composite *shell = (Mat_Composite *)mat->data;
511: Mat_CompositeLink ilink, next = shell->head;
513: PetscFunctionBegin;
514: PetscCall(PetscNew(&ilink));
515: ilink->next = NULL;
516: PetscCall(PetscObjectReference((PetscObject)smat));
517: ilink->mat = smat;
519: if (!next) shell->head = ilink;
520: else {
521: while (next->next) next = next->next;
522: next->next = ilink;
523: ilink->prev = next;
524: }
525: shell->tail = ilink;
526: shell->nmat += 1;
528: /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
529: if (shell->scalings) {
530: PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings));
531: shell->scalings[shell->nmat - 1] = 1.0;
532: }
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: /*@
537: MatCompositeAddMat - Add another matrix to a composite matrix.
539: Collective
541: Input Parameters:
542: + mat - the composite matrix
543: - smat - the partial matrix
545: Level: advanced
547: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
548: @*/
549: PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat)
550: {
551: PetscFunctionBegin;
554: PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type)
559: {
560: Mat_Composite *b = (Mat_Composite *)mat->data;
562: PetscFunctionBegin;
563: b->type = type;
564: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
565: mat->ops->getdiagonal = NULL;
566: mat->ops->mult = MatMult_Composite_Multiplicative;
567: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
568: b->merge_mvctx = PETSC_FALSE;
569: } else {
570: mat->ops->getdiagonal = MatGetDiagonal_Composite;
571: mat->ops->mult = MatMult_Composite;
572: mat->ops->multtranspose = MatMultTranspose_Composite;
573: }
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product.
580: Logically Collective
582: Input Parameters:
583: + mat - the composite matrix
584: - type - the `MatCompositeType` to use for the matrix
586: Level: advanced
588: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE`,
589: `MatCompositeType`
590: @*/
591: PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type)
592: {
593: PetscFunctionBegin;
596: PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type)
601: {
602: Mat_Composite *b = (Mat_Composite *)mat->data;
604: PetscFunctionBegin;
605: *type = b->type;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: MatCompositeGetType - Returns type of composite.
612: Not Collective
614: Input Parameter:
615: . mat - the composite matrix
617: Output Parameter:
618: . type - type of composite
620: Level: advanced
622: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType`
623: @*/
624: PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type)
625: {
626: PetscFunctionBegin;
628: PetscAssertPointer(type, 2);
629: PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str)
634: {
635: Mat_Composite *b = (Mat_Composite *)mat->data;
637: PetscFunctionBegin;
638: b->structure = str;
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: /*@
643: MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.
645: Not Collective
647: Input Parameters:
648: + mat - the composite matrix
649: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN`
651: Level: advanced
653: Note:
654: Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix.
656: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE`
657: @*/
658: PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str)
659: {
660: PetscFunctionBegin;
662: PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str)
667: {
668: Mat_Composite *b = (Mat_Composite *)mat->data;
670: PetscFunctionBegin;
671: *str = b->structure;
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: /*@
676: MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.
678: Not Collective
680: Input Parameter:
681: . mat - the composite matrix
683: Output Parameter:
684: . str - structure of the matrices
686: Level: advanced
688: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE`
689: @*/
690: PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str)
691: {
692: PetscFunctionBegin;
694: PetscAssertPointer(str, 2);
695: PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type)
700: {
701: Mat_Composite *shell = (Mat_Composite *)mat->data;
703: PetscFunctionBegin;
704: shell->mergetype = type;
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*@
709: MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`.
711: Logically Collective
713: Input Parameters:
714: + mat - the composite matrix
715: - type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]),
716: `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1])
718: Level: advanced
720: Note:
721: The resulting matrix is the same regardless of the `MatCompositeMergeType`. Only the order of operation is changed.
722: If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
723: otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].
725: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE`
726: @*/
727: PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type)
728: {
729: PetscFunctionBegin;
732: PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
737: {
738: Mat_Composite *shell = (Mat_Composite *)mat->data;
739: Mat_CompositeLink next = shell->head, prev = shell->tail;
740: Mat tmat, newmat;
741: Vec left, right;
742: PetscScalar scale;
743: PetscInt i;
745: PetscFunctionBegin;
746: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
747: scale = shell->scale;
748: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
749: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
750: i = 0;
751: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
752: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++]));
753: while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure));
754: } else {
755: i = shell->nmat - 1;
756: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
757: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--]));
758: while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure));
759: }
760: } else {
761: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
762: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
763: while ((next = next->next)) {
764: PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
765: PetscCall(MatDestroy(&tmat));
766: tmat = newmat;
767: }
768: } else {
769: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
770: while ((prev = prev->prev)) {
771: PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
772: PetscCall(MatDestroy(&tmat));
773: tmat = newmat;
774: }
775: }
776: if (shell->scalings) {
777: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
778: }
779: }
781: if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left));
782: if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right));
784: PetscCall(MatHeaderReplace(mat, &tmat));
786: PetscCall(MatDiagonalScale(mat, left, right));
787: PetscCall(MatScale(mat, scale));
788: PetscCall(VecDestroy(&left));
789: PetscCall(VecDestroy(&right));
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@
794: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
795: by summing or computing the product of all the matrices inside the composite matrix.
797: Collective
799: Input Parameter:
800: . mat - the composite matrix
802: Options Database Keys:
803: + -mat_composite_merge - merge in `MatAssemblyEnd()`
804: - -mat_composite_merge_type - set merge direction
806: Level: advanced
808: Note:
809: The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix.
811: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE`
812: @*/
813: PetscErrorCode MatCompositeMerge(Mat mat)
814: {
815: PetscFunctionBegin;
817: PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat));
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat)
822: {
823: Mat_Composite *shell = (Mat_Composite *)mat->data;
825: PetscFunctionBegin;
826: *nmat = shell->nmat;
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: /*@
831: MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.
833: Not Collective
835: Input Parameter:
836: . mat - the composite matrix
838: Output Parameter:
839: . nmat - number of matrices in the composite matrix
841: Level: advanced
843: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
844: @*/
845: PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat)
846: {
847: PetscFunctionBegin;
849: PetscAssertPointer(nmat, 2);
850: PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat));
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai)
855: {
856: Mat_Composite *shell = (Mat_Composite *)mat->data;
857: Mat_CompositeLink ilink;
858: PetscInt k;
860: PetscFunctionBegin;
861: PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat);
862: ilink = shell->head;
863: for (k = 0; k < i; k++) ilink = ilink->next;
864: *Ai = ilink->mat;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@
869: MatCompositeGetMat - Returns the ith matrix from the composite matrix.
871: Logically Collective
873: Input Parameters:
874: + mat - the composite matrix
875: - i - the number of requested matrix
877: Output Parameter:
878: . Ai - ith matrix in composite
880: Level: advanced
882: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE`
883: @*/
884: PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai)
885: {
886: PetscFunctionBegin;
889: PetscAssertPointer(Ai, 3);
890: PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai));
891: PetscFunctionReturn(PETSC_SUCCESS);
892: }
894: static PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings)
895: {
896: Mat_Composite *shell = (Mat_Composite *)mat->data;
897: PetscInt nmat;
899: PetscFunctionBegin;
900: PetscCall(MatCompositeGetNumberMat(mat, &nmat));
901: if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings));
902: PetscCall(PetscArraycpy(shell->scalings, scalings, nmat));
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: /*@
907: MatCompositeSetScalings - Sets separate scaling factors for component matrices.
909: Logically Collective
911: Input Parameters:
912: + mat - the composite matrix
913: - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)
915: Level: advanced
917: .seealso: [](ch_matrices), `Mat`, `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE`
918: @*/
919: PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings)
920: {
921: PetscFunctionBegin;
923: PetscAssertPointer(scalings, 2);
925: PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: static struct _MatOps MatOps_Values = {NULL,
930: NULL,
931: NULL,
932: MatMult_Composite,
933: MatMultAdd_Composite,
934: /* 5*/ MatMultTranspose_Composite,
935: MatMultTransposeAdd_Composite,
936: NULL,
937: NULL,
938: NULL,
939: /* 10*/ NULL,
940: NULL,
941: NULL,
942: NULL,
943: NULL,
944: /* 15*/ NULL,
945: NULL,
946: MatGetDiagonal_Composite,
947: MatDiagonalScale_Composite,
948: NULL,
949: /* 20*/ NULL,
950: MatAssemblyEnd_Composite,
951: NULL,
952: NULL,
953: /* 24*/ NULL,
954: NULL,
955: NULL,
956: NULL,
957: NULL,
958: /* 29*/ NULL,
959: NULL,
960: NULL,
961: NULL,
962: NULL,
963: /* 34*/ NULL,
964: NULL,
965: NULL,
966: NULL,
967: NULL,
968: /* 39*/ NULL,
969: NULL,
970: NULL,
971: NULL,
972: NULL,
973: /* 44*/ NULL,
974: MatScale_Composite,
975: MatShift_Basic,
976: NULL,
977: NULL,
978: /* 49*/ NULL,
979: NULL,
980: NULL,
981: NULL,
982: NULL,
983: /* 54*/ NULL,
984: NULL,
985: NULL,
986: NULL,
987: NULL,
988: /* 59*/ NULL,
989: MatDestroy_Composite,
990: NULL,
991: NULL,
992: NULL,
993: /* 64*/ NULL,
994: NULL,
995: NULL,
996: NULL,
997: NULL,
998: /* 69*/ NULL,
999: NULL,
1000: NULL,
1001: NULL,
1002: NULL,
1003: /* 74*/ NULL,
1004: NULL,
1005: MatSetFromOptions_Composite,
1006: NULL,
1007: NULL,
1008: /* 79*/ NULL,
1009: NULL,
1010: NULL,
1011: NULL,
1012: NULL,
1013: /* 84*/ NULL,
1014: NULL,
1015: NULL,
1016: NULL,
1017: NULL,
1018: /* 89*/ NULL,
1019: NULL,
1020: NULL,
1021: NULL,
1022: NULL,
1023: /* 94*/ NULL,
1024: NULL,
1025: NULL,
1026: NULL,
1027: NULL,
1028: /*99*/ NULL,
1029: NULL,
1030: NULL,
1031: NULL,
1032: NULL,
1033: /*104*/ NULL,
1034: NULL,
1035: NULL,
1036: NULL,
1037: NULL,
1038: /*109*/ NULL,
1039: NULL,
1040: NULL,
1041: NULL,
1042: NULL,
1043: /*114*/ NULL,
1044: NULL,
1045: NULL,
1046: NULL,
1047: NULL,
1048: /*119*/ NULL,
1049: NULL,
1050: NULL,
1051: NULL,
1052: NULL,
1053: /*124*/ NULL,
1054: NULL,
1055: NULL,
1056: NULL,
1057: NULL,
1058: /*129*/ NULL,
1059: NULL,
1060: NULL,
1061: NULL,
1062: NULL,
1063: /*134*/ NULL,
1064: NULL,
1065: NULL,
1066: NULL,
1067: NULL,
1068: /*139*/ NULL,
1069: NULL,
1070: NULL,
1071: NULL,
1072: NULL,
1073: /*144*/ NULL,
1074: NULL,
1075: NULL,
1076: NULL,
1077: NULL,
1078: NULL,
1079: /*150*/ NULL,
1080: NULL};
1082: /*MC
1083: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
1084: The matrices need to have a correct size and parallel layout for the sum or product to be valid.
1086: Level: advanced
1088: Note:
1089: To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`);
1091: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`,
1092: `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()`
1093: M*/
1095: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1096: {
1097: Mat_Composite *b;
1099: PetscFunctionBegin;
1100: PetscCall(PetscNew(&b));
1101: A->data = (void *)b;
1102: A->ops[0] = MatOps_Values;
1104: PetscCall(PetscLayoutSetUp(A->rmap));
1105: PetscCall(PetscLayoutSetUp(A->cmap));
1107: A->assembled = PETSC_TRUE;
1108: A->preallocated = PETSC_TRUE;
1109: b->type = MAT_COMPOSITE_ADDITIVE;
1110: b->scale = 1.0;
1111: b->nmat = 0;
1112: b->merge = PETSC_FALSE;
1113: b->mergetype = MAT_COMPOSITE_MERGE_RIGHT;
1114: b->structure = DIFFERENT_NONZERO_PATTERN;
1115: b->merge_mvctx = PETSC_TRUE;
1117: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite));
1118: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite));
1119: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite));
1120: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite));
1121: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite));
1122: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite));
1123: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite));
1124: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite));
1125: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite));
1126: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite));
1128: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE));
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }