Actual source code: mcomposite.c
petsc-3.10.5 2019-03-28
2: #include <petsc/private/matimpl.h>
4: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
5: struct _Mat_CompositeLink {
6: Mat mat;
7: Vec work;
8: Mat_CompositeLink next,prev;
9: };
11: typedef struct {
12: MatCompositeType type;
13: Mat_CompositeLink head,tail;
14: Vec work;
15: PetscScalar scale; /* scale factor supplied with MatScale() */
16: Vec left,right; /* left and right diagonal scaling provided with MatDiagonalScale() */
17: Vec leftwork,rightwork;
18: } Mat_Composite;
20: PetscErrorCode MatDestroy_Composite(Mat mat)
21: {
22: PetscErrorCode ierr;
23: Mat_Composite *shell = (Mat_Composite*)mat->data;
24: Mat_CompositeLink next = shell->head,oldnext;
27: while (next) {
28: MatDestroy(&next->mat);
29: if (next->work && (!next->next || next->work != next->next->work)) {
30: VecDestroy(&next->work);
31: }
32: oldnext = next;
33: next = next->next;
34: PetscFree(oldnext);
35: }
36: VecDestroy(&shell->work);
37: VecDestroy(&shell->left);
38: VecDestroy(&shell->right);
39: VecDestroy(&shell->leftwork);
40: VecDestroy(&shell->rightwork);
41: PetscFree(mat->data);
42: return(0);
43: }
45: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
46: {
47: Mat_Composite *shell = (Mat_Composite*)A->data;
48: Mat_CompositeLink next = shell->head;
49: PetscErrorCode ierr;
50: Vec in,out;
53: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
54: in = x;
55: if (shell->right) {
56: if (!shell->rightwork) {
57: VecDuplicate(shell->right,&shell->rightwork);
58: }
59: VecPointwiseMult(shell->rightwork,shell->right,in);
60: in = shell->rightwork;
61: }
62: while (next->next) {
63: if (!next->work) { /* should reuse previous work if the same size */
64: MatCreateVecs(next->mat,NULL,&next->work);
65: }
66: out = next->work;
67: MatMult(next->mat,in,out);
68: in = out;
69: next = next->next;
70: }
71: MatMult(next->mat,in,y);
72: if (shell->left) {
73: VecPointwiseMult(y,shell->left,y);
74: }
75: VecScale(y,shell->scale);
76: return(0);
77: }
79: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
80: {
81: Mat_Composite *shell = (Mat_Composite*)A->data;
82: Mat_CompositeLink tail = shell->tail;
83: PetscErrorCode ierr;
84: Vec in,out;
87: if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
88: in = x;
89: if (shell->left) {
90: if (!shell->leftwork) {
91: VecDuplicate(shell->left,&shell->leftwork);
92: }
93: VecPointwiseMult(shell->leftwork,shell->left,in);
94: in = shell->leftwork;
95: }
96: while (tail->prev) {
97: if (!tail->prev->work) { /* should reuse previous work if the same size */
98: MatCreateVecs(tail->mat,NULL,&tail->prev->work);
99: }
100: out = tail->prev->work;
101: MatMultTranspose(tail->mat,in,out);
102: in = out;
103: tail = tail->prev;
104: }
105: MatMultTranspose(tail->mat,in,y);
106: if (shell->right) {
107: VecPointwiseMult(y,shell->right,y);
108: }
109: VecScale(y,shell->scale);
110: return(0);
111: }
113: PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
114: {
115: Mat_Composite *shell = (Mat_Composite*)A->data;
116: Mat_CompositeLink next = shell->head;
117: PetscErrorCode ierr;
118: Vec in;
121: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
122: in = x;
123: if (shell->right) {
124: if (!shell->rightwork) {
125: VecDuplicate(shell->right,&shell->rightwork);
126: }
127: VecPointwiseMult(shell->rightwork,shell->right,in);
128: in = shell->rightwork;
129: }
130: MatMult(next->mat,in,y);
131: while ((next = next->next)) {
132: MatMultAdd(next->mat,in,y,y);
133: }
134: if (shell->left) {
135: VecPointwiseMult(y,shell->left,y);
136: }
137: VecScale(y,shell->scale);
138: return(0);
139: }
141: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
142: {
143: Mat_Composite *shell = (Mat_Composite*)A->data;
144: Mat_CompositeLink next = shell->head;
145: PetscErrorCode ierr;
146: Vec in;
149: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
150: in = x;
151: if (shell->left) {
152: if (!shell->leftwork) {
153: VecDuplicate(shell->left,&shell->leftwork);
154: }
155: VecPointwiseMult(shell->leftwork,shell->left,in);
156: in = shell->leftwork;
157: }
158: MatMultTranspose(next->mat,in,y);
159: while ((next = next->next)) {
160: MatMultTransposeAdd(next->mat,in,y,y);
161: }
162: if (shell->right) {
163: VecPointwiseMult(y,shell->right,y);
164: }
165: VecScale(y,shell->scale);
166: return(0);
167: }
169: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
170: {
171: Mat_Composite *shell = (Mat_Composite*)A->data;
172: Mat_CompositeLink next = shell->head;
173: PetscErrorCode ierr;
176: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
177: if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling");
179: MatGetDiagonal(next->mat,v);
180: if (next->next && !shell->work) {
181: VecDuplicate(v,&shell->work);
182: }
183: while ((next = next->next)) {
184: MatGetDiagonal(next->mat,shell->work);
185: VecAXPY(v,1.0,shell->work);
186: }
187: VecScale(v,shell->scale);
188: return(0);
189: }
191: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
192: {
194: PetscBool flg = PETSC_FALSE;
197: PetscOptionsGetBool(((PetscObject)Y)->options,((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);
198: if (flg) {
199: MatCompositeMerge(Y);
200: }
201: return(0);
202: }
204: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
205: {
206: Mat_Composite *a = (Mat_Composite*)inA->data;
209: a->scale *= alpha;
210: return(0);
211: }
213: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
214: {
215: Mat_Composite *a = (Mat_Composite*)inA->data;
219: if (left) {
220: if (!a->left) {
221: VecDuplicate(left,&a->left);
222: VecCopy(left,a->left);
223: } else {
224: VecPointwiseMult(a->left,left,a->left);
225: }
226: }
227: if (right) {
228: if (!a->right) {
229: VecDuplicate(right,&a->right);
230: VecCopy(right,a->right);
231: } else {
232: VecPointwiseMult(a->right,right,a->right);
233: }
234: }
235: return(0);
236: }
238: static struct _MatOps MatOps_Values = {0,
239: 0,
240: 0,
241: MatMult_Composite,
242: 0,
243: /* 5*/ MatMultTranspose_Composite,
244: 0,
245: 0,
246: 0,
247: 0,
248: /* 10*/ 0,
249: 0,
250: 0,
251: 0,
252: 0,
253: /* 15*/ 0,
254: 0,
255: MatGetDiagonal_Composite,
256: MatDiagonalScale_Composite,
257: 0,
258: /* 20*/ 0,
259: MatAssemblyEnd_Composite,
260: 0,
261: 0,
262: /* 24*/ 0,
263: 0,
264: 0,
265: 0,
266: 0,
267: /* 29*/ 0,
268: 0,
269: 0,
270: 0,
271: 0,
272: /* 34*/ 0,
273: 0,
274: 0,
275: 0,
276: 0,
277: /* 39*/ 0,
278: 0,
279: 0,
280: 0,
281: 0,
282: /* 44*/ 0,
283: MatScale_Composite,
284: MatShift_Basic,
285: 0,
286: 0,
287: /* 49*/ 0,
288: 0,
289: 0,
290: 0,
291: 0,
292: /* 54*/ 0,
293: 0,
294: 0,
295: 0,
296: 0,
297: /* 59*/ 0,
298: MatDestroy_Composite,
299: 0,
300: 0,
301: 0,
302: /* 64*/ 0,
303: 0,
304: 0,
305: 0,
306: 0,
307: /* 69*/ 0,
308: 0,
309: 0,
310: 0,
311: 0,
312: /* 74*/ 0,
313: 0,
314: 0,
315: 0,
316: 0,
317: /* 79*/ 0,
318: 0,
319: 0,
320: 0,
321: 0,
322: /* 84*/ 0,
323: 0,
324: 0,
325: 0,
326: 0,
327: /* 89*/ 0,
328: 0,
329: 0,
330: 0,
331: 0,
332: /* 94*/ 0,
333: 0,
334: 0,
335: 0,
336: 0,
337: /*99*/ 0,
338: 0,
339: 0,
340: 0,
341: 0,
342: /*104*/ 0,
343: 0,
344: 0,
345: 0,
346: 0,
347: /*109*/ 0,
348: 0,
349: 0,
350: 0,
351: 0,
352: /*114*/ 0,
353: 0,
354: 0,
355: 0,
356: 0,
357: /*119*/ 0,
358: 0,
359: 0,
360: 0,
361: 0,
362: /*124*/ 0,
363: 0,
364: 0,
365: 0,
366: 0,
367: /*129*/ 0,
368: 0,
369: 0,
370: 0,
371: 0,
372: /*134*/ 0,
373: 0,
374: 0,
375: 0,
376: 0,
377: /*139*/ 0,
378: 0,
379: 0
380: };
382: /*MC
383: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout).
385: Notes:
386: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
388: Level: advanced
390: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
391: M*/
393: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
394: {
395: Mat_Composite *b;
399: PetscNewLog(A,&b);
400: A->data = (void*)b;
401: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
403: PetscLayoutSetUp(A->rmap);
404: PetscLayoutSetUp(A->cmap);
406: A->assembled = PETSC_TRUE;
407: A->preallocated = PETSC_TRUE;
408: b->type = MAT_COMPOSITE_ADDITIVE;
409: b->scale = 1.0;
410: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
411: return(0);
412: }
414: /*@
415: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
417: Collective on MPI_Comm
419: Input Parameters:
420: + comm - MPI communicator
421: . nmat - number of matrices to put in
422: - mats - the matrices
424: Output Parameter:
425: . mat - the matrix
427: Level: advanced
429: Notes:
430: Alternative construction
431: $ MatCreate(comm,&mat);
432: $ MatSetSizes(mat,m,n,M,N);
433: $ MatSetType(mat,MATCOMPOSITE);
434: $ MatCompositeAddMat(mat,mats[0]);
435: $ ....
436: $ MatCompositeAddMat(mat,mats[nmat-1]);
437: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
438: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
440: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
442: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
444: @*/
445: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
446: {
448: PetscInt m,n,M,N,i;
451: if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
454: MatGetLocalSize(mats[0],&m,&n);
455: MatGetSize(mats[0],&M,&N);
456: MatCreate(comm,mat);
457: MatSetSizes(*mat,m,n,M,N);
458: MatSetType(*mat,MATCOMPOSITE);
459: for (i=0; i<nmat; i++) {
460: MatCompositeAddMat(*mat,mats[i]);
461: }
462: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
464: return(0);
465: }
467: /*@
468: MatCompositeAddMat - add another matrix to a composite matrix
470: Collective on Mat
472: Input Parameters:
473: + mat - the composite matrix
474: - smat - the partial matrix
476: Level: advanced
478: .seealso: MatCreateComposite()
479: @*/
480: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
481: {
482: Mat_Composite *shell;
483: PetscErrorCode ierr;
484: Mat_CompositeLink ilink,next;
489: PetscNewLog(mat,&ilink);
490: ilink->next = 0;
491: PetscObjectReference((PetscObject)smat);
492: ilink->mat = smat;
494: shell = (Mat_Composite*)mat->data;
495: next = shell->head;
496: if (!next) shell->head = ilink;
497: else {
498: while (next->next) {
499: next = next->next;
500: }
501: next->next = ilink;
502: ilink->prev = next;
503: }
504: shell->tail = ilink;
505: return(0);
506: }
508: /*@
509: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
511: Collective on MPI_Comm
513: Input Parameters:
514: . mat - the composite matrix
517: Level: advanced
519: Notes:
520: The MatType of the resulting matrix will be the same as the MatType of the FIRST
521: matrix in the composite matrix.
523: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
525: @*/
526: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
527: {
528: Mat_Composite *b = (Mat_Composite*)mat->data;
529: PetscBool flg;
533: PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
534: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
535: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
536: mat->ops->getdiagonal = 0;
537: mat->ops->mult = MatMult_Composite_Multiplicative;
538: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
539: b->type = MAT_COMPOSITE_MULTIPLICATIVE;
540: } else {
541: mat->ops->getdiagonal = MatGetDiagonal_Composite;
542: mat->ops->mult = MatMult_Composite;
543: mat->ops->multtranspose = MatMultTranspose_Composite;
544: b->type = MAT_COMPOSITE_ADDITIVE;
545: }
546: return(0);
547: }
550: /*@
551: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
552: by summing all the matrices inside the composite matrix.
554: Collective on MPI_Comm
556: Input Parameters:
557: . mat - the composite matrix
560: Options Database:
561: . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
563: Level: advanced
565: Notes:
566: The MatType of the resulting matrix will be the same as the MatType of the FIRST
567: matrix in the composite matrix.
569: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
571: @*/
572: PetscErrorCode MatCompositeMerge(Mat mat)
573: {
574: Mat_Composite *shell = (Mat_Composite*)mat->data;
575: Mat_CompositeLink next = shell->head, prev = shell->tail;
576: PetscErrorCode ierr;
577: Mat tmat,newmat;
578: Vec left,right;
579: PetscScalar scale;
582: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
585: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
586: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
587: while ((next = next->next)) {
588: MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
589: }
590: } else {
591: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
592: while ((prev = prev->prev)) {
593: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
594: MatDestroy(&tmat);
595: tmat = newmat;
596: }
597: }
599: scale = shell->scale;
600: if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
601: if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}
603: MatHeaderReplace(mat,&tmat);
605: MatDiagonalScale(mat,left,right);
606: MatScale(mat,scale);
607: VecDestroy(&left);
608: VecDestroy(&right);
609: return(0);
610: }