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: PetscFunctionBegin;
42: while (next) {
43: PetscCall(MatDestroy(&next->mat));
44: if (next->work && (!next->next || next->work != next->next->work)) PetscCall(VecDestroy(&next->work));
45: oldnext = next;
46: next = next->next;
47: PetscCall(PetscFree(oldnext));
48: }
49: PetscCall(VecDestroy(&shell->work));
50: PetscCall(VecDestroy(&shell->left));
51: PetscCall(VecDestroy(&shell->right));
52: PetscCall(VecDestroy(&shell->leftwork));
53: PetscCall(VecDestroy(&shell->rightwork));
54: PetscCall(VecDestroy(&shell->leftwork2));
55: PetscCall(VecDestroy(&shell->rightwork2));
57: if (shell->Mvctx) {
58: for (i = 0; i < shell->nmat; i++) PetscCall(VecDestroy(&shell->lvecs[i]));
59: PetscCall(PetscFree3(shell->location, shell->larray, shell->lvecs));
60: PetscCall(PetscFree(shell->larray));
61: PetscCall(VecDestroy(&shell->gvec));
62: PetscCall(VecScatterDestroy(&shell->Mvctx));
63: }
65: PetscCall(PetscFree(shell->scalings));
66: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeAddMat_C", NULL));
67: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetType_C", NULL));
68: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetType_C", NULL));
69: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMergeType_C", NULL));
70: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetMatStructure_C", NULL));
71: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMatStructure_C", NULL));
72: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeMerge_C", NULL));
73: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetNumberMat_C", NULL));
74: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeGetMat_C", NULL));
75: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatCompositeSetScalings_C", NULL));
76: PetscCall(PetscFree(mat->data));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: PetscErrorCode MatMult_Composite_Multiplicative(Mat A, Vec x, Vec y)
81: {
82: Mat_Composite *shell = (Mat_Composite *)A->data;
83: Mat_CompositeLink next = shell->head;
84: Vec in, out;
85: PetscScalar scale;
86: PetscInt i;
88: PetscFunctionBegin;
89: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
90: in = x;
91: if (shell->right) {
92: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
93: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
94: in = shell->rightwork;
95: }
96: while (next->next) {
97: if (!next->work) { /* should reuse previous work if the same size */
98: PetscCall(MatCreateVecs(next->mat, NULL, &next->work));
99: }
100: out = next->work;
101: PetscCall(MatMult(next->mat, in, out));
102: in = out;
103: next = next->next;
104: }
105: PetscCall(MatMult(next->mat, in, y));
106: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
107: scale = shell->scale;
108: if (shell->scalings) {
109: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
110: }
111: PetscCall(VecScale(y, scale));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A, Vec x, Vec y)
116: {
117: Mat_Composite *shell = (Mat_Composite *)A->data;
118: Mat_CompositeLink tail = shell->tail;
119: Vec in, out;
120: PetscScalar scale;
121: PetscInt i;
123: PetscFunctionBegin;
124: PetscCheck(tail, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
125: in = x;
126: if (shell->left) {
127: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
128: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
129: in = shell->leftwork;
130: }
131: while (tail->prev) {
132: if (!tail->prev->work) { /* should reuse previous work if the same size */
133: PetscCall(MatCreateVecs(tail->mat, NULL, &tail->prev->work));
134: }
135: out = tail->prev->work;
136: PetscCall(MatMultTranspose(tail->mat, in, out));
137: in = out;
138: tail = tail->prev;
139: }
140: PetscCall(MatMultTranspose(tail->mat, in, y));
141: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
143: scale = shell->scale;
144: if (shell->scalings) {
145: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
146: }
147: PetscCall(VecScale(y, scale));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode MatMult_Composite(Mat mat, Vec x, Vec y)
152: {
153: Mat_Composite *shell = (Mat_Composite *)mat->data;
154: Mat_CompositeLink cur = shell->head;
155: Vec in, y2, xin;
156: Mat A, B;
157: PetscInt i, j, k, n, nuniq, lo, hi, mid, *gindices, *buf, *tmp, tot;
158: const PetscScalar *vals;
159: const PetscInt *garray;
160: IS ix, iy;
161: PetscBool match;
163: PetscFunctionBegin;
164: PetscCheck(cur, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
165: in = x;
166: if (shell->right) {
167: if (!shell->rightwork) PetscCall(VecDuplicate(shell->right, &shell->rightwork));
168: PetscCall(VecPointwiseMult(shell->rightwork, shell->right, in));
169: in = shell->rightwork;
170: }
172: /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time
173: we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and
174: it is legal to merge Mvctx, because all component matrices have the same size.
175: */
176: if (shell->merge_mvctx && !shell->Mvctx) {
177: /* Currently only implemented for MATMPIAIJ */
178: for (cur = shell->head; cur; cur = cur->next) {
179: PetscCall(PetscObjectTypeCompare((PetscObject)cur->mat, MATMPIAIJ, &match));
180: if (!match) {
181: shell->merge_mvctx = PETSC_FALSE;
182: goto skip_merge_mvctx;
183: }
184: }
186: /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */
187: tot = 0;
188: for (cur = shell->head; cur; cur = cur->next) {
189: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, NULL));
190: PetscCall(MatGetLocalSize(B, NULL, &n));
191: tot += n;
192: }
193: PetscCall(PetscMalloc3(tot, &shell->location, tot, &shell->larray, shell->nmat, &shell->lvecs));
194: shell->len = tot;
196: /* Go through matrices second time to sort off-diag columns and remove dups */
197: PetscCall(PetscMalloc1(tot, &gindices)); /* No Malloc2() since we will give one to petsc and free the other */
198: PetscCall(PetscMalloc1(tot, &buf));
199: nuniq = 0; /* Number of unique nonzero columns */
200: for (cur = shell->head; cur; cur = cur->next) {
201: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
202: PetscCall(MatGetLocalSize(B, NULL, &n));
203: /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */
204: i = j = k = 0;
205: while (i < n && j < nuniq) {
206: if (garray[i] < gindices[j]) buf[k++] = garray[i++];
207: else if (garray[i] > gindices[j]) buf[k++] = gindices[j++];
208: else {
209: buf[k++] = garray[i++];
210: j++;
211: }
212: }
213: /* Copy leftover in garray[] or gindices[] */
214: if (i < n) {
215: PetscCall(PetscArraycpy(buf + k, garray + i, n - i));
216: nuniq = k + n - i;
217: } else if (j < nuniq) {
218: PetscCall(PetscArraycpy(buf + k, gindices + j, nuniq - j));
219: nuniq = k + nuniq - j;
220: } else nuniq = k;
221: /* Swap gindices and buf to merge garray of the next matrix */
222: tmp = gindices;
223: gindices = buf;
224: buf = tmp;
225: }
226: PetscCall(PetscFree(buf));
228: /* Go through matrices third time to build a map from gindices[] to garray[] */
229: tot = 0;
230: for (cur = shell->head, j = 0; cur; cur = cur->next, j++) { /* j-th matrix */
231: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, NULL, &B, &garray));
232: PetscCall(MatGetLocalSize(B, NULL, &n));
233: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, n, NULL, &shell->lvecs[j]));
234: /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */
235: lo = 0;
236: for (i = 0; i < n; i++) {
237: hi = nuniq;
238: while (hi - lo > 1) {
239: mid = lo + (hi - lo) / 2;
240: if (garray[i] < gindices[mid]) hi = mid;
241: else lo = mid;
242: }
243: shell->location[tot + i] = lo; /* gindices[lo] = garray[i] */
244: lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */
245: }
246: tot += n;
247: }
249: /* Build merged Mvctx */
250: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nuniq, gindices, PETSC_OWN_POINTER, &ix));
251: PetscCall(ISCreateStride(PETSC_COMM_SELF, nuniq, 0, 1, &iy));
252: PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &xin));
253: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nuniq, &shell->gvec));
254: PetscCall(VecScatterCreate(xin, ix, shell->gvec, iy, &shell->Mvctx));
255: PetscCall(VecDestroy(&xin));
256: PetscCall(ISDestroy(&ix));
257: PetscCall(ISDestroy(&iy));
258: }
260: skip_merge_mvctx:
261: PetscCall(VecSet(y, 0));
262: if (!shell->leftwork2) PetscCall(VecDuplicate(y, &shell->leftwork2));
263: y2 = shell->leftwork2;
265: if (shell->Mvctx) { /* Have a merged Mvctx */
266: /* 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
267: in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an opportunity to
268: overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
269: */
270: PetscCall(VecScatterBegin(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
271: PetscCall(VecScatterEnd(shell->Mvctx, in, shell->gvec, INSERT_VALUES, SCATTER_FORWARD));
273: PetscCall(VecGetArrayRead(shell->gvec, &vals));
274: for (i = 0; i < shell->len; i++) shell->larray[i] = vals[shell->location[i]];
275: PetscCall(VecRestoreArrayRead(shell->gvec, &vals));
277: for (cur = shell->head, tot = i = 0; cur; cur = cur->next, i++) { /* i-th matrix */
278: PetscCall(MatMPIAIJGetSeqAIJ(cur->mat, &A, &B, NULL));
279: PetscUseTypeMethod(A, mult, in, y2);
280: PetscCall(MatGetLocalSize(B, NULL, &n));
281: PetscCall(VecPlaceArray(shell->lvecs[i], &shell->larray[tot]));
282: PetscCall((*B->ops->multadd)(B, shell->lvecs[i], y2, y2));
283: PetscCall(VecResetArray(shell->lvecs[i]));
284: PetscCall(VecAXPY(y, (shell->scalings ? shell->scalings[i] : 1.0), y2));
285: tot += n;
286: }
287: } else {
288: if (shell->scalings) {
289: for (cur = shell->head, i = 0; cur; cur = cur->next, i++) {
290: PetscCall(MatMult(cur->mat, in, y2));
291: PetscCall(VecAXPY(y, shell->scalings[i], y2));
292: }
293: } else {
294: for (cur = shell->head; cur; cur = cur->next) PetscCall(MatMultAdd(cur->mat, in, y, y));
295: }
296: }
298: if (shell->left) PetscCall(VecPointwiseMult(y, shell->left, y));
299: PetscCall(VecScale(y, shell->scale));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: PetscErrorCode MatMultTranspose_Composite(Mat A, Vec x, Vec y)
304: {
305: Mat_Composite *shell = (Mat_Composite *)A->data;
306: Mat_CompositeLink next = shell->head;
307: Vec in, y2 = NULL;
308: PetscInt i;
310: PetscFunctionBegin;
311: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
312: in = x;
313: if (shell->left) {
314: if (!shell->leftwork) PetscCall(VecDuplicate(shell->left, &shell->leftwork));
315: PetscCall(VecPointwiseMult(shell->leftwork, shell->left, in));
316: in = shell->leftwork;
317: }
319: PetscCall(MatMultTranspose(next->mat, in, y));
320: if (shell->scalings) {
321: PetscCall(VecScale(y, shell->scalings[0]));
322: if (!shell->rightwork2) PetscCall(VecDuplicate(y, &shell->rightwork2));
323: y2 = shell->rightwork2;
324: }
325: i = 1;
326: while ((next = next->next)) {
327: if (!shell->scalings) PetscCall(MatMultTransposeAdd(next->mat, in, y, y));
328: else {
329: PetscCall(MatMultTranspose(next->mat, in, y2));
330: PetscCall(VecAXPY(y, shell->scalings[i++], y2));
331: }
332: }
333: if (shell->right) PetscCall(VecPointwiseMult(y, shell->right, y));
334: PetscCall(VecScale(y, shell->scale));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscErrorCode MatMultAdd_Composite(Mat A, Vec x, Vec y, Vec z)
339: {
340: Mat_Composite *shell = (Mat_Composite *)A->data;
342: PetscFunctionBegin;
343: if (y != z) {
344: PetscCall(MatMult(A, x, z));
345: PetscCall(VecAXPY(z, 1.0, y));
346: } else {
347: if (!shell->leftwork) PetscCall(VecDuplicate(z, &shell->leftwork));
348: PetscCall(MatMult(A, x, shell->leftwork));
349: PetscCall(VecCopy(y, z));
350: PetscCall(VecAXPY(z, 1.0, shell->leftwork));
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: PetscErrorCode MatMultTransposeAdd_Composite(Mat A, Vec x, Vec y, Vec z)
356: {
357: Mat_Composite *shell = (Mat_Composite *)A->data;
359: PetscFunctionBegin;
360: if (y != z) {
361: PetscCall(MatMultTranspose(A, x, z));
362: PetscCall(VecAXPY(z, 1.0, y));
363: } else {
364: if (!shell->rightwork) PetscCall(VecDuplicate(z, &shell->rightwork));
365: PetscCall(MatMultTranspose(A, x, shell->rightwork));
366: PetscCall(VecCopy(y, z));
367: PetscCall(VecAXPY(z, 1.0, shell->rightwork));
368: }
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: PetscErrorCode MatGetDiagonal_Composite(Mat A, Vec v)
373: {
374: Mat_Composite *shell = (Mat_Composite *)A->data;
375: Mat_CompositeLink next = shell->head;
376: PetscInt i;
378: PetscFunctionBegin;
379: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
380: PetscCheck(!shell->right && !shell->left, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot get diagonal if left or right scaling");
382: PetscCall(MatGetDiagonal(next->mat, v));
383: if (shell->scalings) PetscCall(VecScale(v, shell->scalings[0]));
385: if (next->next && !shell->work) PetscCall(VecDuplicate(v, &shell->work));
386: i = 1;
387: while ((next = next->next)) {
388: PetscCall(MatGetDiagonal(next->mat, shell->work));
389: PetscCall(VecAXPY(v, (shell->scalings ? shell->scalings[i++] : 1.0), shell->work));
390: }
391: PetscCall(VecScale(v, shell->scale));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: PetscErrorCode MatAssemblyEnd_Composite(Mat Y, MatAssemblyType t)
396: {
397: Mat_Composite *shell = (Mat_Composite *)Y->data;
399: PetscFunctionBegin;
400: if (shell->merge) PetscCall(MatCompositeMerge(Y));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: PetscErrorCode MatScale_Composite(Mat inA, PetscScalar alpha)
405: {
406: Mat_Composite *a = (Mat_Composite *)inA->data;
408: PetscFunctionBegin;
409: a->scale *= alpha;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: PetscErrorCode MatDiagonalScale_Composite(Mat inA, Vec left, Vec right)
414: {
415: Mat_Composite *a = (Mat_Composite *)inA->data;
417: PetscFunctionBegin;
418: if (left) {
419: if (!a->left) {
420: PetscCall(VecDuplicate(left, &a->left));
421: PetscCall(VecCopy(left, a->left));
422: } else {
423: PetscCall(VecPointwiseMult(a->left, left, a->left));
424: }
425: }
426: if (right) {
427: if (!a->right) {
428: PetscCall(VecDuplicate(right, &a->right));
429: PetscCall(VecCopy(right, a->right));
430: } else {
431: PetscCall(VecPointwiseMult(a->right, right, a->right));
432: }
433: }
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: PetscErrorCode MatSetFromOptions_Composite(Mat A, PetscOptionItems *PetscOptionsObject)
438: {
439: Mat_Composite *a = (Mat_Composite *)A->data;
441: PetscFunctionBegin;
442: PetscOptionsHeadBegin(PetscOptionsObject, "MATCOMPOSITE options");
443: PetscCall(PetscOptionsBool("-mat_composite_merge", "Merge at MatAssemblyEnd", "MatCompositeMerge", a->merge, &a->merge, NULL));
444: PetscCall(PetscOptionsEnum("-mat_composite_merge_type", "Set composite merge direction", "MatCompositeSetMergeType", MatCompositeMergeTypes, (PetscEnum)a->mergetype, (PetscEnum *)&a->mergetype, NULL));
445: PetscCall(PetscOptionsBool("-mat_composite_merge_mvctx", "Merge MatMult() vecscat contexts", "MatCreateComposite", a->merge_mvctx, &a->merge_mvctx, NULL));
446: PetscOptionsHeadEnd();
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: MatCreateComposite - Creates a matrix as the sum or product of one or more matrices
453: Collective
455: Input Parameters:
456: + comm - MPI communicator
457: . nmat - number of matrices to put in
458: - mats - the matrices
460: Output Parameter:
461: . mat - the matrix
463: Options Database Keys:
464: + -mat_composite_merge - merge in `MatAssemblyEnd()`
465: . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in `MatMult()` for ADDITIVE matrices
466: - -mat_composite_merge_type - set merge direction
468: Level: advanced
470: Note:
471: Alternative construction
472: .vb
473: MatCreate(comm,&mat);
474: MatSetSizes(mat,m,n,M,N);
475: MatSetType(mat,MATCOMPOSITE);
476: MatCompositeAddMat(mat,mats[0]);
477: ....
478: MatCompositeAddMat(mat,mats[nmat-1]);
479: MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
480: MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
481: .ve
483: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
485: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCompositeGetMat()`, `MatCompositeMerge()`, `MatCompositeSetType()`,
486: `MATCOMPOSITE`, `MatCompositeType`
487: @*/
488: PetscErrorCode MatCreateComposite(MPI_Comm comm, PetscInt nmat, const Mat *mats, Mat *mat)
489: {
490: PetscInt m, n, M, N, i;
492: PetscFunctionBegin;
493: PetscCheck(nmat >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must pass in at least one matrix");
496: PetscCall(MatGetLocalSize(mats[0], PETSC_IGNORE, &n));
497: PetscCall(MatGetLocalSize(mats[nmat - 1], &m, PETSC_IGNORE));
498: PetscCall(MatGetSize(mats[0], PETSC_IGNORE, &N));
499: PetscCall(MatGetSize(mats[nmat - 1], &M, PETSC_IGNORE));
500: PetscCall(MatCreate(comm, mat));
501: PetscCall(MatSetSizes(*mat, m, n, M, N));
502: PetscCall(MatSetType(*mat, MATCOMPOSITE));
503: for (i = 0; i < nmat; i++) PetscCall(MatCompositeAddMat(*mat, mats[i]));
504: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
505: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
506: PetscFunctionReturn(PETSC_SUCCESS);
507: }
509: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat, Mat smat)
510: {
511: Mat_Composite *shell = (Mat_Composite *)mat->data;
512: Mat_CompositeLink ilink, next = shell->head;
514: PetscFunctionBegin;
515: PetscCall(PetscNew(&ilink));
516: ilink->next = NULL;
517: PetscCall(PetscObjectReference((PetscObject)smat));
518: ilink->mat = smat;
520: if (!next) shell->head = ilink;
521: else {
522: while (next->next) next = next->next;
523: next->next = ilink;
524: ilink->prev = next;
525: }
526: shell->tail = ilink;
527: shell->nmat += 1;
529: /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
530: if (shell->scalings) {
531: PetscCall(PetscRealloc(sizeof(PetscScalar) * shell->nmat, &shell->scalings));
532: shell->scalings[shell->nmat - 1] = 1.0;
533: }
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: MatCompositeAddMat - Add another matrix to a composite matrix.
540: Collective
542: Input Parameters:
543: + mat - the composite matrix
544: - smat - the partial matrix
546: Level: advanced
548: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
549: @*/
550: PetscErrorCode MatCompositeAddMat(Mat mat, Mat smat)
551: {
552: PetscFunctionBegin;
555: PetscUseMethod(mat, "MatCompositeAddMat_C", (Mat, Mat), (mat, smat));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: static PetscErrorCode MatCompositeSetType_Composite(Mat mat, MatCompositeType type)
560: {
561: Mat_Composite *b = (Mat_Composite *)mat->data;
563: PetscFunctionBegin;
564: b->type = type;
565: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
566: mat->ops->getdiagonal = NULL;
567: mat->ops->mult = MatMult_Composite_Multiplicative;
568: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
569: b->merge_mvctx = PETSC_FALSE;
570: } else {
571: mat->ops->getdiagonal = MatGetDiagonal_Composite;
572: mat->ops->mult = MatMult_Composite;
573: mat->ops->multtranspose = MatMultTranspose_Composite;
574: }
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: /*@
579: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product.
581: Logically Collective
583: Input Parameters:
584: + mat - the composite matrix
585: - type - the `MatCompositeType` to use for the matrix
587: Level: advanced
589: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeGetType()`, `MATCOMPOSITE`,
590: `MatCompositeType`
591: @*/
592: PetscErrorCode MatCompositeSetType(Mat mat, MatCompositeType type)
593: {
594: PetscFunctionBegin;
597: PetscUseMethod(mat, "MatCompositeSetType_C", (Mat, MatCompositeType), (mat, type));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: static PetscErrorCode MatCompositeGetType_Composite(Mat mat, MatCompositeType *type)
602: {
603: Mat_Composite *b = (Mat_Composite *)mat->data;
605: PetscFunctionBegin;
606: *type = b->type;
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: /*@
611: MatCompositeGetType - Returns type of composite.
613: Not Collective
615: Input Parameter:
616: . mat - the composite matrix
618: Output Parameter:
619: . type - type of composite
621: Level: advanced
623: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetType()`, `MATCOMPOSITE`, `MatCompositeType`
624: @*/
625: PetscErrorCode MatCompositeGetType(Mat mat, MatCompositeType *type)
626: {
627: PetscFunctionBegin;
630: PetscUseMethod(mat, "MatCompositeGetType_C", (Mat, MatCompositeType *), (mat, type));
631: PetscFunctionReturn(PETSC_SUCCESS);
632: }
634: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat, MatStructure str)
635: {
636: Mat_Composite *b = (Mat_Composite *)mat->data;
638: PetscFunctionBegin;
639: b->structure = str;
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: /*@
644: MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.
646: Not Collective
648: Input Parameters:
649: + mat - the composite matrix
650: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN` (default) or `SUBSET_NONZERO_PATTERN`
652: Level: advanced
654: Note:
655: Information about the matrices structure is used in `MatCompositeMerge()` for additive composite matrix.
657: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`, `MatCreateComposite()`, `MatCompositeMerge()` `MatCompositeGetMatStructure()`, `MATCOMPOSITE`
658: @*/
659: PetscErrorCode MatCompositeSetMatStructure(Mat mat, MatStructure str)
660: {
661: PetscFunctionBegin;
663: PetscUseMethod(mat, "MatCompositeSetMatStructure_C", (Mat, MatStructure), (mat, str));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat, MatStructure *str)
668: {
669: Mat_Composite *b = (Mat_Composite *)mat->data;
671: PetscFunctionBegin;
672: *str = b->structure;
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: /*@
677: MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.
679: Not Collective
681: Input Parameter:
682: . mat - the composite matrix
684: Output Parameter:
685: . str - structure of the matrices
687: Level: advanced
689: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MATCOMPOSITE`
690: @*/
691: PetscErrorCode MatCompositeGetMatStructure(Mat mat, MatStructure *str)
692: {
693: PetscFunctionBegin;
696: PetscUseMethod(mat, "MatCompositeGetMatStructure_C", (Mat, MatStructure *), (mat, str));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat, MatCompositeMergeType type)
701: {
702: Mat_Composite *shell = (Mat_Composite *)mat->data;
704: PetscFunctionBegin;
705: shell->mergetype = type;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /*@
710: MatCompositeSetMergeType - Sets order of `MatCompositeMerge()`.
712: Logically Collective
714: Input Parameters:
715: + mat - the composite matrix
716: - type - `MAT_COMPOSITE_MERGE RIGHT` (default) to start merge from right with the first added matrix (mat[0]),
717: `MAT_COMPOSITE_MERGE_LEFT` to start merge from left with the last added matrix (mat[nmat-1])
719: Level: advanced
721: Note:
722: The resulting matrix is the same regardless of the `MatCompositeMergeType`. Only the order of operation is changed.
723: If set to `MAT_COMPOSITE_MERGE_RIGHT` the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
724: otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].
726: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeMerge()`, `MATCOMPOSITE`
727: @*/
728: PetscErrorCode MatCompositeSetMergeType(Mat mat, MatCompositeMergeType type)
729: {
730: PetscFunctionBegin;
733: PetscUseMethod(mat, "MatCompositeSetMergeType_C", (Mat, MatCompositeMergeType), (mat, type));
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
738: {
739: Mat_Composite *shell = (Mat_Composite *)mat->data;
740: Mat_CompositeLink next = shell->head, prev = shell->tail;
741: Mat tmat, newmat;
742: Vec left, right;
743: PetscScalar scale;
744: PetscInt i;
746: PetscFunctionBegin;
747: PetscCheck(next, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must provide at least one matrix with MatCompositeAddMat()");
748: scale = shell->scale;
749: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
750: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
751: i = 0;
752: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
753: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i++]));
754: while ((next = next->next)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i++] : 1.0), next->mat, shell->structure));
755: } else {
756: i = shell->nmat - 1;
757: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
758: if (shell->scalings) PetscCall(MatScale(tmat, shell->scalings[i--]));
759: while ((prev = prev->prev)) PetscCall(MatAXPY(tmat, (shell->scalings ? shell->scalings[i--] : 1.0), prev->mat, shell->structure));
760: }
761: } else {
762: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
763: PetscCall(MatDuplicate(next->mat, MAT_COPY_VALUES, &tmat));
764: while ((next = next->next)) {
765: PetscCall(MatMatMult(next->mat, tmat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
766: PetscCall(MatDestroy(&tmat));
767: tmat = newmat;
768: }
769: } else {
770: PetscCall(MatDuplicate(prev->mat, MAT_COPY_VALUES, &tmat));
771: while ((prev = prev->prev)) {
772: PetscCall(MatMatMult(tmat, prev->mat, MAT_INITIAL_MATRIX, PETSC_DECIDE, &newmat));
773: PetscCall(MatDestroy(&tmat));
774: tmat = newmat;
775: }
776: }
777: if (shell->scalings) {
778: for (i = 0; i < shell->nmat; i++) scale *= shell->scalings[i];
779: }
780: }
782: if ((left = shell->left)) PetscCall(PetscObjectReference((PetscObject)left));
783: if ((right = shell->right)) PetscCall(PetscObjectReference((PetscObject)right));
785: PetscCall(MatHeaderReplace(mat, &tmat));
787: PetscCall(MatDiagonalScale(mat, left, right));
788: PetscCall(MatScale(mat, scale));
789: PetscCall(VecDestroy(&left));
790: PetscCall(VecDestroy(&right));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@
795: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
796: by summing or computing the product of all the matrices inside the composite matrix.
798: Collective
800: Input Parameter:
801: . mat - the composite matrix
803: Options Database Keys:
804: + -mat_composite_merge - merge in `MatAssemblyEnd()`
805: - -mat_composite_merge_type - set merge direction
807: Level: advanced
809: Note:
810: The `MatType` of the resulting matrix will be the same as the `MatType` of the FIRST matrix in the composite matrix.
812: .seealso: [](ch_matrices), `Mat`, `MatDestroy()`, `MatMult()`, `MatCompositeAddMat()`, `MatCreateComposite()`, `MatCompositeSetMatStructure()`, `MatCompositeSetMergeType()`, `MATCOMPOSITE`
813: @*/
814: PetscErrorCode MatCompositeMerge(Mat mat)
815: {
816: PetscFunctionBegin;
818: PetscUseMethod(mat, "MatCompositeMerge_C", (Mat), (mat));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat, PetscInt *nmat)
823: {
824: Mat_Composite *shell = (Mat_Composite *)mat->data;
826: PetscFunctionBegin;
827: *nmat = shell->nmat;
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: /*@
832: MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.
834: Not Collective
836: Input Parameter:
837: . mat - the composite matrix
839: Output Parameter:
840: . nmat - number of matrices in the composite matrix
842: Level: advanced
844: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetMat()`, `MATCOMPOSITE`
845: @*/
846: PetscErrorCode MatCompositeGetNumberMat(Mat mat, PetscInt *nmat)
847: {
848: PetscFunctionBegin;
851: PetscUseMethod(mat, "MatCompositeGetNumberMat_C", (Mat, PetscInt *), (mat, nmat));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat, PetscInt i, Mat *Ai)
856: {
857: Mat_Composite *shell = (Mat_Composite *)mat->data;
858: Mat_CompositeLink ilink;
859: PetscInt k;
861: PetscFunctionBegin;
862: PetscCheck(i < shell->nmat, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_OUTOFRANGE, "index out of range: %" PetscInt_FMT " >= %" PetscInt_FMT, i, shell->nmat);
863: ilink = shell->head;
864: for (k = 0; k < i; k++) ilink = ilink->next;
865: *Ai = ilink->mat;
866: PetscFunctionReturn(PETSC_SUCCESS);
867: }
869: /*@
870: MatCompositeGetMat - Returns the ith matrix from the composite matrix.
872: Logically Collective
874: Input Parameters:
875: + mat - the composite matrix
876: - i - the number of requested matrix
878: Output Parameter:
879: . Ai - ith matrix in composite
881: Level: advanced
883: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeGetNumberMat()`, `MatCompositeAddMat()`, `MATCOMPOSITE`
884: @*/
885: PetscErrorCode MatCompositeGetMat(Mat mat, PetscInt i, Mat *Ai)
886: {
887: PetscFunctionBegin;
891: PetscUseMethod(mat, "MatCompositeGetMat_C", (Mat, PetscInt, Mat *), (mat, i, Ai));
892: PetscFunctionReturn(PETSC_SUCCESS);
893: }
895: PetscErrorCode MatCompositeSetScalings_Composite(Mat mat, const PetscScalar *scalings)
896: {
897: Mat_Composite *shell = (Mat_Composite *)mat->data;
898: PetscInt nmat;
900: PetscFunctionBegin;
901: PetscCall(MatCompositeGetNumberMat(mat, &nmat));
902: if (!shell->scalings) PetscCall(PetscMalloc1(nmat, &shell->scalings));
903: PetscCall(PetscArraycpy(shell->scalings, scalings, nmat));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: /*@
908: MatCompositeSetScalings - Sets separate scaling factors for component matrices.
910: Logically Collective
912: Input Parameters:
913: + mat - the composite matrix
914: - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)
916: Level: advanced
918: .seealso: [](ch_matrices), `Mat`, `MatScale()`, `MatDiagonalScale()`, `MATCOMPOSITE`
919: @*/
920: PetscErrorCode MatCompositeSetScalings(Mat mat, const PetscScalar *scalings)
921: {
922: PetscFunctionBegin;
926: PetscUseMethod(mat, "MatCompositeSetScalings_C", (Mat, const PetscScalar *), (mat, scalings));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static struct _MatOps MatOps_Values = {NULL,
931: NULL,
932: NULL,
933: MatMult_Composite,
934: MatMultAdd_Composite,
935: /* 5*/ MatMultTranspose_Composite,
936: MatMultTransposeAdd_Composite,
937: NULL,
938: NULL,
939: NULL,
940: /* 10*/ NULL,
941: NULL,
942: NULL,
943: NULL,
944: NULL,
945: /* 15*/ NULL,
946: NULL,
947: MatGetDiagonal_Composite,
948: MatDiagonalScale_Composite,
949: NULL,
950: /* 20*/ NULL,
951: MatAssemblyEnd_Composite,
952: NULL,
953: NULL,
954: /* 24*/ NULL,
955: NULL,
956: NULL,
957: NULL,
958: NULL,
959: /* 29*/ NULL,
960: NULL,
961: NULL,
962: NULL,
963: NULL,
964: /* 34*/ NULL,
965: NULL,
966: NULL,
967: NULL,
968: NULL,
969: /* 39*/ NULL,
970: NULL,
971: NULL,
972: NULL,
973: NULL,
974: /* 44*/ NULL,
975: MatScale_Composite,
976: MatShift_Basic,
977: NULL,
978: NULL,
979: /* 49*/ NULL,
980: NULL,
981: NULL,
982: NULL,
983: NULL,
984: /* 54*/ NULL,
985: NULL,
986: NULL,
987: NULL,
988: NULL,
989: /* 59*/ NULL,
990: MatDestroy_Composite,
991: NULL,
992: NULL,
993: NULL,
994: /* 64*/ NULL,
995: NULL,
996: NULL,
997: NULL,
998: NULL,
999: /* 69*/ NULL,
1000: NULL,
1001: NULL,
1002: NULL,
1003: NULL,
1004: /* 74*/ NULL,
1005: NULL,
1006: MatSetFromOptions_Composite,
1007: NULL,
1008: NULL,
1009: /* 79*/ NULL,
1010: NULL,
1011: NULL,
1012: NULL,
1013: NULL,
1014: /* 84*/ NULL,
1015: NULL,
1016: NULL,
1017: NULL,
1018: NULL,
1019: /* 89*/ NULL,
1020: NULL,
1021: NULL,
1022: NULL,
1023: NULL,
1024: /* 94*/ NULL,
1025: NULL,
1026: NULL,
1027: NULL,
1028: NULL,
1029: /*99*/ NULL,
1030: NULL,
1031: NULL,
1032: NULL,
1033: NULL,
1034: /*104*/ NULL,
1035: NULL,
1036: NULL,
1037: NULL,
1038: NULL,
1039: /*109*/ NULL,
1040: NULL,
1041: NULL,
1042: NULL,
1043: NULL,
1044: /*114*/ NULL,
1045: NULL,
1046: NULL,
1047: NULL,
1048: NULL,
1049: /*119*/ NULL,
1050: NULL,
1051: NULL,
1052: NULL,
1053: NULL,
1054: /*124*/ NULL,
1055: NULL,
1056: NULL,
1057: NULL,
1058: NULL,
1059: /*129*/ NULL,
1060: NULL,
1061: NULL,
1062: NULL,
1063: NULL,
1064: /*134*/ NULL,
1065: NULL,
1066: NULL,
1067: NULL,
1068: NULL,
1069: /*139*/ NULL,
1070: NULL,
1071: NULL,
1072: NULL,
1073: NULL,
1074: /*144*/ NULL,
1075: NULL,
1076: NULL,
1077: NULL,
1078: NULL,
1079: NULL,
1080: /*150*/ NULL,
1081: NULL};
1083: /*MC
1084: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
1085: The matrices need to have a correct size and parallel layout for the sum or product to be valid.
1087: Level: advanced
1089: Note:
1090: To use the product of the matrices call `MatCompositeSetType`(mat,`MAT_COMPOSITE_MULTIPLICATIVE`);
1092: .seealso: [](ch_matrices), `Mat`, `MatCreateComposite()`, `MatCompositeSetScalings()`, `MatCompositeAddMat()`, `MatSetType()`, `MatCompositeSetType()`, `MatCompositeGetType()`,
1093: `MatCompositeSetMatStructure()`, `MatCompositeGetMatStructure()`, `MatCompositeMerge()`, `MatCompositeSetMergeType()`, `MatCompositeGetNumberMat()`, `MatCompositeGetMat()`
1094: M*/
1096: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1097: {
1098: Mat_Composite *b;
1100: PetscFunctionBegin;
1101: PetscCall(PetscNew(&b));
1102: A->data = (void *)b;
1103: PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps)));
1105: PetscCall(PetscLayoutSetUp(A->rmap));
1106: PetscCall(PetscLayoutSetUp(A->cmap));
1108: A->assembled = PETSC_TRUE;
1109: A->preallocated = PETSC_TRUE;
1110: b->type = MAT_COMPOSITE_ADDITIVE;
1111: b->scale = 1.0;
1112: b->nmat = 0;
1113: b->merge = PETSC_FALSE;
1114: b->mergetype = MAT_COMPOSITE_MERGE_RIGHT;
1115: b->structure = DIFFERENT_NONZERO_PATTERN;
1116: b->merge_mvctx = PETSC_TRUE;
1118: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeAddMat_C", MatCompositeAddMat_Composite));
1119: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetType_C", MatCompositeSetType_Composite));
1120: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetType_C", MatCompositeGetType_Composite));
1121: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMergeType_C", MatCompositeSetMergeType_Composite));
1122: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetMatStructure_C", MatCompositeSetMatStructure_Composite));
1123: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMatStructure_C", MatCompositeGetMatStructure_Composite));
1124: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeMerge_C", MatCompositeMerge_Composite));
1125: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetNumberMat_C", MatCompositeGetNumberMat_Composite));
1126: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeGetMat_C", MatCompositeGetMat_Composite));
1127: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatCompositeSetScalings_C", MatCompositeSetScalings_Composite));
1129: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATCOMPOSITE));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }