Actual source code: adamat.c

petsc-3.8.4 2018-03-24
Report Typos and Errors
  1:  #include <petscmat.h>

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

 14: extern PetscErrorCode MatCreateADA(Mat,Vec, Vec, Mat*);

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

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

 36: static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
 37: {

 41:   MatMult_ADA(mat,a,y);
 42:   return(0);
 43: }

 45: static PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
 46: {
 47:   TaoMatADACtx   ctx;
 48:   PetscReal      zero=0.0,one = 1.0;

 52:   MatShellGetContext(M,(void **)&ctx);
 53:   if (!ctx->D2){
 54:     VecDuplicate(D,&ctx->D2);
 55:     VecSet(ctx->D2, zero);
 56:   }
 57:   VecAXPY(ctx->D2, one, D);
 58:   return(0);
 59: }

 61: static PetscErrorCode MatDestroy_ADA(Mat mat)
 62: {
 64:   TaoMatADACtx   ctx;

 67:   ierr=MatShellGetContext(mat,(void **)&ctx);
 68:   ierr=VecDestroy(&ctx->W);
 69:   ierr=VecDestroy(&ctx->W2);
 70:   ierr=VecDestroy(&ctx->ADADiag);
 71:   ierr=MatDestroy(&ctx->A);
 72:   ierr=VecDestroy(&ctx->D1);
 73:   ierr=VecDestroy(&ctx->D2);
 74:   PetscFree(ctx);
 75:   return(0);
 76: }

 78: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
 79: {
 81:   return(0);
 82: }

 84: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
 85: {
 87:   TaoMatADACtx   ctx;

 90:   MatShellGetContext(Y,(void **)&ctx);
 91:   VecShift(ctx->D2,a);
 92:   return(0);
 93: }

 95: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
 96: {
 97:   PetscErrorCode    ierr;
 98:   TaoMatADACtx      ctx;
 99:   Mat               A2;
100:   Vec               D1b=NULL,D2b;

103:   MatShellGetContext(mat,(void **)&ctx);
104:   MatDuplicate(ctx->A,op,&A2);
105:   if (ctx->D1){
106:     VecDuplicate(ctx->D1,&D1b);
107:     VecCopy(ctx->D1,D1b);
108:   }
109:   VecDuplicate(ctx->D2,&D2b);
110:   VecCopy(ctx->D2,D2b);
111:   MatCreateADA(A2,D1b,D2b,M);
112:   if (ctx->D1){
113:     PetscObjectDereference((PetscObject)D1b);
114:   }
115:   PetscObjectDereference((PetscObject)D2b);
116:   PetscObjectDereference((PetscObject)A2);
117:   return(0);
118: }

120: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
121: {
123:   TaoMatADACtx   ctx1,ctx2;

126:   MatShellGetContext(A,(void **)&ctx1);
127:   MatShellGetContext(B,(void **)&ctx2);
128:   VecEqual(ctx1->D2,ctx2->D2,flg);
129:   if (*flg==PETSC_TRUE){
130:     VecEqual(ctx1->D1,ctx2->D1,flg);
131:   }
132:   if (*flg==PETSC_TRUE){
133:     MatEqual(ctx1->A,ctx2->A,flg);
134:   }
135:   return(0);
136: }

138: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
139: {
141:   TaoMatADACtx   ctx;

144:   MatShellGetContext(mat,(void **)&ctx);
145:   VecScale(ctx->D1,a);
146:   if (ctx->D2){
147:     VecScale(ctx->D2,a);
148:   }
149:   return(0);
150: }

