Actual source code: adamat.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/matimpl.h>              /*I  "mat.h"  I*/
  2: #include <petsc/private/vecimpl.h>

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

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

 19: static PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
 20: {
 21:   TaoMatADACtx   ctx;
 22:   PetscReal      one = 1.0;

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

 41: static PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
 42: {

 46:   MatMult_ADA(mat,a,y);
 47:   return(0);
 48: }

 52: static PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
 53: {
 54:   TaoMatADACtx   ctx;
 55:   PetscReal      zero=0.0,one = 1.0;

 59:   MatShellGetContext(M,(void **)&ctx);
 60:   if (!ctx->D2){
 61:     VecDuplicate(D,&ctx->D2);
 62:     VecSet(ctx->D2, zero);
 63:   }
 64:   VecAXPY(ctx->D2, one, D);
 65:   return(0);
 66: }

 70: static PetscErrorCode MatDestroy_ADA(Mat mat)
 71: {
 73:   TaoMatADACtx   ctx;

 76:   ierr=MatShellGetContext(mat,(void **)&ctx);
 77:   ierr=VecDestroy(&ctx->W);
 78:   ierr=VecDestroy(&ctx->W2);
 79:   ierr=VecDestroy(&ctx->ADADiag);
 80:   ierr=MatDestroy(&ctx->A);
 81:   ierr=VecDestroy(&ctx->D1);
 82:   ierr=VecDestroy(&ctx->D2);
 83:   PetscFree(ctx);
 84:   return(0);
 85: }

 89: static PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
 90: {
 92:   return(0);
 93: }

 97: static PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
 98: {
100:   TaoMatADACtx   ctx;

103:   MatShellGetContext(Y,(void **)&ctx);
104:   VecShift(ctx->D2,a);
105:   return(0);
106: }

110: static PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
111: {
112:   PetscErrorCode    ierr;
113:   TaoMatADACtx      ctx;
114:   Mat               A2;
115:   Vec               D1b=NULL,D2b;

118:   MatShellGetContext(mat,(void **)&ctx);
119:   MatDuplicate(ctx->A,op,&A2);
120:   if (ctx->D1){
121:     VecDuplicate(ctx->D1,&D1b);
122:     VecCopy(ctx->D1,D1b);
123:   }
124:   VecDuplicate(ctx->D2,&D2b);
125:   VecCopy(ctx->D2,D2b);
126:   MatCreateADA(A2,D1b,D2b,M);
127:   if (ctx->D1){
128:     PetscObjectDereference((PetscObject)D1b);
129:   }
130:   PetscObjectDereference((PetscObject)D2b);
131:   PetscObjectDereference((PetscObject)A2);
132:   return(0);
133: }

137: static PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
138: {
140:   TaoMatADACtx   ctx1,ctx2;

143:   MatShellGetContext(A,(void **)&ctx1);
144:   MatShellGetContext(B,(void **)&ctx2);
145:   VecEqual(ctx1->D2,ctx2->D2,flg);
146:   if (*flg==PETSC_TRUE){
147:     VecEqual(ctx1->D1,ctx2->D1,flg);
148:   }
149:   if (*flg==PETSC_TRUE){
150:     MatEqual(ctx1->A,ctx2->A,flg);
151:   }
152:   return(0);
153: }

157: static PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
158: {
160:   TaoMatADACtx   ctx;

163:   MatShellGetContext(mat,(void **)&ctx);
164:   VecScale(ctx->D1,a);
165:   if (ctx->D2){
166:     VecScale(ctx->D2,a);
167:   }
168:   return(0);
169: }

173: static PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
174: {
176:   TaoMatADACtx   ctx;

179:   if (*B){
180:     MatShellGetContext(mat,(void **)&ctx);
181:     MatDuplicate(mat,MAT_COPY_VALUES,B);
182:   }
183:   return(0);
184: }

188: PetscErrorCode MatADAComputeDiagonal(Mat mat)
189: {
191:   PetscInt       i,m,n,low,high;
192:   PetscScalar    *dtemp,*dptr;
193:   TaoMatADACtx   ctx;

196:   MatShellGetContext(mat,(void **)&ctx);
197:   MatGetOwnershipRange(mat, &low, &high);
198:   MatGetSize(mat,&m,&n);

200:   PetscMalloc1(n,&dtemp );

202:   for (i=0; i<n; i++){
203:     MatGetColumnVector(ctx->A, ctx->W, i);
204:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
205:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
206:   }
207:   for (i=0; i<n; i++){
208:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
209:   }

211:   VecGetArray(ctx->ADADiag,&dptr);
212:   for (i=low; i<high; i++){
213:     dptr[i-low]= dtemp[i];
214:   }
215:   VecRestoreArray(ctx->ADADiag,&dptr);
216:   PetscFree(dtemp);
217:   return(0);
218: }

222: static PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
223: {
224:   PetscErrorCode  ierr;
225:   PetscReal       one=1.0;
226:   TaoMatADACtx    ctx;

229:   MatShellGetContext(mat,(void **)&ctx);
230:   MatADAComputeDiagonal(mat);
231:   ierr=VecCopy(ctx->ADADiag,v);
232:   if (ctx->D2){
233:     ierr=VecAXPY(v, one, ctx->D2);
234:   }
235:   return(0);
236: }

240: static PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
241: {
242:   PetscErrorCode    ierr;
243:   PetscInt          low,high;
244:   IS                ISrow;
245:   Vec               D1,D2;
246:   Mat               Atemp;
247:   TaoMatADACtx      ctx;
248:   PetscBool         isequal;

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

255:   MatGetOwnershipRange(ctx->A,&low,&high);
256:   ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
257:   MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
258:   ISDestroy(&ISrow);

260:   if (ctx->D1){
261:     ierr=VecDuplicate(ctx->D1,&D1);
262:     ierr=VecCopy(ctx->D1,D1);
263:   } else {
264:     D1 = NULL;
265:   }

267:   if (ctx->D2){
268:     Vec D2sub;

270:     ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);
271:     ierr=VecDuplicate(D2sub,&D2);
272:     ierr=VecCopy(D2sub,D2);
273:     ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);
274:   } else {
275:     D2 = NULL;
276:   }

