Actual source code: adamat.c

  1: #include "adamat.h"                /*I  "mat.h"  I*/

  5: /*@C
  6:    MatCreateADA - Creates a matrix M=A^T D1 A + D2 where D1, D2 are diagonal
  7:    
  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;
 35:   PetscErrorCode          ierr;
 36:   PetscInt nloc,n;


 40:   PetscNew(_p_TaoMatADACtx,&ctx);

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

 60:   ctx->GotDiag=0;
 61:    PetscObjectReference((PetscObject)mat);

 63:   ierr=VecGetLocalSize(d2,&nloc);
 64:   ierr=VecGetSize(d2,&n);

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

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

 83:   PetscLogObjectParent(*J,ctx->W); 
 84:   PetscLogObjectParent(mat,*J); 

 86:   MatSetOption(*J,MAT_SYMMETRIC,PETSC_TRUE);
 87:   return(0);  
 88: }

 92: PetscErrorCode MatMult_ADA(Mat mat,Vec a,Vec y)
 93: {
 94:   TaoMatADACtx ctx;
 95:   PetscReal        one = 1.0;
 96:   PetscErrorCode           ierr;

 99:   MatShellGetContext(mat,(void **)&ctx);

101:   MatMult(ctx->A,a,ctx->W);
102:   if (ctx->D1){
103:     VecPointwiseMult(ctx->W,ctx->D1,ctx->W);
104:   }
105:   MatMultTranspose(ctx->A,ctx->W,y);
106:   if (ctx->D2){
107:     VecPointwiseMult(ctx->W2, ctx->D2, a);
108:     VecAXPY(y, one, ctx->W2);
109:   }
110:   return(0);
111: } 

115: PetscErrorCode MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
116: {

120:   MatMult_ADA(mat,a,y);
121:   return(0);
122: } 

126: PetscErrorCode MatDiagonalSet_ADA(Vec D, Mat M)
127: {
128:   TaoMatADACtx ctx;
129:   PetscReal        zero=0.0,one = 1.0;
130:   PetscErrorCode   ierr;

133:   MatShellGetContext(M,(void **)&ctx);

135:   if (ctx->D2==PETSC_NULL){
136:     VecDuplicate(D,&ctx->D2);
137:     VecSet(ctx->D2, zero);
138:   }
139:   VecAXPY(ctx->D2, one, D);

141:   return(0);
142: } 

146: PetscErrorCode MatDestroy_ADA(Mat mat)
147: {
148:   PetscErrorCode          ierr;
149:   TaoMatADACtx ctx;

152:   ierr=MatShellGetContext(mat,(void **)&ctx);
153:   if (ctx->W) {
154:     ierr=VecDestroy(&ctx->W);
155:   }
156:   if (ctx->W2) {
157:     ierr=VecDestroy(&ctx->W2);
158:   }
159:   if (ctx->ADADiag) {
160:     ierr=VecDestroy(&ctx->ADADiag);
161:   }
162:   if (ctx->A) {
163:     ierr=MatDestroy(&ctx->A);
164:   }
165:   if (ctx->D1) {
166:     ierr=VecDestroy(&ctx->D1);
167:   }
168:   if (ctx->D2) {
169:     ierr=VecDestroy(&ctx->D2);
170:   }
171:   PetscFree(ctx); 
172:   return(0);
173: }

177: PetscErrorCode MatView_ADA(Mat mat,PetscViewer viewer)
178: {

181:   return(0);
182: }

186: PetscErrorCode MatShift_ADA(Mat Y, PetscReal a)
187: {
189:   TaoMatADACtx   ctx;

192:   MatShellGetContext(Y,(void **)&ctx);
193:   VecShift(ctx->D2,a);
194:   return(0);
195: }

199: PetscErrorCode MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
200: {
201:   PetscErrorCode    ierr;
202:   TaoMatADACtx      ctx;
203:   Mat               A2;
204:   Vec               D1b=NULL,D2b;

207:   MatShellGetContext(mat,(void **)&ctx);
208:   MatDuplicate(ctx->A,op,&A2);
209:   if (ctx->D1){
210:     VecDuplicate(ctx->D1,&D1b);
211:     VecCopy(ctx->D1,D1b);
212:   }
213:   VecDuplicate(ctx->D2,&D2b);
214:   VecCopy(ctx->D2,D2b);
215:   MatCreateADA(A2,D1b,D2b,M);
216:   if (ctx->D1){
217:     PetscObjectDereference((PetscObject)D1b);
218:   }
219:   PetscObjectDereference((PetscObject)D2b);
220:   PetscObjectDereference((PetscObject)A2);

222:   return(0);
223: }

227: PetscErrorCode MatEqual_ADA(Mat A,Mat B,PetscBool *flg)
228: {
229:   PetscErrorCode    ierr;
230:   TaoMatADACtx  ctx1,ctx2;

233:   MatShellGetContext(A,(void **)&ctx1);
234:   MatShellGetContext(B,(void **)&ctx2);
235:   VecEqual(ctx1->D2,ctx2->D2,flg);
236:   if (*flg==PETSC_TRUE){
237:     VecEqual(ctx1->D1,ctx2->D1,flg);
238:   }
239:   if (*flg==PETSC_TRUE){
240:     MatEqual(ctx1->A,ctx2->A,flg);
241:   }
242:   return(0);
243: }

247: PetscErrorCode MatScale_ADA(Mat mat, PetscReal a)
248: {
249:   PetscErrorCode   ierr;
250:   TaoMatADACtx ctx;

253:   MatShellGetContext(mat,(void **)&ctx);
254:   VecScale(ctx->D1,a);
255:   if (ctx->D2){
256:     VecScale(ctx->D2,a);
257:   }
258:   return(0);
259: }

263: PetscErrorCode MatTranspose_ADA(Mat mat,Mat *B)
264: {
266:   TaoMatADACtx ctx;

269:   if (*B){
270:     MatShellGetContext(mat,(void **)&ctx);
271:     MatDuplicate(mat,MAT_COPY_VALUES,B);
272:   }
273:   return(0);
274: }

278: PetscErrorCode MatADAComputeDiagonal(Mat mat)
279: {
281:   PetscInt          i,m,n,low,high;
282:   PetscReal       *dtemp,*dptr;
283:   TaoMatADACtx ctx;

286:   MatShellGetContext(mat,(void **)&ctx);

288:   MatGetOwnershipRange(mat, &low, &high);
289:   MatGetSize(mat,&m,&n);
290:   
291:   PetscMalloc( n*sizeof(PetscReal),&dtemp ); 

293:   for (i=0; i<n; i++){
294:     MatGetColumnVector(ctx->A, ctx->W, i);
295:     VecPointwiseMult(ctx->W,ctx->W,ctx->W);
296:     VecDotBegin(ctx->D1, ctx->W,dtemp+i);
297:   }
298:   for (i=0; i<n; i++){
299:     VecDotEnd(ctx->D1, ctx->W,dtemp+i);
300:   } 

302:   VecGetArray(ctx->ADADiag,&dptr);
303:   for (i=low; i<high; i++){
304:     dptr[i-low]= dtemp[i];
305:   }
306:   VecRestoreArray(ctx->ADADiag,&dptr);
307:   if (dtemp) {
308:     PetscFree(dtemp); 
309:   }
310:   return(0);
311: }

315: PetscErrorCode MatGetDiagonal_ADA(Mat mat,Vec v)
316: {
317:   PetscErrorCode      ierr;
318:   PetscReal       one=1.0;
319:   TaoMatADACtx ctx;

322:   MatShellGetContext(mat,(void **)&ctx);
323:   MatADAComputeDiagonal(mat); 
324:   ierr=VecCopy(ctx->ADADiag,v);
325:   if (ctx->D2){
326:     ierr=VecAXPY(v, one, ctx->D2);
327:   }

329:   return(0);
330: }

334: PetscErrorCode MatGetSubMatrices_ADA(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
335: {
337:   PetscInt i;

340:   if (scall == MAT_INITIAL_MATRIX) {
341:     PetscMalloc( (n+1)*sizeof(Mat),B );
342:   }

344:   for ( i=0; i<n; i++ ) {
345:     MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);
346:   }
347:   return(0);
348: }

352: PetscErrorCode MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll,
353:                         Mat *newmat)
354: {
355:   PetscErrorCode          ierr;
356:   PetscInt low,high;
357:   PetscInt          n,nlocal,i;
358:   const PetscInt          *iptr;
359:   PetscReal       *dptr,*ddptr,zero=0.0;
360:   const VecType type_name;
361:   IS           ISrow;
362:   Vec          D1,D2;
363:   Mat          Atemp;
364:   TaoMatADACtx ctx;


368:   MatShellGetContext(mat,(void **)&ctx);

370:   MatGetOwnershipRange(ctx->A,&low,&high);
371:   ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);
372:   MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);
373:   ISDestroy(&ISrow);

375:   if (ctx->D1){
376:     ierr=VecDuplicate(ctx->D1,&D1);
377:     ierr=VecCopy(ctx->D1,D1);
378:   } else {
379:     D1=PETSC_NULL;
380:   }

382:   if (ctx->D2){
383:     ierr=ISGetSize(isrow,&n);
384:     ierr=ISGetLocalSize(isrow,&nlocal);
385:     ierr=VecCreate(((PetscObject)(ctx->D2))->comm,&D2);
386:     ierr=VecGetType(ctx->D2,&type_name);
387:     ierr=VecSetSizes(D2,nlocal,n);
388:     ierr=VecSetType(D2,type_name);
389:     ierr=VecSet(D2, zero);
390:     ierr=VecGetArray(ctx->D2, &dptr); 
391:     ierr=VecGetArray(D2, &ddptr); 
392:     ierr=ISGetIndices(isrow,&iptr); 
393:     for (i=0;i<nlocal;i++){
394:       ddptr[i] = dptr[iptr[i]-low];
395:     }
396:     ierr=ISRestoreIndices(isrow,&iptr); 
397:     ierr=VecRestoreArray(D2, &ddptr); 
398:     ierr=VecRestoreArray(ctx->D2, &dptr); 
399:    
400:   } else {
401:     D2=PETSC_NULL;
402:   }

404:   MatCreateADA(Atemp,D1,D2,newmat);
405:   MatShellGetContext(*newmat,(void **)&ctx);
406:   PetscObjectDereference((PetscObject)Atemp);
407:   if (ctx->D1){
408:     PetscObjectDereference((PetscObject)D1);
409:   }
410:   if (ctx->D2){
411:     PetscObjectDereference((PetscObject)D2);
412:   }
413:   return(0);
414: }

418: PetscErrorCode MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
419: {
421:   PetscInt m,n;

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

426:   if (*ncols>0){
427:     PetscMalloc( (*ncols)*sizeof(PetscInt),cols );
428:     PetscMalloc( (*ncols)*sizeof(PetscReal),vals );
429:   } else {
430:     *cols=PETSC_NULL;
431:     *vals=PETSC_NULL;
432:   }
433:   
434:   return(0);
435: }

439: PetscErrorCode MatRestoreRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscReal **vals)
440: {
443:   if (*ncols>0){
444:     PetscFree(*cols);  
445:     PetscFree(*vals);  
446:   }
447:   *cols=PETSC_NULL;
448:   *vals=PETSC_NULL;
449:   return(0);
450: }

454: PetscErrorCode MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
455: {
456:   PetscErrorCode    ierr;
457:   PetscInt low,high;
458:   PetscReal zero=0.0,one=1.0;

461:   ierr=VecSet(Y, zero);
462:   ierr=VecGetOwnershipRange(Y,&low,&high);
463:   if (col>=low && col<high){
464:     ierr=VecSetValue(Y,col,one,INSERT_VALUES);
465:   }
466:   ierr=VecAssemblyBegin(Y);
467:   ierr=VecAssemblyEnd(Y);
468:   ierr=MatMult_ADA(mat,Y,Y);

470:   return(0);
471: }

473: PetscErrorCode MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
474: {
476:   PetscMPIInt size;
477:   PetscBool sametype, issame, ismpidense, isseqdense;
478:   TaoMatADACtx  ctx;

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


485:   PetscTypeCompare((PetscObject)mat,newtype,&sametype);
486:   PetscTypeCompare((PetscObject)mat,MATSAME,&issame); 
487:   PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense); 
488:   PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense); 

490:   if (sametype || issame) {

492:     ierr=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);

494:   } else if (ismpidense) {

496:     PetscInt i,j,low,high,m,n,M,N;
497:     PetscReal *dptr;
498:     Vec X;

500:     VecDuplicate(ctx->D2,&X);
501:     ierr=MatGetSize(mat,&M,&N);
502:     ierr=MatGetLocalSize(mat,&m,&n);
503:     MatCreateMPIDense(((PetscObject)mat)->comm,m,m,N,N,PETSC_NULL,NewMat);
504:     
505:     MatGetOwnershipRange(*NewMat,&low,&high);
506:     for (i=0;i<M;i++){
507:       MatGetColumnVector_ADA(mat,X,i);
508:       VecGetArray(X,&dptr);
509:       for (j=0; j<high-low; j++){
510:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
511:       }
512:       ierr=VecRestoreArray(X,&dptr);
513:     }
514:     ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
515:     ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
516:     VecDestroy(&X);

518:   } else if (isseqdense && size==1){

520:     PetscInt i,j,low,high,m,n,M,N;
521:     PetscReal *dptr;
522:     Vec X;

524:     VecDuplicate(ctx->D2,&X);
525:     MatGetSize(mat,&M,&N);
526:     MatGetLocalSize(mat,&m,&n);
527:     MatCreateSeqDense(((PetscObject)mat)->comm,N,N,PETSC_NULL,NewMat);
528:     
529:     MatGetOwnershipRange(*NewMat,&low,&high);
530:     for (i=0;i<M;i++){
531:       MatGetColumnVector_ADA(mat,X,i);
532:       VecGetArray(X,&dptr);
533:       for (j=0; j<high-low; j++){
534:         MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);
535:       }
536:       ierr=VecRestoreArray(X,&dptr);
537:     }
538:     ierr=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);
539:     ierr=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);
540:     ierr=VecDestroy(&X);

542:   } else {
543:     SETERRQ(PETSC_COMM_SELF,1,"No support to convert objects to that type");
544:   }
545:   return(0);
546: }

550: PetscErrorCode MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
551: {
553:   TaoMatADACtx  ctx;

556:   MatShellGetContext(mat,(void **)&ctx);

558:   if (type == NORM_FROBENIUS) {
559:     *norm = 1.0;
560:   } else if (type == NORM_1 || type == NORM_INFINITY) {
561:     *norm = 1.0;
562:   } else {
563:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
564:   }
565:   return(0);
566: }