Actual source code: mcomposite.c
2: #include <petsc/private/matimpl.h>
4: const char *const MatCompositeMergeTypes[] = {"left","right","MatCompositeMergeType","MAT_COMPOSITE_",NULL};
6: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
7: struct _Mat_CompositeLink {
8: Mat mat;
9: Vec work;
10: Mat_CompositeLink next,prev;
11: };
13: typedef struct {
14: MatCompositeType type;
15: Mat_CompositeLink head,tail;
16: Vec work;
17: PetscScalar scale; /* scale factor supplied with MatScale() */
18: Vec left,right; /* left and right diagonal scaling provided with MatDiagonalScale() */
19: Vec leftwork,rightwork,leftwork2,rightwork2; /* Two pairs of working vectors */
20: PetscInt nmat;
21: PetscBool merge;
22: MatCompositeMergeType mergetype;
23: MatStructure structure;
25: PetscScalar *scalings;
26: PetscBool merge_mvctx; /* Whether need to merge mvctx of component matrices */
27: Vec *lvecs; /* [nmat] Basically, they are Mvctx->lvec of each component matrix */
28: PetscScalar *larray; /* [len] Data arrays of lvecs[] are stored consecutively in larray */
29: PetscInt len; /* Length of larray[] */
30: Vec gvec; /* Union of lvecs[] without duplicated entries */
31: PetscInt *location; /* A map that maps entries in garray[] to larray[] */
32: VecScatter Mvctx;
33: } Mat_Composite;
35: PetscErrorCode MatDestroy_Composite(Mat mat)
36: {
37: Mat_Composite *shell = (Mat_Composite*)mat->data;
38: Mat_CompositeLink next = shell->head,oldnext;
39: PetscInt i;
41: while (next) {
42: MatDestroy(&next->mat);
43: if (next->work && (!next->next || next->work != next->next->work)) {
44: VecDestroy(&next->work);
45: }
46: oldnext = next;
47: next = next->next;
48: PetscFree(oldnext);
49: }
50: VecDestroy(&shell->work);
51: VecDestroy(&shell->left);
52: VecDestroy(&shell->right);
53: VecDestroy(&shell->leftwork);
54: VecDestroy(&shell->rightwork);
55: VecDestroy(&shell->leftwork2);
56: VecDestroy(&shell->rightwork2);
58: if (shell->Mvctx) {
59: for (i=0; i<shell->nmat; i++) VecDestroy(&shell->lvecs[i]);
60: PetscFree3(shell->location,shell->larray,shell->lvecs);
61: PetscFree(shell->larray);
62: VecDestroy(&shell->gvec);
63: VecScatterDestroy(&shell->Mvctx);
64: }
66: PetscFree(shell->scalings);
67: PetscFree(mat->data);
68: return 0;
69: }
71: PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y)
72: {
73: Mat_Composite *shell = (Mat_Composite*)A->data;
74: Mat_CompositeLink next = shell->head;
75: Vec in,out;
76: PetscScalar scale;
77: PetscInt i;
80: in = x;
81: if (shell->right) {
82: if (!shell->rightwork) {
83: VecDuplicate(shell->right,&shell->rightwork);
84: }
85: VecPointwiseMult(shell->rightwork,shell->right,in);
86: in = shell->rightwork;
87: }
88: while (next->next) {
89: if (!next->work) { /* should reuse previous work if the same size */
90: MatCreateVecs(next->mat,NULL,&next->work);
91: }
92: out = next->work;
93: MatMult(next->mat,in,out);
94: in = out;
95: next = next->next;
96: }
97: MatMult(next->mat,in,y);
98: if (shell->left) {
99: VecPointwiseMult(y,shell->left,y);
100: }
101: scale = shell->scale;
102: if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];}
103: VecScale(y,scale);
104: return 0;
105: }
107: PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y)
108: {
109: Mat_Composite *shell = (Mat_Composite*)A->data;
110: Mat_CompositeLink tail = shell->tail;
111: Vec in,out;
112: PetscScalar scale;
113: PetscInt i;
116: in = x;
117: if (shell->left) {
118: if (!shell->leftwork) {
119: VecDuplicate(shell->left,&shell->leftwork);
120: }
121: VecPointwiseMult(shell->leftwork,shell->left,in);
122: in = shell->leftwork;
123: }
124: while (tail->prev) {
125: if (!tail->prev->work) { /* should reuse previous work if the same size */
126: MatCreateVecs(tail->mat,NULL,&tail->prev->work);
127: }
128: out = tail->prev->work;
129: MatMultTranspose(tail->mat,in,out);
130: in = out;
131: tail = tail->prev;
132: }
133: MatMultTranspose(tail->mat,in,y);
134: if (shell->right) {
135: VecPointwiseMult(y,shell->right,y);
136: }
138: scale = shell->scale;
139: if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];}
140: VecScale(y,scale);
141: return 0;
142: }
144: PetscErrorCode MatMult_Composite(Mat mat,Vec x,Vec y)
145: {
146: Mat_Composite *shell = (Mat_Composite*)mat->data;
147: Mat_CompositeLink cur = shell->head;
148: Vec in,y2,xin;
149: Mat A,B;
150: PetscInt i,j,k,n,nuniq,lo,hi,mid,*gindices,*buf,*tmp,tot;
151: const PetscScalar *vals;
152: const PetscInt *garray;
153: IS ix,iy;
154: PetscBool match;
157: in = x;
158: if (shell->right) {
159: if (!shell->rightwork) {
160: VecDuplicate(shell->right,&shell->rightwork);
161: }
162: VecPointwiseMult(shell->rightwork,shell->right,in);
163: in = shell->rightwork;
164: }
166: /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time
167: we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and
168: it is legal to merge Mvctx, because all component matrices have the same size.
169: */
170: if (shell->merge_mvctx && !shell->Mvctx) {
171: /* Currently only implemented for MATMPIAIJ */
172: for (cur=shell->head; cur; cur=cur->next) {
173: PetscObjectTypeCompare((PetscObject)cur->mat,MATMPIAIJ,&match);
174: if (!match) {
175: shell->merge_mvctx = PETSC_FALSE;
176: goto skip_merge_mvctx;
177: }
178: }
180: /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */
181: tot = 0;
182: for (cur=shell->head; cur; cur=cur->next) {
183: MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,NULL);
184: MatGetLocalSize(B,NULL,&n);
185: tot += n;
186: }
187: PetscMalloc3(tot,&shell->location,tot,&shell->larray,shell->nmat,&shell->lvecs);
188: shell->len = tot;
190: /* Go through matrices second time to sort off-diag columns and remove dups */
191: PetscMalloc1(tot,&gindices); /* No Malloc2() since we will give one to petsc and free the other */
192: PetscMalloc1(tot,&buf);
193: nuniq = 0; /* Number of unique nonzero columns */
194: for (cur=shell->head; cur; cur=cur->next) {
195: MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,&garray);
196: MatGetLocalSize(B,NULL,&n);
197: /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */
198: i = j = k = 0;
199: while (i < n && j < nuniq) {
200: if (garray[i] < gindices[j]) buf[k++] = garray[i++];
201: else if (garray[i] > gindices[j]) buf[k++] = gindices[j++];
202: else {buf[k++] = garray[i++]; j++;}
203: }
204: /* Copy leftover in garray[] or gindices[] */
205: if (i < n) {
206: PetscArraycpy(buf+k,garray+i,n-i);
207: nuniq = k + n-i;
208: } else if (j < nuniq) {
209: PetscArraycpy(buf+k,gindices+j,nuniq-j);
210: nuniq = k + nuniq-j;
211: } else nuniq = k;
212: /* Swap gindices and buf to merge garray of the next matrix */
213: tmp = gindices;
214: gindices = buf;
215: buf = tmp;
216: }
217: PetscFree(buf);
219: /* Go through matrices third time to build a map from gindices[] to garray[] */
220: tot = 0;
221: for (cur=shell->head,j=0; cur; cur=cur->next,j++) { /* j-th matrix */
222: MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,&garray);
223: MatGetLocalSize(B,NULL,&n);
224: VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,NULL,&shell->lvecs[j]);
225: /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */
226: lo = 0;
227: for (i=0; i<n; i++) {
228: hi = nuniq;
229: while (hi - lo > 1) {
230: mid = lo + (hi - lo)/2;
231: if (garray[i] < gindices[mid]) hi = mid;
232: else lo = mid;
233: }
234: shell->location[tot+i] = lo; /* gindices[lo] = garray[i] */
235: lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */
236: }
237: tot += n;
238: }
240: /* Build merged Mvctx */
241: ISCreateGeneral(PETSC_COMM_SELF,nuniq,gindices,PETSC_OWN_POINTER,&ix);
242: ISCreateStride(PETSC_COMM_SELF,nuniq,0,1,&iy);
243: VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mat->cmap->n,mat->cmap->N,NULL,&xin);
244: VecCreateSeq(PETSC_COMM_SELF,nuniq,&shell->gvec);
245: VecScatterCreate(xin,ix,shell->gvec,iy,&shell->Mvctx);
246: VecDestroy(&xin);
247: ISDestroy(&ix);
248: ISDestroy(&iy);
249: }
251: skip_merge_mvctx:
252: VecSet(y,0);
253: if (!shell->leftwork2) VecDuplicate(y,&shell->leftwork2);
254: y2 = shell->leftwork2;
256: if (shell->Mvctx) { /* Have a merged Mvctx */
257: /* Suppose we want to compute y = sMx, where s is the scaling factor and A, B are matrix M's diagonal/off-diagonal part. We could do
258: in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an oppertunity to
259: overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach.
260: */
261: VecScatterBegin(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);
262: VecScatterEnd(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);
264: VecGetArrayRead(shell->gvec,&vals);
265: for (i=0; i<shell->len; i++) shell->larray[i] = vals[shell->location[i]];
266: VecRestoreArrayRead(shell->gvec,&vals);
268: for (cur=shell->head,tot=i=0; cur; cur=cur->next,i++) { /* i-th matrix */
269: MatMPIAIJGetSeqAIJ(cur->mat,&A,&B,NULL);
270: (*A->ops->mult)(A,in,y2);
271: MatGetLocalSize(B,NULL,&n);
272: VecPlaceArray(shell->lvecs[i],&shell->larray[tot]);
273: (*B->ops->multadd)(B,shell->lvecs[i],y2,y2);
274: VecResetArray(shell->lvecs[i]);
275: VecAXPY(y,(shell->scalings ? shell->scalings[i] : 1.0),y2);
276: tot += n;
277: }
278: } else {
279: if (shell->scalings) {
280: for (cur=shell->head,i=0; cur; cur=cur->next,i++) {
281: MatMult(cur->mat,in,y2);
282: VecAXPY(y,shell->scalings[i],y2);
283: }
284: } else {
285: for (cur=shell->head; cur; cur=cur->next) MatMultAdd(cur->mat,in,y,y);
286: }
287: }
289: if (shell->left) VecPointwiseMult(y,shell->left,y);
290: VecScale(y,shell->scale);
291: return 0;
292: }
294: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
295: {
296: Mat_Composite *shell = (Mat_Composite*)A->data;
297: Mat_CompositeLink next = shell->head;
298: Vec in,y2 = NULL;
299: PetscInt i;
302: in = x;
303: if (shell->left) {
304: if (!shell->leftwork) {
305: VecDuplicate(shell->left,&shell->leftwork);
306: }
307: VecPointwiseMult(shell->leftwork,shell->left,in);
308: in = shell->leftwork;
309: }
311: MatMultTranspose(next->mat,in,y);
312: if (shell->scalings) {
313: VecScale(y,shell->scalings[0]);
314: if (!shell->rightwork2) VecDuplicate(y,&shell->rightwork2);
315: y2 = shell->rightwork2;
316: }
317: i = 1;
318: while ((next = next->next)) {
319: if (!shell->scalings) MatMultTransposeAdd(next->mat,in,y,y);
320: else {
321: MatMultTranspose(next->mat,in,y2);
322: VecAXPY(y,shell->scalings[i++],y2);
323: }
324: }
325: if (shell->right) {
326: VecPointwiseMult(y,shell->right,y);
327: }
328: VecScale(y,shell->scale);
329: return 0;
330: }
332: PetscErrorCode MatMultAdd_Composite(Mat A,Vec x,Vec y,Vec z)
333: {
334: Mat_Composite *shell = (Mat_Composite*)A->data;
336: if (y != z) {
337: MatMult(A,x,z);
338: VecAXPY(z,1.0,y);
339: } else {
340: if (!shell->leftwork) {
341: VecDuplicate(z,&shell->leftwork);
342: }
343: MatMult(A,x,shell->leftwork);
344: VecCopy(y,z);
345: VecAXPY(z,1.0,shell->leftwork);
346: }
347: return 0;
348: }
350: PetscErrorCode MatMultTransposeAdd_Composite(Mat A,Vec x,Vec y, Vec z)
351: {
352: Mat_Composite *shell = (Mat_Composite*)A->data;
354: if (y != z) {
355: MatMultTranspose(A,x,z);
356: VecAXPY(z,1.0,y);
357: } else {
358: if (!shell->rightwork) {
359: VecDuplicate(z,&shell->rightwork);
360: }
361: MatMultTranspose(A,x,shell->rightwork);
362: VecCopy(y,z);
363: VecAXPY(z,1.0,shell->rightwork);
364: }
365: return 0;
366: }
368: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
369: {
370: Mat_Composite *shell = (Mat_Composite*)A->data;
371: Mat_CompositeLink next = shell->head;
372: PetscInt i;
377: MatGetDiagonal(next->mat,v);
378: if (shell->scalings) VecScale(v,shell->scalings[0]);
380: if (next->next && !shell->work) {
381: VecDuplicate(v,&shell->work);
382: }
383: i = 1;
384: while ((next = next->next)) {
385: MatGetDiagonal(next->mat,shell->work);
386: VecAXPY(v,(shell->scalings ? shell->scalings[i++] : 1.0),shell->work);
387: }
388: VecScale(v,shell->scale);
389: return 0;
390: }
392: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
393: {
394: Mat_Composite *shell = (Mat_Composite*)Y->data;
396: if (shell->merge) {
397: MatCompositeMerge(Y);
398: }
399: return 0;
400: }
402: PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha)
403: {
404: Mat_Composite *a = (Mat_Composite*)inA->data;
406: a->scale *= alpha;
407: return 0;
408: }
410: PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right)
411: {
412: Mat_Composite *a = (Mat_Composite*)inA->data;
414: if (left) {
415: if (!a->left) {
416: VecDuplicate(left,&a->left);
417: VecCopy(left,a->left);
418: } else {
419: VecPointwiseMult(a->left,left,a->left);
420: }
421: }
422: if (right) {
423: if (!a->right) {
424: VecDuplicate(right,&a->right);
425: VecCopy(right,a->right);
426: } else {
427: VecPointwiseMult(a->right,right,a->right);
428: }
429: }
430: return 0;
431: }
433: PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A)
434: {
435: Mat_Composite *a = (Mat_Composite*)A->data;
437: PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");
438: PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);
439: PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);
440: PetscOptionsBool("-mat_composite_merge_mvctx","Merge MatMult() vecscat contexts","MatCreateComposite",a->merge_mvctx,&a->merge_mvctx,NULL);
441: PetscOptionsTail();
442: return 0;
443: }
445: /*@
446: MatCreateComposite - Creates a matrix as the sum or product of one or more matrices
448: Collective
450: Input Parameters:
451: + comm - MPI communicator
452: . nmat - number of matrices to put in
453: - mats - the matrices
455: Output Parameter:
456: . mat - the matrix
458: Options Database Keys:
459: + -mat_composite_merge - merge in MatAssemblyEnd()
460: . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in MatMult() for ADDITIVE matrices
461: - -mat_composite_merge_type - set merge direction
463: Level: advanced
465: Notes:
466: Alternative construction
467: $ MatCreate(comm,&mat);
468: $ MatSetSizes(mat,m,n,M,N);
469: $ MatSetType(mat,MATCOMPOSITE);
470: $ MatCompositeAddMat(mat,mats[0]);
471: $ ....
472: $ MatCompositeAddMat(mat,mats[nmat-1]);
473: $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);
474: $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);
476: For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0]
478: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeGetMat(), MatCompositeMerge(), MatCompositeSetType(), MATCOMPOSITE
480: @*/
481: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
482: {
483: PetscInt m,n,M,N,i;
488: MatGetLocalSize(mats[0],PETSC_IGNORE,&n);
489: MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);
490: MatGetSize(mats[0],PETSC_IGNORE,&N);
491: MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);
492: MatCreate(comm,mat);
493: MatSetSizes(*mat,m,n,M,N);
494: MatSetType(*mat,MATCOMPOSITE);
495: for (i=0; i<nmat; i++) {
496: MatCompositeAddMat(*mat,mats[i]);
497: }
498: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
499: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
500: return 0;
501: }
503: static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat)
504: {
505: Mat_Composite *shell = (Mat_Composite*)mat->data;
506: Mat_CompositeLink ilink,next = shell->head;
508: PetscNewLog(mat,&ilink);
509: ilink->next = NULL;
510: PetscObjectReference((PetscObject)smat);
511: ilink->mat = smat;
513: if (!next) shell->head = ilink;
514: else {
515: while (next->next) {
516: next = next->next;
517: }
518: next->next = ilink;
519: ilink->prev = next;
520: }
521: shell->tail = ilink;
522: shell->nmat += 1;
524: /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */
525: if (shell->scalings) {
526: PetscRealloc(sizeof(PetscScalar)*shell->nmat,&shell->scalings);
527: shell->scalings[shell->nmat-1] = 1.0;
528: }
529: return 0;
530: }
532: /*@
533: MatCompositeAddMat - Add another matrix to a composite matrix.
535: Collective on Mat
537: Input Parameters:
538: + mat - the composite matrix
539: - smat - the partial matrix
541: Level: advanced
543: .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE
544: @*/
545: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
546: {
549: PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));
550: return 0;
551: }
553: static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type)
554: {
555: Mat_Composite *b = (Mat_Composite*)mat->data;
557: b->type = type;
558: if (type == MAT_COMPOSITE_MULTIPLICATIVE) {
559: mat->ops->getdiagonal = NULL;
560: mat->ops->mult = MatMult_Composite_Multiplicative;
561: mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative;
562: b->merge_mvctx = PETSC_FALSE;
563: } else {
564: mat->ops->getdiagonal = MatGetDiagonal_Composite;
565: mat->ops->mult = MatMult_Composite;
566: mat->ops->multtranspose = MatMultTranspose_Composite;
567: }
568: return 0;
569: }
571: /*@
572: MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product.
574: Logically Collective on Mat
576: Input Parameters:
577: . mat - the composite matrix
579: Level: advanced
581: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE
583: @*/
584: PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type)
585: {
588: PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));
589: return 0;
590: }
592: static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type)
593: {
594: Mat_Composite *b = (Mat_Composite*)mat->data;
596: *type = b->type;
597: return 0;
598: }
600: /*@
601: MatCompositeGetType - Returns type of composite.
603: Not Collective
605: Input Parameter:
606: . mat - the composite matrix
608: Output Parameter:
609: . type - type of composite
611: Level: advanced
613: .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE
615: @*/
616: PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type)
617: {
620: PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));
621: return 0;
622: }
624: static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str)
625: {
626: Mat_Composite *b = (Mat_Composite*)mat->data;
628: b->structure = str;
629: return 0;
630: }
632: /*@
633: MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix.
635: Not Collective
637: Input Parameters:
638: + mat - the composite matrix
639: - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN
641: Level: advanced
643: Notes:
644: Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix.
646: .seealso: MatAXPY(), MatCreateComposite(), MatCompositeMerge() MatCompositeGetMatStructure(), MATCOMPOSITE
648: @*/
649: PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str)
650: {
652: PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));
653: return 0;
654: }
656: static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str)
657: {
658: Mat_Composite *b = (Mat_Composite*)mat->data;
660: *str = b->structure;
661: return 0;
662: }
664: /*@
665: MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix.
667: Not Collective
669: Input Parameter:
670: . mat - the composite matrix
672: Output Parameter:
673: . str - structure of the matrices
675: Level: advanced
677: .seealso: MatCreateComposite(), MatCompositeSetMatStructure(), MATCOMPOSITE
679: @*/
680: PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str)
681: {
684: PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));
685: return 0;
686: }
688: static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type)
689: {
690: Mat_Composite *shell = (Mat_Composite*)mat->data;
692: shell->mergetype = type;
693: return 0;
694: }
696: /*@
697: MatCompositeSetMergeType - Sets order of MatCompositeMerge().
699: Logically Collective on Mat
701: Input Parameters:
702: + mat - the composite matrix
703: - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]),
704: MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1])
706: Level: advanced
708: Notes:
709: The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed.
710: If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0])))
711: otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0].
713: .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE
715: @*/
716: PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type)
717: {
720: PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));
721: return 0;
722: }
724: static PetscErrorCode MatCompositeMerge_Composite(Mat mat)
725: {
726: Mat_Composite *shell = (Mat_Composite*)mat->data;
727: Mat_CompositeLink next = shell->head, prev = shell->tail;
728: Mat tmat,newmat;
729: Vec left,right;
730: PetscScalar scale;
731: PetscInt i;
734: scale = shell->scale;
735: if (shell->type == MAT_COMPOSITE_ADDITIVE) {
736: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
737: i = 0;
738: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
739: if (shell->scalings) MatScale(tmat,shell->scalings[i++]);
740: while ((next = next->next)) {
741: MatAXPY(tmat,(shell->scalings ? shell->scalings[i++] : 1.0),next->mat,shell->structure);
742: }
743: } else {
744: i = shell->nmat-1;
745: MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
746: if (shell->scalings) MatScale(tmat,shell->scalings[i--]);
747: while ((prev = prev->prev)) {
748: MatAXPY(tmat,(shell->scalings ? shell->scalings[i--] : 1.0),prev->mat,shell->structure);
749: }
750: }
751: } else {
752: if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) {
753: MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);
754: while ((next = next->next)) {
755: MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
756: MatDestroy(&tmat);
757: tmat = newmat;
758: }
759: } else {
760: MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);
761: while ((prev = prev->prev)) {
762: MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);
763: MatDestroy(&tmat);
764: tmat = newmat;
765: }
766: }
767: if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];}
768: }
770: if ((left = shell->left)) PetscObjectReference((PetscObject)left);
771: if ((right = shell->right)) PetscObjectReference((PetscObject)right);
773: MatHeaderReplace(mat,&tmat);
775: MatDiagonalScale(mat,left,right);
776: MatScale(mat,scale);
777: VecDestroy(&left);
778: VecDestroy(&right);
779: return 0;
780: }
782: /*@
783: MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix
784: by summing or computing the product of all the matrices inside the composite matrix.
786: Collective
788: Input Parameter:
789: . mat - the composite matrix
791: Options Database Keys:
792: + -mat_composite_merge - merge in MatAssemblyEnd()
793: - -mat_composite_merge_type - set merge direction
795: Level: advanced
797: Notes:
798: The MatType of the resulting matrix will be the same as the MatType of the FIRST
799: matrix in the composite matrix.
801: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMatStructure(), MatCompositeSetMergeType(), MATCOMPOSITE
803: @*/
804: PetscErrorCode MatCompositeMerge(Mat mat)
805: {
807: PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));
808: return 0;
809: }
811: static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat)
812: {
813: Mat_Composite *shell = (Mat_Composite*)mat->data;
815: *nmat = shell->nmat;
816: return 0;
817: }
819: /*@
820: MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix.
822: Not Collective
824: Input Parameter:
825: . mat - the composite matrix
827: Output Parameter:
828: . nmat - number of matrices in the composite matrix
830: Level: advanced
832: .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE
834: @*/
835: PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat)
836: {
839: PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));
840: return 0;
841: }
843: static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai)
844: {
845: Mat_Composite *shell = (Mat_Composite*)mat->data;
846: Mat_CompositeLink ilink;
847: PetscInt k;
850: ilink = shell->head;
851: for (k=0; k<i; k++) {
852: ilink = ilink->next;
853: }
854: *Ai = ilink->mat;
855: return 0;
856: }
858: /*@
859: MatCompositeGetMat - Returns the ith matrix from the composite matrix.
861: Logically Collective on Mat
863: Input Parameters:
864: + mat - the composite matrix
865: - i - the number of requested matrix
867: Output Parameter:
868: . Ai - ith matrix in composite
870: Level: advanced
872: .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE
874: @*/
875: PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai)
876: {
880: PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));
881: return 0;
882: }
884: PetscErrorCode MatCompositeSetScalings_Composite(Mat mat,const PetscScalar *scalings)
885: {
886: Mat_Composite *shell = (Mat_Composite*)mat->data;
887: PetscInt nmat;
889: MatCompositeGetNumberMat(mat,&nmat);
890: if (!shell->scalings) PetscMalloc1(nmat,&shell->scalings);
891: PetscArraycpy(shell->scalings,scalings,nmat);
892: return 0;
893: }
895: /*@
896: MatCompositeSetScalings - Sets separate scaling factors for component matrices.
898: Logically Collective on Mat
900: Input Parameters:
901: + mat - the composite matrix
902: - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat)
904: Level: advanced
906: .seealso: MatScale(), MatDiagonalScale(), MATCOMPOSITE
908: @*/
909: PetscErrorCode MatCompositeSetScalings(Mat mat,const PetscScalar *scalings)
910: {
914: PetscUseMethod(mat,"MatCompositeSetScalings_C",(Mat,const PetscScalar*),(mat,scalings));
915: return 0;
916: }
918: static struct _MatOps MatOps_Values = {NULL,
919: NULL,
920: NULL,
921: MatMult_Composite,
922: MatMultAdd_Composite,
923: /* 5*/ MatMultTranspose_Composite,
924: MatMultTransposeAdd_Composite,
925: NULL,
926: NULL,
927: NULL,
928: /* 10*/ NULL,
929: NULL,
930: NULL,
931: NULL,
932: NULL,
933: /* 15*/ NULL,
934: NULL,
935: MatGetDiagonal_Composite,
936: MatDiagonalScale_Composite,
937: NULL,
938: /* 20*/ NULL,
939: MatAssemblyEnd_Composite,
940: NULL,
941: NULL,
942: /* 24*/ NULL,
943: NULL,
944: NULL,
945: NULL,
946: NULL,
947: /* 29*/ NULL,
948: NULL,
949: NULL,
950: NULL,
951: NULL,
952: /* 34*/ NULL,
953: NULL,
954: NULL,
955: NULL,
956: NULL,
957: /* 39*/ NULL,
958: NULL,
959: NULL,
960: NULL,
961: NULL,
962: /* 44*/ NULL,
963: MatScale_Composite,
964: MatShift_Basic,
965: NULL,
966: NULL,
967: /* 49*/ NULL,
968: NULL,
969: NULL,
970: NULL,
971: NULL,
972: /* 54*/ NULL,
973: NULL,
974: NULL,
975: NULL,
976: NULL,
977: /* 59*/ NULL,
978: MatDestroy_Composite,
979: NULL,
980: NULL,
981: NULL,
982: /* 64*/ NULL,
983: NULL,
984: NULL,
985: NULL,
986: NULL,
987: /* 69*/ NULL,
988: NULL,
989: NULL,
990: NULL,
991: NULL,
992: /* 74*/ NULL,
993: NULL,
994: MatSetFromOptions_Composite,
995: NULL,
996: NULL,
997: /* 79*/ NULL,
998: NULL,
999: NULL,
1000: NULL,
1001: NULL,
1002: /* 84*/ NULL,
1003: NULL,
1004: NULL,
1005: NULL,
1006: NULL,
1007: /* 89*/ NULL,
1008: NULL,
1009: NULL,
1010: NULL,
1011: NULL,
1012: /* 94*/ NULL,
1013: NULL,
1014: NULL,
1015: NULL,
1016: NULL,
1017: /*99*/ NULL,
1018: NULL,
1019: NULL,
1020: NULL,
1021: NULL,
1022: /*104*/ NULL,
1023: NULL,
1024: NULL,
1025: NULL,
1026: NULL,
1027: /*109*/ NULL,
1028: NULL,
1029: NULL,
1030: NULL,
1031: NULL,
1032: /*114*/ NULL,
1033: NULL,
1034: NULL,
1035: NULL,
1036: NULL,
1037: /*119*/ NULL,
1038: NULL,
1039: NULL,
1040: NULL,
1041: NULL,
1042: /*124*/ NULL,
1043: NULL,
1044: NULL,
1045: NULL,
1046: NULL,
1047: /*129*/ NULL,
1048: NULL,
1049: NULL,
1050: NULL,
1051: NULL,
1052: /*134*/ NULL,
1053: NULL,
1054: NULL,
1055: NULL,
1056: NULL,
1057: /*139*/ NULL,
1058: NULL,
1059: NULL,
1060: NULL,
1061: NULL,
1062: /*144*/NULL,
1063: NULL,
1064: NULL,
1065: NULL
1066: };
1068: /*MC
1069: MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices.
1070: The matrices need to have a correct size and parallel layout for the sum or product to be valid.
1072: Notes:
1073: to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE);
1075: Level: advanced
1077: .seealso: MatCreateComposite(), MatCompositeSetScalings(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeSetMatStructure(), MatCompositeGetMatStructure(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat()
1078: M*/
1080: PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A)
1081: {
1082: Mat_Composite *b;
1084: PetscNewLog(A,&b);
1085: A->data = (void*)b;
1086: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1088: PetscLayoutSetUp(A->rmap);
1089: PetscLayoutSetUp(A->cmap);
1091: A->assembled = PETSC_TRUE;
1092: A->preallocated = PETSC_TRUE;
1093: b->type = MAT_COMPOSITE_ADDITIVE;
1094: b->scale = 1.0;
1095: b->nmat = 0;
1096: b->merge = PETSC_FALSE;
1097: b->mergetype = MAT_COMPOSITE_MERGE_RIGHT;
1098: b->structure = DIFFERENT_NONZERO_PATTERN;
1099: b->merge_mvctx = PETSC_TRUE;
1101: PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);
1102: PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);
1103: PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);
1104: PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);
1105: PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);
1106: PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);
1107: PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);
1108: PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);
1109: PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);
1110: PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetScalings_C",MatCompositeSetScalings_Composite);
1112: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
1113: return 0;
1114: }