152: static PetscErrorCode MatTranspose_ADA(Mat mat,MatReuse reuse,Mat *B)
153: {
155:   TaoMatADACtx   ctx;

158:   MatShellGetContext(mat,(void **)&ctx);
159:   if (reuse == MAT_INITIAL_MATRIX){
160:     MatDuplicate(mat,MAT_COPY_VALUES,B);
161:   } else if (reuse == MAT_REUSE_MATRIX){
162:     MatCopy(mat,*B,SAME_NONZERO_PATTERN);
163:   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Does not support inplace transpose");
164:   return(0);
165: }

167: PetscErrorCode MatADAComputeDiagonal(Mat mat)
168: {
170:   PetscInt       i,m,n,low,high;
171:   PetscScalar    *dtemp,*dptr;
172:   TaoMatADACtx   ctx;

175:   MatShellGetContext(mat,(void **)&ctx);
176:   MatGetOwnershipRange(mat, &low, &high);
177:   MatGetSize(mat,&m,&n);

179:   PetscMalloc1(n,&dtemp );

181:   for (i=0; i<n; i++){
182:     MatGetColumnVector(ctx->A, ctx->W, i);
183:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
184:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
185:   }
186:   for (i=0; i<n; i++){
187:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
188:   }

190:   VecGetArray(ctx->ADADiag,&dptr);
191:   for (i=low; i<high; i++){
192:     dptr[i-low]= dtemp[i];
193:   }
194:   VecRestoreArray(ctx->ADADiag,&dptr);
195:   PetscFree(dtemp);
196:   return(0);
197: }

199: static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
200: {
201:   PetscErrorCode  ierr;
202:   PetscReal       one=1.0;
203:   TaoMatADACtx    ctx;

206:   MatShellGetContext(mat,(void **)&ctx);
207:   MatADAComputeDiagonal(mat);
208:   ierr=VecCopy(ctx->ADADiag,v);
209:   if (ctx->D2){
210:     ierr=VecAXPY(v, one, ctx->D2);
211:   }
212:   return(0);
213: }

215: static PetscErrorCode MatCreateSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
216: {
217:   PetscErrorCode    ierr;
218:   PetscInt          low,high;
219:   IS                ISrow;
220:   Vec               D1,D2;
221:   Mat               Atemp;
222:   TaoMatADACtx      ctx;
223:   PetscBool         isequal;

226:   ISEqual(isrow,iscol,&isequal);
227:   if (!isequal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices");
228:   MatShellGetContext(mat,(void **)&ctx);

230:   MatGetOwnershipRange(ctx->A,&low,&high);
231:   ISCreateStride(PetscObjectComm((PetscObject)mat),high-low,low,1,&ISrow);
232:   MatCreateSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
233:   ISDestroy(&ISrow);

235:   if (ctx->D1){
236:     ierr=VecDuplicate(ctx->D1,&D1);
237:     ierr=VecCopy(ctx->D1,D1);
238:   } else {
239:     D1 = NULL;
240:   }

242:   if (ctx->D2){
243:     Vec D2sub;

245:     ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);
246:     ierr=VecDuplicate(D2sub,&D2);
247:     ierr=VecCopy(D2sub,D2);
248:     ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);
249:   } else {
250:     D2 = NULL;
251:   }

253:   MatCreateADA(Atemp,D1,D2,newmat);
254:   MatShellGetContext(*newmat,(void **)&ctx);
255:   PetscObjectDereference((PetscObject)Atemp);
256:   if (ctx->D1){
257:     PetscObjectDereference((PetscObject)D1);
258:   }
259:   if (ctx->D2){
260:     PetscObjectDereference((PetscObject)D2);
261:   }
262:   return(0);
263: }

265: static PetscErrorCode MatCreateSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
266: {
268:   PetscInt       i;

271:   if (scall == MAT_INITIAL_MATRIX) {
272:     PetscCalloc1(n+1,B );
273:   }
274:   for ( i=0; i<n; i++ ) {
275:     MatCreateSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
276:   }
277:   return(0);
278: }

280: static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
281: {
283:   PetscInt       low,high;
284:   PetscScalar    zero=0.0,one=1.0;

287:   ierr=VecSet(Y, zero);
288:   ierr=VecGetOwnershipRange(Y,&low,&high);
289:   if (col>=low && col<high){
290:     ierr=VecSetValue(Y,col,one,INSERT_VALUES);
291:   }
292:   ierr=VecAssemblyBegin(Y);
293:   ierr=VecAssemblyEnd(Y);
294:   ierr=MatMult_ADA(mat,Y,Y);
295:   return(0);
296: }

298: PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
299: {
301:   PetscMPIInt    size;
302:   PetscBool      sametype, issame, isdense, isseqdense;
303:   TaoMatADACtx   ctx;

306:   MatShellGetContext(mat,(void **)&ctx);
307:   MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);

309:   PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
310:   PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
311:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
312:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);

