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: }