Actual source code: adamat.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <../src/tao/matrix/adamat.h>                /*I  "mat.h"  I*/

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

  8:    Collective on matrix

 10:    Input Parameters:
 11: +  mat - matrix of arbitrary type
 12: .  d1 - A vector with diagonal elements of D1
 13: -  d2 - A vector

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

 18:    Notes:
 19:    The user provides the input data and is responsible for destroying
 20:    this data after matrix J has been destroyed.
 21:    The operation MatMult(A,D2,D1) must be well defined.
 22:    Before calling the operation MatGetDiagonal(), the function
 23:    MatADAComputeDiagonal() must be called.  The matrices A and D1 must
 24:    be the same during calls to MatADAComputeDiagonal() and
 25:    MatGetDiagonal().

 27:    Level: developer

 29: .seealso: MatCreate()
 30: @*/
 31: PetscErrorCode MatCreateADA(Mat mat,Vec d1, Vec d2, Mat *J)
 32: {
 33:   MPI_Comm       comm=((PetscObject)mat)->comm;
 34:   TaoMatADACtx   ctx;
 36:   PetscInt       nloc,n;

 39:   PetscNew(&ctx);
 40:   ctx->A=mat;
 41:   ctx->D1=d1;
 42:   ctx->D2=d2;
 43:   if (d1){
 44:     VecDuplicate(d1,&ctx->W);
 45:     PetscObjectReference((PetscObject)d1);
 46:   } else {
 47:     ctx->W = NULL;
 48:   }
 49:   if (d2){
 50:     VecDuplicate(d2,&ctx->W2);
 51:     VecDuplicate(d2,&ctx->ADADiag);
 52:      PetscObjectReference((PetscObject)d2);
 53:   } else {
 54:     ctx->W2      = NULL;
 55:     ctx->ADADiag = NULL;
 56:   }

 58:   ctx->GotDiag=0;
 59:    PetscObjectReference((PetscObject)mat);

 61:   ierr=VecGetLocalSize(d2,&nloc);
 62:   ierr=VecGetSize(d2,&n);

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

 66:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_ADA);
 67:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_ADA);
 68:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_ADA);
 69:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_ADA);
 70:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_ADA);
 71:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_ADA);
 72:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_ADA);
 73:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_ADA);
 74:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_ADA);
 75:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_ADA);
 76:   MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_ADA);
 77:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_ADA);
 78:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_ADA);
 79:   MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_ADA);

 81:   PetscLogObjectParent((PetscObject)(*J),(PetscObject)ctx->W);
 82:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));

 84:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
 85:   return(0);
 86: }

 90: PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
 91: {
 92:   TaoMatADACtx   ctx;
 93:   PetscReal      one = 1.0;

 97:   MatShellGetContext(mat,(void **)&ctx);
 98:   MatMult(ctx->A,a,ctx->W);
 99:   if (ctx->D1){
100:     VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
101:   }
102:   MatMultTranspose(ctx->A,ctx->W,y);
103:   if (ctx->D2){
104:     VecPointwiseMult(ctx->W2, ctx->D2, a);
105:     VecAXPY(y, one, ctx->W2);
106:   }
107:   return(0);
108: }

112: PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
113: {

117:   MatMult_ADA(mat,a,y);
118:   return(0);
119: }

123: PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
124: {
125:   TaoMatADACtx   ctx;
126:   PetscReal      zero=0.0,one = 1.0;

130:   MatShellGetContext(M,(void **)&ctx);
131:   if (!ctx->D2){
132:     VecDuplicate(D,&ctx->D2);
133:     VecSet(ctx->D2, zero);
134:   }
135:   VecAXPY(ctx->D2, one, D);
136:   return(0);
137: }

141: PetscErrorCode MatDestroy_ADA(Mat mat)
142: {
144:   TaoMatADACtx   ctx;

147:   ierr=MatShellGetContext(mat,(void **)&ctx);
148:   ierr=VecDestroy(&ctx->W);
149:   ierr=VecDestroy(&ctx->W2);
150:   ierr=VecDestroy(&ctx->ADADiag);
151:   ierr=MatDestroy(&ctx->A);
152:   ierr=VecDestroy(&ctx->D1);
153:   ierr=VecDestroy(&ctx->D2);
154:   PetscFree(ctx);
155:   return(0);
156: }

160: PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
161: {
163:   return(0);
164: }

168: PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
169: {
171:   TaoMatADACtx   ctx;

174:   MatShellGetContext(Y,(void **)&ctx);
175:   VecShift(ctx->D2,a);
176:   return(0);
177: }