314:   if (sametype || issame) {
315:     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
316:   } else if (isdense) {
317:     PetscInt          i,j,low,high,m,n,M,N;
318:     const PetscScalar *dptr;
319:     Vec               X;

321:     VecDuplicate(ctx->D2,&X);
322:     ierr=MatGetSize(mat,&M,&N);
323:     ierr=MatGetLocalSize(mat,&m,&n);
324:     MatCreateDense(PetscObjectComm((PetscObject)mat),m,m,N,N,NULL,NewMat);
325:     MatGetOwnershipRange(*NewMat,&low,&high);
326:     for (i=0;i<M;i++){
327:       MatGetColumnVector_ADA(mat,X,i);
328:       VecGetArrayRead(X,&dptr);
329:       for (j=0; j<high-low; j++){
330:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
331:       }
332:       VecRestoreArrayRead(X,&dptr);
333:     }
334:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
335:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
336:     VecDestroy(&X);
337:   } else if (isseqdense && size==1){
338:     PetscInt          i,j,low,high,m,n,M,N;
339:     const PetscScalar *dptr;
340:     Vec               X;

342:     VecDuplicate(ctx->D2,&X);
343:     MatGetSize(mat,&M,&N);
344:     MatGetLocalSize(mat,&m,&n);
345:     MatCreateSeqDense(PetscObjectComm((PetscObject)mat),N,N,NULL,NewMat);
346:     MatGetOwnershipRange(*NewMat,&low,&high);
347:     for (i=0;i<M;i++){
348:       MatGetColumnVector_ADA(mat,X,i);
349:       VecGetArrayRead(X,&dptr);
350:       for (j=0; j<high-low; j++){
351:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
352:       }
353:       VecRestoreArrayRead(X,&dptr);
354:     }
355:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
356:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
357:     VecDestroy(&X);
358:   } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
359:   return(0);
360: }

362: static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
363: {
365:   TaoMatADACtx   ctx;

368:   MatShellGetContext(mat,(void **)&ctx);
369:   if (type == NORM_FROBENIUS) {
370:     *norm = 1.0;
371:   } else if (type == NORM_1 || type == NORM_INFINITY) {
372:     *norm = 1.0;
373:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
374:   return(0);
375: }

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

380:    Collective on matrix

382:    Input Parameters:
383: +  mat - matrix of arbitrary type
384: .  d1 - A vector with diagonal elements of D1
385: -  d2 - A vector

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

390:    Notes:
391:    The user provides the input data and is responsible for destroying
392:    this data after matrix J has been destroyed.
393:    The operation MatMult(A,D2,D1) must be well defined.
394:    Before calling the operation MatGetDiagonal(), the function
395:    MatADAComputeDiagonal() must be called.  The matrices A and D1 must
396:    be the same during calls to MatADAComputeDiagonal() and
397:    MatGetDiagonal().

399:    Level: developer

401: .seealso: MatCreate()
402: @*/
403: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
404: {
405:   MPI_Comm       comm=PetscObjectComm((PetscObject)mat);
406:   TaoMatADACtx   ctx;
408:   PetscInt       nloc,n;

411:   PetscNew(&ctx);
412:   ctx->A=mat;
413:   ctx->D1=d1;
414:   ctx->D2=d2;
415:   if (d1){
416:     VecDuplicate(d1,&ctx->W);
417:     PetscObjectReference((PetscObject)d1);
418:   } else {
419:     ctx->W = NULL;
420:   }
421:   if (d2){
422:     VecDuplicate(d2,&ctx->W2);
423:     VecDuplicate(d2,&ctx->ADADiag);
424:      PetscObjectReference((PetscObject)d2);
425:   } else {
426:     ctx->W2      = NULL;
427:     ctx->ADADiag = NULL;
428:   }

430:   ctx->GotDiag=0;
431:    PetscObjectReference((PetscObject)mat);

433:   ierr=VecGetLocalSize(d2,&nloc);
434:   ierr=VecGetSize(d2,&n);

436:   MatCreateShell(comm,nloc,nloc,n,n,ctx,J);

438:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
439:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
440:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
441:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
442:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
443:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
444:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
445:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
446:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
447:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
448:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_ADA);
449:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
450:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
451:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_ADA);

453:   PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
454:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));

456:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
457:   return(0);
458: }