Actual source code: mcomposite.c
petsc-3.3-p7 2013-05-11
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: };
10:
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,PETSC_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,PETSC_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,PETSC_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};
355: /*MC
356: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices (all matrices are of same size and parallel layout).
358: Notes: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
360: Level: advanced
362: .seealso: MatCreateComposite(), MatCompositeAddMat(), MatSetType(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
363: M*/
365: EXTERN_C_BEGIN
368: PetscErrorCode MatCreate_Composite(Mat A)
369: {
370: Mat_Composite *b;
374: PetscNewLog(A,Mat_Composite,&b);
375: A->data = (void*)b;
376: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
378: PetscLayoutSetUp(A->rmap);
379: PetscLayoutSetUp(A->cmap);
381: A->assembled = PETSC_TRUE;
382: A->preallocated = PETSC_TRUE;
383: b->type = MAT_COMPOSITE_ADDITIVE;
384: b->scale = 1.0;
385: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
386: return(0);
387: }
388: EXTERN_C_END
392: /*@C
393: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
395: Collective on MPI_Comm
397: Input Parameters:
398: + comm - MPI communicator
399: . nmat - number of matrices to put in
400: - mats - the matrices
402: Output Parameter:
403: . A - the matrix
405: Level: advanced
407: Notes:
408: Alternative construction
409: $ MatCreate(comm,&mat);
410: $ MatSetSizes(mat,m,n,M,N);
411: $ MatSetType(mat,MATCOMPOSITE);
412: $ MatCompositeAddMat(mat,mats[0]);
413: $ ....
414: $ MatCompositeAddMat(mat,mats[nmat-1]);
415: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
416: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
418: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
420: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeMerge(), MatCompositeSetType(), MatCompositeType
422: @*/
423: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
424: {
426: PetscInt m,n,M,N,i;
427:
429: if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
432: MatGetLocalSize(mats[0],&m,&n);
433: MatGetSize(mats[0],&M,&N);
434: MatCreate(comm,mat);
435: MatSetSizes(*mat,m,n,M,N);
436: MatSetType(*mat,MATCOMPOSITE);
437: for (i=0; i<nmat; i++) {
438: MatCompositeAddMat(*mat,mats[i]);
439: }
440: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
441: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
442: return(0);
443: }
447: /*@
448: MatCompositeAddMat - add another matrix to a composite matrix
450: Collective on Mat
452: Input Parameters:
453: + mat - the composite matrix
454: - smat - the partial matrix
456: Level: advanced
458: .seealso: MatCreateComposite()
459: @*/
460: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
461: {
462: Mat_Composite *shell;
463: PetscErrorCode ierr;
464: Mat_CompositeLink ilink,next;
469: PetscNewLog(mat,struct _Mat_CompositeLink,&ilink);
470: ilink->next = 0;
471: PetscObjectReference((PetscObject)smat);
472: ilink->mat = smat;
474: shell = (Mat_Composite*)mat->data;
475: next = shell->head;
476: if (!next) {
477: shell->head = ilink;
478: } else {
479: while (next->next) {
480: next = next->next;
481: }
482: next->next = ilink;
483: ilink->prev = next;
484: }
485: shell->tail = ilink;
486: return(0);
487: }
491: /*@C
492: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product
494: Collective on MPI_Comm
496: Input Parameters:
497: . mat - the composite matrix
500: Level: advanced
502: Notes:
503: The MatType of the resulting matrix will be the same as the MatType of the FIRST
504: matrix in the composite matrix.
506: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
508: @*/
509: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
510: {
511: Mat_Composite *b = (Mat_Composite*)mat->data;
512: PetscBool flg;
516: PetscObjectTypeCompare((PetscObject)mat,MATCOMPOSITE,&flg);
517: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Can only use with composite matrix");
518: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
519: mat->ops->getdiagonal = 0;
520: mat->ops->mult = MatMult_Composite_Multiplicative;
521: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
522: b->type = MAT_COMPOSITE_MULTIPLICATIVE;
523: } else {
524: mat->ops->getdiagonal = MatGetDiagonal_Composite;
525: mat->ops->mult = MatMult_Composite;
526: mat->ops->multtranspose = MatMultTranspose_Composite;
527: b->type = MAT_COMPOSITE_ADDITIVE;
528: }
529: return(0);
530: }
535: /*@C
536: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
537: by summing all the matrices inside the composite matrix.
539: Collective on MPI_Comm
541: Input Parameters:
542: . mat - the composite matrix
545: Options Database:
546: . -mat_composite_merge (you must call MatAssemblyBegin()/MatAssemblyEnd() to have this checked)
548: Level: advanced
550: Notes:
551: The MatType of the resulting matrix will be the same as the MatType of the FIRST
552: matrix in the composite matrix.
554: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MATCOMPOSITE
556: @*/
557: PetscErrorCode MatCompositeMerge(Mat mat)
558: {
559: Mat_Composite *shell = (Mat_Composite*)mat->data;
560: Mat_CompositeLink next = shell->head, prev = shell->tail;
561: PetscErrorCode ierr;
562: Mat tmat,newmat;
565: if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
568: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
569: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
570: while ((next = next->next)) {
571: MatAXPY(tmat,1.0,next->mat,DIFFERENT_NONZERO_PATTERN);
572: }
573: } else {
574: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
575: while ((prev = prev->prev)) {
576: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
577: MatDestroy(&tmat);
578: tmat = newmat;
579: }
580: }
581: MatHeaderReplace(mat,tmat);
582: MatDiagonalScale(mat,shell->left,shell->right);
583: MatScale(mat,shell->scale);
584: return(0);
585: }