Actual source code: mcomposite.c
petsc-3.4.5 2014-06-29
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: MatGetVecs(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: MatGetVecs(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: 0,
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: };
399: /*MC
400: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout).
402: Notes: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
404: Level: advanced
406: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
407: M*/
411: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
412: {
413: Mat_Composite *b;
417: PetscNewLog(A,Mat_Composite,&b);
418: A->data = (void*)b;
419: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
421: PetscLayoutSetUp(A->rmap);
422: PetscLayoutSetUp(A->cmap);
424: A->assembled = PETSC_TRUE;
425: A->preallocated = PETSC_TRUE;
426: b->type = MAT_COMPOSITE_ADDITIVE;
427: b->scale = 1.0;
428: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
429: return(0);
430: }
434: /*@C
435: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
437: Collective on MPI_Comm
439: Input Parameters:
440: + comm - MPI communicator
441: . nmat - number of matrices to put in
442: - mats - the matrices
444: Output Parameter:
445: . A - the matrix
447: Level: advanced
449: Notes:
450: Alternative construction
451: $ MatCreate(comm,&mat);
452: $ MatSetSizes(mat,m,n,M,N);
453: $ MatSetType(mat,MATCOMPOSITE);
454: $ MatCompositeAddMat(mat,mats[0]);
455: $ ....
456: $ MatCompositeAddMat(mat,mats[nmat-1]);
457: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
458: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
460: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
462: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
464: @*/
465: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
466: {
468: PetscInt m,n,M,N,i;
471: if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
474: MatGetLocalSize(mats[0],&m,&n);
475: MatGetSize(mats[0],&M,&N);
476: MatCreate(comm,mat);
477: MatSetSizes(*mat,m,n,M,N);
478: MatSetType(*mat,MATCOMPOSITE);
479: for (i=0; i<nmat; i++) {
480: MatCompositeAddMat(*mat,mats[i]);
481: }
482: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
483: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
484: return(0);
485: }
489: /*@
490: MatCompositeAddMat - add another matrix to a composite matrix
492: Collective on Mat
494: Input Parameters:
495: + mat - the composite matrix
496: - smat - the partial matrix
498: Level: advanced
500: .seealso: MatCreateComposite()
501: @*/
502: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
503: {
504: Mat_Composite *shell;
505: PetscErrorCode ierr;
506: Mat_CompositeLink ilink,next;
511: PetscNewLog(mat,struct _Mat_CompositeLink,&ilink);
512: ilink->next = 0;
513: PetscObjectReference((PetscObject)smat);
514: ilink->mat = smat;
516: shell = (Mat_Composite*)mat->data;
517: next = shell->head;
518: if (!next) shell->head = ilink;
519: else {
520: while (next->next) {
521: next = next->next;
522: }
523: next->next = ilink;
524: ilink->prev = next;
525: }
526: shell->tail = ilink;
527: return(0);
528: }
532: /*@C
533: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
535: Collective on MPI_Comm
537: Input Parameters:
538: . mat - the composite matrix
541: Level: advanced
543: Notes:
544: The MatType of the resulting matrix will be the same as the MatType of the FIRST
545: matrix in the composite matrix.
547: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
549: @*/
550: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
551: {
552: Mat_Composite *b = (Mat_Composite*)mat->data;
553: PetscBool flg;
557: PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
558: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
559: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
560: mat->ops->getdiagonal = 0;
561: mat->ops->mult = MatMult_Composite_Multiplicative;
562: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
563: b->type = MAT_COMPOSITE_MULTIPLICATIVE;
564: } else {
565: mat->ops->getdiagonal = MatGetDiagonal_Composite;
566: mat->ops->mult = MatMult_Composite;
567: mat->ops->multtranspose = MatMultTranspose_Composite;
568: b->type = MAT_COMPOSITE_ADDITIVE;
569: }
570: return(0);
571: }
576: /*@C
577: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
578: by summing all the matrices inside the composite matrix.
580: Collective on MPI_Comm
582: Input Parameters:
583: . mat - the composite matrix
586: Options Database:
587: . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
589: Level: advanced
591: Notes:
592: The MatType of the resulting matrix will be the same as the MatType of the FIRST
593: matrix in the composite matrix.
595: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
597: @*/
598: PetscErrorCode MatCompositeMerge(Mat mat)
599: {
600: Mat_Composite *shell = (Mat_Composite*)mat->data;
601: Mat_CompositeLink next = shell->head, prev = shell->tail;
602: PetscErrorCode ierr;
603: Mat tmat,newmat;
604: Vec left,right;
605: PetscScalar scale;
608: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
611: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
612: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
613: while ((next = next->next)) {
614: MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
615: }
616: } else {
617: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
618: while ((prev = prev->prev)) {
619: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
620: MatDestroy(&tmat);
621: tmat = newmat;
622: }
623: }
625: scale = shell->scale;
626: if ((left = shell->left)) {PetscObjectReference((PetscObject)left);}
627: if ((right = shell->right)) {PetscObjectReference((PetscObject)right);}
629: MatHeaderReplace(mat,tmat);
631: MatDiagonalScale(mat,left,right);
632: MatScale(mat,scale);
633: VecDestroy(&left);
634: VecDestroy(&right);
635: return(0);
636: }