278:   MatCreateADA(Atemp,D1,D2,newmat);
279:   MatShellGetContext(*newmat,(void **)&ctx);
280:   PetscObjectDereference((PetscObject)Atemp);
281:   if (ctx->D1){
282:     PetscObjectDereference((PetscObject)D1);
283:   }
284:   if (ctx->D2){
285:     PetscObjectDereference((PetscObject)D2);
286:   }
287:   return(0);
288: }

292: static PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
293: {
295:   PetscInt       i;

298:   if (scall == MAT_INITIAL_MATRIX) {
299:     PetscMalloc1(n+1,B );
300:   }
301:   for ( i=0; i<n; i++ ) {
302:     MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
303:   }
304:   return(0);
305: }

309: static PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
310: {
312:   PetscInt       low,high;
313:   PetscScalar    zero=0.0,one=1.0;

316:   ierr=VecSet(Y, zero);
317:   ierr=VecGetOwnershipRange(Y,&low,&high);
318:   if (col>=low && col<high){
319:     ierr=VecSetValue(Y,col,one,INSERT_VALUES);
320:   }
321:   ierr=VecAssemblyBegin(Y);
322:   ierr=VecAssemblyEnd(Y);
323:   ierr=MatMult_ADA(mat,Y,Y);
324:   return(0);
325: }

329: PETSC_INTERN PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
330: {
332:   PetscMPIInt    size;
333:   PetscBool      sametype, issame, isdense, isseqdense;
334:   TaoMatADACtx   ctx;

337:   MatShellGetContext(mat,(void **)&ctx);
338:   MPI_Comm_size(((PetscObject)mat)->comm,&size);

340:   PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
341:   PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
342:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
343:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);

