Actual source code: mcomposite.c
petsc-3.6.1 2015-08-06
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
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;
22: PetscErrorCode MatDestroy_Composite(Mat mat)
23: {
24: PetscErrorCode ierr;
25: Mat_Composite *shell = (Mat_Composite*)mat->data;
26: Mat_CompositeLink next = shell->head,oldnext;
29: while (next) {
30: MatDestroy(&next->mat);
31: if (next->work && (!next->next || next->work != next->next->work)) {
32: VecDestroy(&next->work);
33: }
34: oldnext = next;
35: next = next->next;
36: PetscFree(oldnext);
37: }
38: VecDestroy(&shell->work);
39: VecDestroy(&shell->left);
40: VecDestroy(&shell->right);
41: VecDestroy(&shell->leftwork);
42: VecDestroy(&shell->rightwork);
43: PetscFree(mat->data);
44: return(0);
45: }
49: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
50: {
51: Mat_Composite *shell = (Mat_Composite*)A->data;
52: Mat_CompositeLink next = shell->head;
53: PetscErrorCode ierr;
54: Vec in,out;
57: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
58: in = x;
59: if (shell->right) {
60: if (!shell->rightwork) {
61: VecDuplicate(shell->right,&shell->rightwork);
62: }
63: VecPointwiseMult(shell->rightwork,shell->right,in);
64: in = shell->rightwork;
65: }
66: while (next->next) {
67: if (!next->work) { /* should reuse previous work if the same size */
68: MatCreateVecs(next->mat,NULL,&next->work);
69: }
70: out = next->work;
71: MatMult(next->mat,in,out);
72: in = out;
73: next = next->next;
74: }
75: MatMult(next->mat,in,y);
76: if (shell->left) {
77: VecPointwiseMult(y,shell->left,y);
78: }
79: VecScale(y,shell->scale);
80: return(0);
81: }
85: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
86: {
87: Mat_Composite *shell = (Mat_Composite*)A->data;
88: Mat_CompositeLink tail = shell->tail;
89: PetscErrorCode ierr;
90: Vec in,out;
93: if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
94: in = x;
95: if (shell->left) {
96: if (!shell->leftwork) {
97: VecDuplicate(shell->left,&shell->leftwork);
98: }
99: VecPointwiseMult(shell->leftwork,shell->left,in);
100: in = shell->leftwork;
101: }
102: while (tail->prev) {
103: if (!tail->prev->work) { /* should reuse previous work if the same size */
104: MatCreateVecs(tail->mat,NULL,&tail->prev->work);
105: }
106: out = tail->prev->work;
107: MatMultTranspose(tail->mat,in,out);
108: in = out;
109: tail = tail->prev;
110: }
111: MatMultTranspose(tail->mat,in,y);
112: if (shell->right) {
113: VecPointwiseMult(y,shell->right,y);
114: }
115: VecScale(y,shell->scale);
116: return(0);
117: }
121: PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
122: {
123: Mat_Composite *shell = (Mat_Composite*)A->data;
124: Mat_CompositeLink next = shell->head;
125: PetscErrorCode ierr;
126: Vec in;
129: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
130: in = x;
131: if (shell->right) {
132: if (!shell->rightwork) {
133: VecDuplicate(shell->right,&shell->rightwork);
134: }
135: VecPointwiseMult(shell->rightwork,shell->right,in);
136: in = shell->rightwork;
137: }
138: MatMult(next->mat,in,y);
139: while ((next = next->next)) {
140: MatMultAdd(next->mat,in,y,y);
141: }
142: if (shell->left) {
143: VecPointwiseMult(y,shell->left,y);
144: }
145: VecScale(y,shell->scale);
146: return(0);
147: }
151: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
152: {
153: Mat_Composite *shell = (Mat_Composite*)A->data;
154: Mat_CompositeLink next = shell->head;
155: PetscErrorCode ierr;
156: Vec in;
159: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
160: in = x;
161: if (shell->left) {
162: if (!shell->leftwork) {
163: VecDuplicate(shell->left,&shell->leftwork);
164: }
165: VecPointwiseMult(shell->leftwork,shell->left,in);
166: in = shell->leftwork;
167: }
168: MatMultTranspose(next->mat,in,y);
169: while ((next = next->next)) {
170: MatMultTransposeAdd(next->mat,in,y,y);
171: }
172: if (shell->right) {
173: VecPointwiseMult(y,shell->right,y);
174: }
175: VecScale(y,shell->scale);
176: return(0);
177: }
181: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
182: {
183: Mat_Composite *shell = (Mat_Composite*)A->data;
184: Mat_CompositeLink next = shell->head;
185: PetscErrorCode ierr;
188: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
189: if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling");
191: MatGetDiagonal(next->mat,v);
192: if (next->next && !shell->work) {
193: VecDuplicate(v,&shell->work);
194: }
195: while ((next = next->next)) {
196: MatGetDiagonal(next->mat,shell->work);
197: VecAXPY(v,1.0,shell->work);
198: }
199: VecScale(v,shell->scale);
200: return(0);
201: }
205: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
206: {
208: PetscBool flg = PETSC_FALSE;
211: PetscOptionsGetBool(((PetscObject)Y)->prefix,"-mat_composite_merge",&flg,NULL);
212: if (flg) {
213: MatCompositeMerge(Y);
214: }
215: return(0);
216: }
220: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
221: {
222: Mat_Composite *a = (Mat_Composite*)inA->data;
225: a->scale *= alpha;
226: return(0);
227: }
231: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
232: {
233: Mat_Composite *a = (Mat_Composite*)inA->data;
237: if (left) {
238: if (!a->left) {
239: VecDuplicate(left,&a->left);
240: VecCopy(left,a->left);
241: } else {
242: VecPointwiseMult(a->left,left,a->left);
243: }
244: }
245: if (right) {
246: if (!a->right) {
247: VecDuplicate(right,&a->right);
248: VecCopy(right,a->right);
249: } else {
250: VecPointwiseMult(a->right,right,a->right);
251: }
252: }
253: return(0);
254: }
256: static struct _MatOps MatOps_Values = {0,
257: 0,
258: 0,
259: MatMult_Composite,
260: 0,
261: /* 5*/ MatMultTranspose_Composite,
262: 0,
263: 0,
264: 0,
265: 0,
266: /* 10*/ 0,
267: 0,
268: 0,
269: 0,
270: 0,
271: /* 15*/ 0,
272: 0,
273: MatGetDiagonal_Composite,
274: MatDiagonalScale_Composite,
275: 0,
276: /* 20*/ 0,
277: MatAssemblyEnd_Composite,
278: 0,
279: 0,
280: /* 24*/ 0,
281: 0,
282: 0,
283: 0,
284: 0,
285: /* 29*/ 0,
286: 0,
287: 0,
288: 0,
289: 0,
290: /* 34*/ 0,
291: 0,
292: 0,
293: 0,
294: 0,
295: /* 39*/ 0,
296: 0,
297: 0,
298: 0,
299: 0,
300: /* 44*/ 0,
301: MatScale_Composite,
302: MatShift_Basic,
303: 0,
304: 0,
305: /* 49*/ 0,
306: 0,
307: 0,
308: 0,
309: 0,
310: /* 54*/ 0,
311: 0,
312: 0,
313: 0,
314: 0,
315: /* 59*/ 0,
316: MatDestroy_Composite,
317: 0,
318: 0,
319: 0,
320: /* 64*/ 0,
321: 0,
322: 0,
323: 0,
324: 0,
325: /* 69*/ 0,
326: 0,
327: 0,
328: 0,
329: 0,
330: /* 74*/ 0,
331: 0,
332: 0,
333: 0,
334: 0,
335: /* 79*/ 0,
336: 0,
337: 0,
338: 0,
339: 0,
340: /* 84*/ 0,
341: 0,
342: 0,
343: 0,
344: 0,
345: /* 89*/ 0,
346: 0,
347: 0,
348: 0,
349: 0,
350: /* 94*/ 0,
351: 0,
352: 0,
353: 0,
354: 0,
355: /*99*/ 0,
356: 0,
357: 0,
358: 0,
359: 0,
360: /*104*/ 0,
361: 0,
362: 0,
363: 0,
364: 0,
365: /*109*/ 0,
366: 0,
367: 0,
368: 0,
369: 0,
370: /*114*/ 0,
371: 0,
372: 0,
373: 0,
374: 0,
375: /*119*/ 0,
376: 0,
377: 0,
378: 0,
379: 0,
380: /*124*/ 0,
381: 0,
382: 0,
383: 0,
384: 0,
385: /*129*/ 0,
386: 0,
387: 0,
388: 0,
389: 0,
390: /*134*/ 0,
391: 0,
392: 0,
393: 0,
394: 0,
395: /*139*/ 0,
396: 0,
397: 0
398: };
400: /*MC
401: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout).
403: Notes: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
405: Level: advanced
407: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
408: M*/
412: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
413: {
414: Mat_Composite *b;
418: PetscNewLog(A,&b);
419: A->data = (void*)b;
420: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
422: PetscLayoutSetUp(A->rmap);
423: PetscLayoutSetUp(A->cmap);
425: A->assembled = PETSC_TRUE;
426: A->preallocated = PETSC_TRUE;
427: b->type = MAT_COMPOSITE_ADDITIVE;
428: b->scale = 1.0;
429: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
430: return(0);
431: }
435: /*@C
436: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
438: Collective on MPI_Comm
440: Input Parameters:
441: + comm - MPI communicator
442: . nmat - number of matrices to put in
443: - mats - the matrices
445: Output Parameter:
446: . A - the matrix
448: Level: advanced
450: Notes:
451: Alternative construction
452: $ MatCreate(comm,&mat);
453: $ MatSetSizes(mat,m,n,M,N);
454: $ MatSetType(mat,MATCOMPOSITE);
455: $ MatCompositeAddMat(mat,mats[0]);
456: $ ....
457: $ MatCompositeAddMat(mat,mats[nmat-1]);
458: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
459: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
461: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
463: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
465: @*/
466: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
467: {
469: PetscInt m,n,M,N,i;
472: if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
475: MatGetLocalSize(mats[0],&m,&n);
476: MatGetSize(mats[0],&M,&N);
477: MatCreate(comm,mat);
478: MatSetSizes(*mat,m,n,M,N);
479: MatSetType(*mat,MATCOMPOSITE);
480: for (i=0; i<nmat; i++) {
481: MatCompositeAddMat(*mat,mats[i]);
482: }
483: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
484: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
485: return(0);
486: }
490: /*@
491: MatCompositeAddMat - add another matrix to a composite matrix
493: Collective on Mat
495: Input Parameters:
496: + mat - the composite matrix
497: - smat - the partial matrix
499: Level: advanced
501: .seealso: MatCreateComposite()
502: @*/
503: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
504: {
505: Mat_Composite *shell;
506: PetscErrorCode ierr;
507: Mat_CompositeLink ilink,next;
512: PetscNewLog(mat,&ilink);
513: ilink->next = 0;
514: PetscObjectReference((PetscObject)smat);
515: ilink->mat = smat;
517: shell = (Mat_Composite*)mat->data;
518: next = shell->head;
519: if (!next) shell->head = ilink;
520: else {
521: while (next->next) {
522: next = next->next;
523: }
524: next->next = ilink;
525: ilink->prev = next;
526: }
527: shell->tail = ilink;
528: return(0);
529: }
533: /*@C
534: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
536: Collective on MPI_Comm
538: Input Parameters:
539: . mat - the composite matrix
542: Level: advanced
544: Notes:
545: The MatType of the resulting matrix will be the same as the MatType of the FIRST
546: matrix in the composite matrix.
548: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
550: @*/
551: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
552: {
553: Mat_Composite *b = (Mat_Composite*)mat->data;
554: PetscBool flg;
558: PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
559: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
560: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
561: mat->ops->getdiagonal = 0;
562: mat->ops->mult = MatMult_Composite_Multiplicative;
563: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
564: b->type = MAT_COMPOSITE_MULTIPLICATIVE;
565: } else {
566: mat->ops->getdiagonal = MatGetDiagonal_Composite;
567: mat->ops->mult = MatMult_Composite;
568: mat->ops->multtranspose = MatMultTranspose_Composite;
569: b->type = MAT_COMPOSITE_ADDITIVE;
570: }
571: return(0);
572: }
577: /*@C
578: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
579: by summing all the matrices inside the composite matrix.
581: Collective on MPI_Comm
583: Input Parameters:
584: . mat - the composite matrix
587: Options Database:
588: . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
590: Level: advanced
592: Notes:
593: The MatType of the resulting matrix will be the same as the MatType of the FIRST
594: matrix in the composite matrix.
596: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
598: @*/
599: PetscErrorCode MatCompositeMerge(Mat mat)
600: {
601: Mat_Composite *shell = (Mat_Composite*)mat->data;
602: Mat_CompositeLink next = shell->head, prev = shell->tail;
603: PetscErrorCode ierr;
604: Mat tmat,newmat;
605: Vec left,right;
606: PetscScalar scale;
609: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
612: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
613: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
614: while ((next = next->next)) {
615: MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
616: }
617: } else {
618: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
619: while ((prev = prev->prev)) {
620: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
621: MatDestroy(&tmat);
622: tmat = newmat;
623: }
624: }
626: scale = shell->scale;
627: if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
628: if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}
630: MatHeaderReplace(mat,tmat);
632: MatDiagonalScale(mat,left,right);
633: MatScale(mat,scale);
634: VecDestroy(&left);
635: VecDestroy(&right);
636: return(0);
637: }