Actual source code: adamat.c

  1: #include <petscmat.h>

  3: PETSC_INTERN PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);

  5: typedef struct {
  6:   Mat      A;
  7:   Vec      D1;
  8:   Vec      D2;
  9:   Vec      W;
 10:   Vec      W2;
 11:   Vec      ADADiag;
 12:   PetscInt GotDiag;
 13: } _p_TaoMatADACtx;
 14: typedef  _p_TaoMatADACtx* TaoMatADACtx;

 16: static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
 17: {
 18:   TaoMatADACtx   ctx;
 19:   PetscReal      one = 1.0;

 21:   MatShellGetContext(mat,&ctx);
 22:   MatMult(ctx->A,a,ctx->W);
 23:   if (ctx->D1) {
 24:     VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
 25:   }
 26:   MatMultTranspose(ctx->A,ctx->W,y);
 27:   if (ctx->D2) {
 28:     VecPointwiseMult(ctx->W2, ctx->D2, a);
 29:     VecAXPY(y, one, ctx->W2);
 30:   }
 31:   return 0;
 32: }

 34: static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
 35: {
 36:   MatMult_ADA(mat,a,y);
 37:   return 0;
 38: }

 40: static PetscErrorCode MatDiagonalSet_ADA(Mat M,Vec D, InsertMode mode)
 41: {
 42:   TaoMatADACtx   ctx;
 43:   PetscReal      zero=0.0,one = 1.0;

 46:   MatShellGetContext(M,&ctx);
 47:   if (!ctx->D2) {
 48:     VecDuplicate(D,&ctx->D2);
 49:     VecSet(ctx->D2, zero);
 50:   }
 51:   VecAXPY(ctx->D2, one, D);
 52:   return 0;
 53: }

 55: static PetscErrorCode MatDestroy_ADA(Mat mat)
 56: {
 57:   TaoMatADACtx   ctx;

 59:   MatShellGetContext(mat,&ctx);
 60:   VecDestroy(&ctx->W);
 61:   VecDestroy(&ctx->W2);
 62:   VecDestroy(&ctx->ADADiag);
 63:   MatDestroy(&ctx->A);
 64:   VecDestroy(&ctx->D1);
 65:   VecDestroy(&ctx->D2);
 66:   PetscFree(ctx);
 67:   return 0;
 68: }

 70: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
 71: {
 72:   return 0;
 73: }

 75: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
 76: {
 77:   TaoMatADACtx   ctx;

 79:   MatShellGetContext(Y,&ctx);
 80:   VecShift(ctx->D2,a);
 81:   return 0;
 82: }

 84: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
 85: {
 86:   TaoMatADACtx      ctx;
 87:   Mat               A2;
 88:   Vec               D1b=NULL,D2b;

 90:   MatShellGetContext(mat,&ctx);
 91:   MatDuplicate(ctx->A,op,&A2);
 92:   if (ctx->D1) {
 93:     VecDuplicate(ctx->D1,&D1b);
 94:     VecCopy(ctx->D1,D1b);
 95:   }
 96:   VecDuplicate(ctx->D2,&D2b);
 97:   VecCopy(ctx->D2,D2b);
 98:   MatCreateADA(A2,D1b,D2b,M);
 99:   if (ctx->D1) {
100:     PetscObjectDereference((PetscObject)D1b);
101:   }
102:   PetscObjectDereference((PetscObject)D2b);
103:   PetscObjectDereference((PetscObject)A2);
104:   return 0;
105: }

107: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
108: {
109:   TaoMatADACtx   ctx1,ctx2;

111:   MatShellGetContext(A,&ctx1);
112:   MatShellGetContext(B,&ctx2);
113:   VecEqual(ctx1->D2,ctx2->D2,flg);
114:   if (*flg==PETSC_TRUE) {
115:     VecEqual(ctx1->D1,ctx2->D1,flg);
116:   }
117:   if (*flg==PETSC_TRUE) {
118:     MatEqual(ctx1->A,ctx2->A,flg);
119:   }
120:   return 0;
121: }

123: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
124: {
125:   TaoMatADACtx   ctx;

127:   MatShellGetContext(mat,&ctx);
128:   VecScale(ctx->D1,a);
129:   if (ctx->D2) {
130:     VecScale(ctx->D2,a);
131:   }
132:   return 0;
133: }