181: PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
182: {
183:   PetscErrorCode    ierr;
184:   TaoMatADACtx      ctx;
185:   Mat               A2;
186:   Vec               D1b=NULL,D2b;

189:   MatShellGetContext(mat,(void **)&ctx);
190:   MatDuplicate(ctx->A,op,&A2);
191:   if (ctx->D1){
192:     VecDuplicate(ctx->D1,&D1b);
193:     VecCopy(ctx->D1,D1b);
194:   }
195:   VecDuplicate(ctx->D2,&D2b);
196:   VecCopy(ctx->D2,D2b);
197:   MatCreateADA(A2,D1b,D2b,M);
198:   if (ctx->D1){
199:     PetscObjectDereference((PetscObject)D1b);
200:   }
201:   PetscObjectDereference((PetscObject)D2b);
202:   PetscObjectDereference((PetscObject)A2);
203:   return(0);
204: }

208: PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
209: {
211:   TaoMatADACtx   ctx1,ctx2;

214:   MatShellGetContext(A,(void **)&ctx1);
215:   MatShellGetContext(B,(void **)&ctx2);
216:   VecEqual(ctx1->D2,ctx2->D2,flg);
217:   if (*flg==PETSC_TRUE){
218:     VecEqual(ctx1->D1,ctx2->D1,flg);
219:   }
220:   if (*flg==PETSC_TRUE){
221:     MatEqual(ctx1->A,ctx2->A,flg);
222:   }
223:   return(0);
224: }

228: PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
229: {
231:   TaoMatADACtx   ctx;

234:   MatShellGetContext(mat,(void **)&ctx);
235:   VecScale(ctx->D1,a);
236:   if (ctx->D2){
237:     VecScale(ctx->D2,a);
238:   }
239:   return(0);
240: }

244: PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
245: {
247:   TaoMatADACtx   ctx;

250:   if (*B){
251:     MatShellGetContext(mat,(void **)&ctx);
252:     MatDuplicate(mat,MAT_COPY_VALUES,B);
253:   }
254:   return(0);
255: }

259: PetscErrorCode MatADAComputeDiagonal(Mat mat)
260: {
262:   PetscInt       i,m,n,low,high;
263:   PetscScalar    *dtemp,*dptr;
264:   TaoMatADACtx   ctx;

267:   MatShellGetContext(mat,(void **)&ctx);
268:   MatGetOwnershipRange(mat, &low, &high);
269:   MatGetSize(mat,&m,&n);

271:   PetscMalloc1(n,&dtemp );

273:   for (i=0; i<n; i++){
274:     MatGetColumnVector(ctx->A, ctx->W, i);
275:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
276:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
277:   }
278:   for (i=0; i<n; i++){
279:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
280:   }

282:   VecGetArray(ctx->ADADiag,&dptr);
283:   for (i=low; i<high; i++){
284:     dptr[i-low]= dtemp[i];
285:   }
286:   VecRestoreArray(ctx->ADADiag,&dptr);
287:   PetscFree(dtemp);
288:   return(0);
289: }

293: PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
294: {
295:   PetscErrorCode  ierr;
296:   PetscReal       one=1.0;
297:   TaoMatADACtx    ctx;

300:   MatShellGetContext(mat,(void **)&ctx);
301:   MatADAComputeDiagonal(mat);
302:   ierr=VecCopy(ctx->ADADiag,v);
303:   if (ctx->D2){
304:     ierr=VecAXPY(v, one, ctx->D2);
305:   }
306:   return(0);
307: }

311: PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
312: {
314:   PetscInt       i;

317:   if (scall == MAT_INITIAL_MATRIX) {
318:     PetscMalloc1(n+1,B );
319:   }
320:   for ( i=0; i<n; i++ ) {
321:     MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
322:   }
323:   return(0);
324: }

328: PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll, Mat *newmat)
329: {
330:   PetscErrorCode    ierr;
331:   PetscInt          low,high;
332:   IS                ISrow;
333:   Vec               D1,D2;
334:   Mat               Atemp;
335:   TaoMatADACtx      ctx;
336:   PetscBool         isequal;

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

343:   MatGetOwnershipRange(ctx->A,&low,&high);
344:   ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
345:   MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
346:   ISDestroy(&ISrow);

348:   if (ctx->D1){
349:     ierr=VecDuplicate(ctx->D1,&D1);
350:     ierr=VecCopy(ctx->D1,D1);
351:   } else {
352:     D1 = NULL;
353:   }

355:   if (ctx->D2){
356:     Vec D2sub;

358:     ierr=VecGetSubVector(ctx->D2,isrow,&D2sub);
359:     ierr=VecDuplicate(D2sub,&D2);
360:     ierr=VecCopy(D2sub,D2);
361:     ierr=VecRestoreSubVector(ctx->D2,isrow,&D2sub);
362:   } else {
363:     D2 = NULL;
364:   }

366:   MatCreateADA(Atemp,D1,D2,newmat);
367:   MatShellGetContext(*newmat,(void **)&ctx);
368:   PetscObjectDereference((PetscObject)Atemp);
369:   if (ctx->D1){
370:     PetscObjectDereference((PetscObject)D1);
371:   }
372:   if (ctx->D2){
373:     PetscObjectDereference((PetscObject)D2);
374:   }
375:   return(0);
376: }

