Actual source code: mcomposite.c
petsc-3.9.4 2018-09-11
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: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
387: Level: advanced
389: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
390: M*/
392: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
393: {
394: Mat_Composite *b;
398: PetscNewLog(A,&b);
399: A->data = (void*)b;
400: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
402: PetscLayoutSetUp(A->rmap);
403: PetscLayoutSetUp(A->cmap);
405: A->assembled = PETSC_TRUE;
406: A->preallocated = PETSC_TRUE;
407: b->type = MAT_COMPOSITE_ADDITIVE;
408: b->scale = 1.0;
409: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
410: return(0);
411: }
413: /*@
414: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
416: Collective on MPI_Comm
418: Input Parameters:
419: + comm - MPI communicator
420: . nmat - number of matrices to put in
421: - mats - the matrices
423: Output Parameter:
424: . mat - the matrix
426: Level: advanced
428: Notes:
429: Alternative construction
430: $ MatCreate(comm,&mat);
431: $ MatSetSizes(mat,m,n,M,N);
432: $ MatSetType(mat,MATCOMPOSITE);
433: $ MatCompositeAddMat(mat,mats[0]);
434: $ ....
435: $ MatCompositeAddMat(mat,mats[nmat-1]);
436: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
437: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
439: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
441: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
443: @*/
444: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
445: {
447: PetscInt m,n,M,N,i;
450: if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
453: MatGetLocalSize(mats[0],&m,&n);
454: MatGetSize(mats[0],&M,&N);
455: MatCreate(comm,mat);
456: MatSetSizes(*mat,m,n,M,N);
457: MatSetType(*mat,MATCOMPOSITE);
458: for (i=0; i<nmat; i++) {
459: MatCompositeAddMat(*mat,mats[i]);
460: }
461: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
462: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
463: return(0);
464: }
466: /*@
467: MatCompositeAddMat - add another matrix to a composite matrix
469: Collective on Mat
471: Input Parameters:
472: + mat - the composite matrix
473: - smat - the partial matrix
475: Level: advanced
477: .seealso: MatCreateComposite()
478: @*/
479: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
480: {
481: Mat_Composite *shell;
482: PetscErrorCode ierr;
483: Mat_CompositeLink ilink,next;
488: PetscNewLog(mat,&ilink);
489: ilink->next = 0;
490: PetscObjectReference((PetscObject)smat);
491: ilink->mat = smat;
493: shell = (Mat_Composite*)mat->data;
494: next = shell->head;
495: if (!next) shell->head = ilink;
496: else {
497: while (next->next) {
498: next = next->next;
499: }
500: next->next = ilink;
501: ilink->prev = next;
502: }
503: shell->tail = ilink;
504: return(0);
505: }
507: /*@
508: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
510: Collective on MPI_Comm
512: Input Parameters:
513: . mat - the composite matrix
516: Level: advanced
518: Notes:
519: The MatType of the resulting matrix will be the same as the MatType of the FIRST
520: matrix in the composite matrix.
522: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
524: @*/
525: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
526: {
527: Mat_Composite *b = (Mat_Composite*)mat->data;
528: PetscBool flg;
532: PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
533: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
534: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
535: mat->ops->getdiagonal = 0;
536: mat->ops->mult = MatMult_Composite_Multiplicative;
537: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
538: b->type = MAT_COMPOSITE_MULTIPLICATIVE;
539: } else {
540: mat->ops->getdiagonal = MatGetDiagonal_Composite;
541: mat->ops->mult = MatMult_Composite;
542: mat->ops->multtranspose = MatMultTranspose_Composite;
543: b->type = MAT_COMPOSITE_ADDITIVE;
544: }
545: return(0);
546: }
549: /*@
550: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
551: by summing all the matrices inside the composite matrix.
553: Collective on MPI_Comm
555: Input Parameters:
556: . mat - the composite matrix
559: Options Database:
560: . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
562: Level: advanced
564: Notes:
565: The MatType of the resulting matrix will be the same as the MatType of the FIRST
566: matrix in the composite matrix.
568: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
570: @*/
571: PetscErrorCode MatCompositeMerge(Mat mat)
572: {
573: Mat_Composite *shell = (Mat_Composite*)mat->data;
574: Mat_CompositeLink next = shell->head, prev = shell->tail;
575: PetscErrorCode ierr;
576: Mat tmat,newmat;
577: Vec left,right;
578: PetscScalar scale;
581: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
584: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
585: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
586: while ((next = next->next)) {
587: MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
588: }
589: } else {
590: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
591: while ((prev = prev->prev)) {
592: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
593: MatDestroy(&tmat);
594: tmat = newmat;
595: }
596: }
598: scale = shell->scale;
599: if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
600: if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}
602: MatHeaderReplace(mat,&tmat);
604: MatDiagonalScale(mat,left,right);
605: MatScale(mat,scale);
606: VecDestroy(&left);
607: VecDestroy(&right);
608: return(0);
609: }