135: static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B)
136: {
137:   TaoMatADACtx   ctx;

139:   MatShellGetContext(mat,&ctx);
140:   if (reuse == MAT_INITIAL_MATRIX) {
141:     MatDuplicate(mat,MAT_COPY_VALUES,B);
142:   } else if (reuse == MAT_REUSE_MATRIX) {
143:     MatCopy(mat,*B,SAME_NONZERO_PATTERN);
144:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose");
145:   return 0;
146: }

148: static PetscErrorCode MatADAComputeDiagonal(Mat mat)
149: {
150:   PetscInt       i,m,n,low,high;
151:   PetscScalar    *dtemp,*dptr;
152:   TaoMatADACtx   ctx;

154:   MatShellGetContext(mat,&ctx);
155:   MatGetOwnershipRange(mat, &low, &high);
156:   MatGetSize(mat,&m,&n);

158:   PetscMalloc1(n,&dtemp);
159:   for (i=0; i<n; i++) {
160:     MatGetColumnVector(ctx->A, ctx->W, i);
161:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
162:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
163:   }
164:   for (i=0; i<n; i++) {
165:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
166:   }

168:   VecGetArray(ctx->ADADiag,&dptr);
169:   for (i=low; i<high; i++) {
170:     dptr[i-low]= dtemp[i];
171:   }
172:   VecRestoreArray(ctx->ADADiag,&dptr);
173:   PetscFree(dtemp);
174:   return 0;
175: }

177: static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
178: {
179:   PetscReal       one=1.0;
180:   TaoMatADACtx    ctx;

182:   MatShellGetContext(mat,&ctx);
183:   MatADAComputeDiagonal(mat);
184:   VecCopy(ctx->ADADiag,v);
185:   if (ctx->D2) {
186:     VecAXPY(v, one, ctx->D2);
187:   }
188:   return 0;
189: }

191: static PetscErrorCode MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
192: {
193:   PetscInt          low,high;
194:   IS                ISrow;
195:   Vec               D1,D2;
196:   Mat               Atemp;
197:   TaoMatADACtx      ctx;
198:   PetscBool         isequal;

200:   ISEqual(isrow,iscol,&isequal);
202:   MatShellGetContext(mat,&ctx);

204:   MatGetOwnershipRange(ctx->A,&low,&high);
205:   ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);
206:   MatCreateSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
207:   ISDestroy(&ISrow);

209:   if (ctx->D1) {
210:     VecDuplicate(ctx->D1,&D1);
211:     VecCopy(ctx->D1,D1);
212:   } else {
213:     D1 = NULL;
214:   }

216:   if (ctx->D2) {
217:     Vec D2sub;

219:     VecGetSubVector(ctx->D2,isrow,&D2sub);
220:     VecDuplicate(D2sub,&D2);
221:     VecCopy(D2sub,D2);
222:     VecRestoreSubVector(ctx->D2,isrow,&D2sub);
223:   } else {
224:     D2 = NULL;
225:   }

227:   MatCreateADA(Atemp,D1,D2,newmat);
228:   MatShellGetContext(*newmat,&ctx);
229:   PetscObjectDereference((PetscObject)Atemp);
230:   if (ctx->D1) {
231:     PetscObjectDereference((PetscObject)D1);
232:   }
233:   if (ctx->D2) {
234:     PetscObjectDereference((PetscObject)D2);
235:   }
236:   return 0;
237: }

239: static PetscErrorCode MatCreateSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
240: {
241:   PetscInt       i;

243:   if (scall == MAT_INITIAL_MATRIX) {
244:     PetscCalloc1(n+1,B);
245:   }
246:   for (i=0; i<n; i++) {
247:     MatCreateSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
248:   }
249:   return 0;
250: }

252: static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
253: {
254:   PetscInt       low,high;
255:   PetscScalar    zero=0.0,one=1.0;

257:   VecSet(Y, zero);
258:   VecGetOwnershipRange(Y,&low,&high);
259:   if (col>=low && col<high) {
260:     VecSetValue(Y,col,one,INSERT_VALUES);
261:   }
262:   VecAssemblyBegin(Y);
263:   VecAssemblyEnd(Y);
264:   MatMult_ADA(mat,Y,Y);
265:   return 0;
266: }

