Actual source code: mcomposite.c
1: #define PETSCMAT_DLL
3: #include include/private/matimpl.h
5: typedef struct _Mat_CompositeLink *Mat_CompositeLink;
6: struct _Mat_CompositeLink {
7: Mat mat;
8: Mat_CompositeLink next;
9: };
10:
11: typedef struct {
12: Mat_CompositeLink head;
13: Vec work;
14: } Mat_Composite;
18: PetscErrorCode MatDestroy_Composite(Mat mat)
19: {
20: PetscErrorCode ierr;
21: Mat_Composite *shell = (Mat_Composite*)mat->data;
22: Mat_CompositeLink next = shell->head;
25: while (next) {
26: MatDestroy(next->mat);
27: next = next->next;
28: }
29: if (shell->work) {VecDestroy(shell->work);}
30: PetscFree(shell);
31: mat->data = 0;
32: return(0);
33: }
37: PetscErrorCode MatMult_Composite(Mat A,Vec x,Vec y)
38: {
39: Mat_Composite *shell = (Mat_Composite*)A->data;
40: Mat_CompositeLink next = shell->head;
41: PetscErrorCode ierr;
44: if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
45: MatMult(next->mat,x,y);
46: while ((next = next->next)) {
47: MatMultAdd(next->mat,x,y,y);
48: }
49: return(0);
50: }
54: PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y)
55: {
56: Mat_Composite *shell = (Mat_Composite*)A->data;
57: Mat_CompositeLink next = shell->head;
58: PetscErrorCode ierr;
61: if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
62: MatMultTranspose(next->mat,x,y);
63: while ((next = next->next)) {
64: MatMultTransposeAdd(next->mat,x,y,y);
65: }
66: return(0);
67: }
71: PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v)
72: {
73: Mat_Composite *shell = (Mat_Composite*)A->data;
74: Mat_CompositeLink next = shell->head;
75: PetscErrorCode ierr;
78: if (!next) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()");
79: MatGetDiagonal(next->mat,v);
80: if (next->next && !shell->work) {
81: VecDuplicate(v,&shell->work);
82: }
83: while ((next = next->next)) {
84: MatGetDiagonal(next->mat,shell->work);
85: VecAXPY(v,1.0,shell->work);
86: }
87: return(0);
88: }
92: PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t)
93: {
95: return(0);
96: }
99: static struct _MatOps MatOps_Values = {0,
100: 0,
101: 0,
102: MatMult_Composite,
103: 0,
104: /* 5*/ MatMultTranspose_Composite,
105: 0,
106: 0,
107: 0,
108: 0,
109: /*10*/ 0,
110: 0,
111: 0,
112: 0,
113: 0,
114: /*15*/ 0,
115: 0,
116: MatGetDiagonal_Composite,
117: 0,
118: 0,
119: /*20*/ 0,
120: MatAssemblyEnd_Composite,
121: 0,
122: 0,
123: 0,
124: /*25*/ 0,
125: 0,
126: 0,
127: 0,
128: 0,
129: /*30*/ 0,
130: 0,
131: 0,
132: 0,
133: 0,
134: /*35*/ 0,
135: 0,
136: 0,
137: 0,
138: 0,
139: /*40*/ 0,
140: 0,
141: 0,
142: 0,
143: 0,
144: /*45*/ 0,
145: 0,
146: 0,
147: 0,
148: 0,
149: /*50*/ 0,
150: 0,
151: 0,
152: 0,
153: 0,
154: /*55*/ 0,
155: 0,
156: 0,
157: 0,
158: 0,
159: /*60*/ 0,
160: MatDestroy_Composite,
161: 0,
162: 0,
163: 0,
164: /*65*/ 0,
165: 0,
166: 0,
167: 0,
168: 0,
169: /*70*/ 0,
170: 0,
171: 0,
172: 0,
173: 0,
174: /*75*/ 0,
175: 0,
176: 0,
177: 0,
178: 0,
179: /*80*/ 0,
180: 0,
181: 0,
182: 0,
183: 0,
184: /*85*/ 0,
185: 0,
186: 0,
187: 0,
188: 0,
189: /*90*/ 0,
190: 0,
191: 0,
192: 0,
193: 0,
194: /*95*/ 0,
195: 0,
196: 0,
197: 0};
199: /*MC
200: MATCOMPOSITE - A matrix defined by the sum of one or more matrices (all matrices are of same size and parallel layout)
202: Level: advanced
204: .seealso: MatCreateComposite
205: M*/
210: PetscErrorCode MatCreate_Composite(Mat A)
211: {
212: Mat_Composite *b;
216: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
218: PetscNew(Mat_Composite,&b);
219: A->data = (void*)b;
221: PetscMapInitialize(A->comm,&A->rmap);
222: PetscMapInitialize(A->comm,&A->cmap);
224: A->assembled = PETSC_TRUE;
225: A->preallocated = PETSC_FALSE;
226: PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);
227: return(0);
228: }
233: /*@C
234: MatCreateComposite - Creates a matrix as the sum of zero or more matrices
236: Collective on MPI_Comm
238: Input Parameters:
239: + comm - MPI communicator
240: . nmat - number of matrices to put in
241: - mats - the matrices
243: Output Parameter:
244: . A - the matrix
246: Level: advanced
248: Notes:
249: Alternative construction
250: $ MatCreate(comm,&mat);
251: $ MatSetSizes(mat,m,n,M,N);
252: $ MatSetType(mat,MATCOMPOSITE);
253: $ MatCompositeAddMat(mat,mats[0]);
254: $ ....
255: $ MatCompositeAddMat(mat,mats[nmat-1]);
257: .seealso: MatDestroy(), MatMult(), MatCompositeAddMat()
259: @*/
260: PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat)
261: {
263: PetscInt m,n,M,N,i;
266: if (nmat < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix");
269: MatGetLocalSize(mats[0],&m,&n);
270: MatGetSize(mats[0],&M,&N);
271: MatCreate(comm,mat);
272: MatSetSizes(*mat,m,n,M,N);
273: MatSetType(*mat,MATCOMPOSITE);
274: for (i=0; i<nmat; i++) {
275: MatCompositeAddMat(*mat,mats[i]);
276: }
277: return(0);
278: }
282: /*@
283: MatCompositeAddMat - add another matrix to a composite matrix
285: Collective on Mat
287: Input Parameters:
288: + mat - the composite matrix
289: - smat - the partial matrix
291: Level: advanced
293: .seealso: MatCreateComposite()
294: @*/
295: PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat)
296: {
297: Mat_Composite *shell = (Mat_Composite*)mat->data;
298: PetscErrorCode ierr;
299: Mat_CompositeLink ilink,next;
304: PetscNew(struct _Mat_CompositeLink,&ilink);
305: ilink->next = 0;
306: ilink->mat = smat;
307: PetscObjectReference((PetscObject)smat);
309: next = shell->head;
310: if (!next) {
311: shell->head = ilink;
312: } else {
313: while (next->next) {
314: next = next->next;
315: }
316: next->next = ilink;
317: }
318: return(0);
319: }