380: PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
381: {
383:   PetscInt       m,n;

386:   MatGetSize(mat,&m,&n);

388:   if (*ncols>0){
389:     PetscMalloc1(*ncols,cols );
390:     PetscMalloc1(*ncols,vals );
391:   } else {
392:     *cols=NULL;
393:     *vals=NULL;
394:   }
395:   return(0);
396: }

400: PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
401: {

405:   if (*ncols>0){
406:     PetscFree(*cols);
407:     PetscFree(*vals);
408:   }
409:   *cols=NULL;
410:   *vals=NULL;
411:   return(0);
412: }

416: PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
417: {
419:   PetscInt       low,high;
420:   PetscScalar    zero=0.0,one=1.0;

423:   ierr=VecSet(Y, zero);
424:   ierr=VecGetOwnershipRange(Y,&low,&high);
425:   if (col>=low && col<high){
426:     ierr=VecSetValue(Y,col,one,INSERT_VALUES);
427:   }
428:   ierr=VecAssemblyBegin(Y);
429:   ierr=VecAssemblyEnd(Y);
430:   ierr=MatMult_ADA(mat,Y,Y);
431:   return(0);
432: }

434: PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
435: {
437:   PetscMPIInt    size;
438:   PetscBool      sametype, issame, isdense, isseqdense;
439:   TaoMatADACtx   ctx;

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

445:   PetscObjectTypeCompare((PetscObject)mat,newtype,&sametype);
446:   PetscObjectTypeCompare((PetscObject)mat,MATSAME,&issame);
447:   PetscObjectTypeCompare((PetscObject)mat,MATMPIDENSE,&isdense);
448:   PetscObjectTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);

450:   if (sametype || issame) {
451:     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);
452:   } else if (isdense) {
453:     PetscInt          i,j,low,high,m,n,M,N;
454:     const PetscScalar *dptr;
455:     Vec               X;

457:     VecDuplicate(ctx->D2,&X);
458:     ierr=MatGetSize(mat,&M,&N);
459:     ierr=MatGetLocalSize(mat,&m,&n);
460:     MatCreateDense(((PetscObject)mat)->comm,m,m,N,N,NULL,NewMat);
461:     MatGetOwnershipRange(*NewMat,&low,&high);
462:     for (i=0;i<M;i++){
463:       MatGetColumnVector_ADA(mat,X,i);
464:       VecGetArrayRead(X,&dptr);
465:       for (j=0; j<high-low; j++){
466:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
467:       }
468:       VecRestoreArrayRead(X,&dptr);
469:     }
470:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
471:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
472:     VecDestroy(&X);
473:   } else if (isseqdense && size==1){
474:     PetscInt          i,j,low,high,m,n,M,N;
475:     const PetscScalar *dptr;
476:     Vec               X;

478:     VecDuplicate(ctx->D2,&X);
479:     MatGetSize(mat,&M,&N);
480:     MatGetLocalSize(mat,&m,&n);
481:     MatCreateSeqDense(((PetscObject)mat)->comm,N,N,NULL,NewMat);
482:     MatGetOwnershipRange(*NewMat,&low,&high);
483:     for (i=0;i<M;i++){
484:       MatGetColumnVector_ADA(mat,X,i);
485:       VecGetArrayRead(X,&dptr);
486:       for (j=0; j<high-low; j++){
487:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
488:       }
489:       VecRestoreArrayRead(X,&dptr);
490:     }
491:     MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
492:     MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
493:     VecDestroy(&X);
494:   } else SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
495:   return(0);
496: }

500: PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
501: {
503:   TaoMatADACtx   ctx;

506:   MatShellGetContext(mat,(void **)&ctx);
507:   if (type == NORM_FROBENIUS) {
508:     *norm = 1.0;
509:   } else if (type == NORM_1 || type == NORM_INFINITY) {
510:     *norm = 1.0;
511:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
512:   return(0);
513: }