268: PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
269: {
270:   PetscMPIInt    size;
271:   PetscBool      sametype, issame, isdense, isseqdense;
272:   TaoMatADACtx   ctx;

274:   MatShellGetContext(mat,&ctx);
275:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);

277:   PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
278:   PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
279:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
280:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);

282:   if (sametype || issame) {
283:     MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
284:   } else if (isdense) {
285:     PetscInt          i,j,low,high,m,n,M,N;
286:     const PetscScalar *dptr;
287:     Vec               X;

289:     VecDuplicate(ctx->D2,&X);
290:     MatGetSize(mat,&M,&N);
291:     MatGetLocalSize(mat,&m,&n);
292:     MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);
293:     MatGetOwnershipRange(*NewMat,&low,&high);
294:     for (i=0;i<M;i++) {
295:       MatGetColumnVector_ADA(mat,X,i);
296:       VecGetArrayRead(X,&dptr);
297:       for (j=0; j<high-low; j++) {
298:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
299:       }
300:       VecRestoreArrayRead(X,&dptr);
301:     }
302:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
303:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
304:     VecDestroy(&X);
305:   } else if (isseqdense && size==1) {
306:     PetscInt          i,j,low,high,m,n,M,N;
307:     const PetscScalar *dptr;
308:     Vec               X;

310:     VecDuplicate(ctx->D2,&X);
311:     MatGetSize(mat,&M,&N);
312:     MatGetLocalSize(mat,&m,&n);
313:     MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);
314:     MatGetOwnershipRange(*NewMat,&low,&high);
315:     for (i=0;i<M;i++) {
316:       MatGetColumnVector_ADA(mat,X,i);
317:       VecGetArrayRead(X,&dptr);
318:       for (j=0; j<high-low; j++) {
319:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
320:       }
321:       VecRestoreArrayRead(X,&dptr);
322:     }
323:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
324:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
325:     VecDestroy(&X);
326:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"No support to convert objects to that type");
327:   return 0;
328: }

330: static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
331: {
332:   TaoMatADACtx   ctx;

334:   MatShellGetContext(mat,&ctx);
335:   if (type == NORM_FROBENIUS) {
336:     *norm = 1.0;
337:   } else if (type == NORM_1 || type == NORM_INFINITY) {
338:     *norm = 1.0;
339:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
340:   return 0;
341: }

343: /*@C
344:    MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal

346:    Collective on matrix

348:    Input Parameters:
349: +  mat - matrix of arbitrary type
350: .  d1 - A vector defining a diagonal matrix
351: -  d2 - A vector defining a diagonal matrix

353:    Output Parameters:
354: .  J - New matrix whose operations are defined in terms of mat, D1, and D2.

356:    Notes:
357:    The user provides the input data and is responsible for destroying
358:    this data after matrix J has been destroyed.

360:    Level: developer

362: .seealso: MatCreate()
363: @*/
364: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
365: {
366:   MPI_Comm       comm = PetscObjectComm((PetscObject)mat);
367:   TaoMatADACtx   ctx;
368:   PetscInt       nloc,n;

370:   PetscNew(&ctx);
371:   ctx->A=mat;
372:   ctx->D1=d1;
373:   ctx->D2=d2;
374:   if (d1) {
375:     VecDuplicate(d1,&ctx->W);
376:     PetscObjectReference((PetscObject)d1);
377:   } else {
378:     ctx->W = NULL;
379:   }
380:   if (d2) {
381:     VecDuplicate(d2,&ctx->W2);
382:     VecDuplicate(d2,&ctx->ADADiag);
383:     PetscObjectReference((PetscObject)d2);
384:   } else {
385:     ctx->W2      = NULL;
386:     ctx->ADADiag = NULL;
387:   }

389:   ctx->GotDiag = 0;
390:   PetscObjectReference((PetscObject)mat);

392:   VecGetLocalSize(d2,&nloc);
393:   VecGetSize(d2,&n);

395:   MatCreateShell(comm,nloc,nloc,n,n,ctx,J);
396:   MatShellSetManageScalingShifts(*J);
397:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
398:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
399:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
400:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
401:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
402:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
403:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
404:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
405:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
406:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
407:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);
408:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
409:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
410:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);

412:   PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
413:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));

415:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
416:   return 0;
417: }