345:   if (sametype || issame) {
346:     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
347:   } else if (isdense) {
348:     PetscInt          i,j,low,high,m,n,M,N;
349:     const PetscScalar *dptr;
350:     Vec               X;

352:     VecDuplicate(ctx->D2,&X);
353:     ierr=MatGetSize(mat,&M,&N);
354:     ierr=MatGetLocalSize(mat,&m,&n);
355:     MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);
356:     MatGetOwnershipRange(*NewMat,&low,&high);
357:     for (i=0;i<M;i++){
358:       MatGetColumnVector_ADA(mat,X,i);
359:       VecGetArrayRead(X,&dptr);
360:       for (j=0; j<high-low; j++){
361:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
362:       }
363:       VecRestoreArrayRead(X,&dptr);
364:     }
365:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
366:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
367:     VecDestroy(&X);
368:   } else if (isseqdense && size==1){
369:     PetscInt          i,j,low,high,m,n,M,N;
370:     const PetscScalar *dptr;
371:     Vec               X;

373:     VecDuplicate(ctx->D2,&X);
374:     MatGetSize(mat,&M,&N);
375:     MatGetLocalSize(mat,&m,&n);
376:     MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);
377:     MatGetOwnershipRange(*NewMat,&low,&high);
378:     for (i=0;i<M;i++){
379:       MatGetColumnVector_ADA(mat,X,i);
380:       VecGetArrayRead(X,&dptr);
381:       for (j=0; j<high-low; j++){
382:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
383:       }
384:       VecRestoreArrayRead(X,&dptr);
385:     }
386:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
387:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
388:     VecDestroy(&X);
389:   } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
390:   return(0);
391: }

395: static PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
396: {
398:   TaoMatADACtx   ctx;

401:   MatShellGetContext(mat,(void **)&ctx);
402:   if (type == NORM_FROBENIUS) {
403:     *norm = 1.0;
404:   } else if (type == NORM_1 || type == NORM_INFINITY) {
405:     *norm = 1.0;
406:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
407:   return(0);
408: }

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

415:    Collective on matrix

417:    Input Parameters:
418: +  mat - matrix of arbitrary type
419: .  d1 - A vector with diagonal elements of D1
420: -  d2 - A vector

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

425:    Notes:
426:    The user provides the input data and is responsible for destroying
427:    this data after matrix J has been destroyed.
428:    The operation MatMult(A,D2,D1) must be well defined.
429:    Before calling the operation MatGetDiagonal(), the function
430:    MatADAComputeDiagonal() must be called.  The matrices A and D1 must
431:    be the same during calls to MatADAComputeDiagonal() and
432:    MatGetDiagonal().

434:    Level: developer

436: .seealso: MatCreate()
437: @*/
438: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
439: {
440:   MPI_Comm       comm=((PetscObject)mat)->comm;
441:   TaoMatADACtx   ctx;
443:   PetscInt       nloc,n;

446:   PetscNew(&ctx);
447:   ctx->A=mat;
448:   ctx->D1=d1;
449:   ctx->D2=d2;
450:   if (d1){
451:     VecDuplicate(d1,&ctx->W);
452:     PetscObjectReference((PetscObject)d1);
453:   } else {
454:     ctx->W = NULL;
455:   }
456:   if (d2){
457:     VecDuplicate(d2,&ctx->W2);
458:     VecDuplicate(d2,&ctx->ADADiag);
459:      PetscObjectReference((PetscObject)d2);
460:   } else {
461:     ctx->W2      = NULL;
462:     ctx->ADADiag = NULL;
463:   }

465:   ctx->GotDiag=0;
466:    PetscObjectReference((PetscObject)mat);

468:   ierr=VecGetLocalSize(d2,&nloc);
469:   ierr=VecGetSize(d2,&n);

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

473:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
474:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
475:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
476:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
477:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
478:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
479:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
480:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
481:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
482:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
483:   MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);
484:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
485:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
486:   MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);

488:   PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
489:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));

491:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
492:   return(0);
493: }