Actual source code: normalmat.c

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

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

  8:    Collective on matrix

 10:    Input Parameters:
 11: +  mat - matrix of arbitrary type
 12: .  D1 - A vector
 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: int MatCreateADA(Mat mat,Vec D1, Vec D2, Mat *J)
 32: {
 33:   MPI_Comm     comm=((PetscObject)mat)->comm;
 34:   TaoMatADACtx ctx;
 35:   int          info;
 36:   PetscInt       nloc,n;

 39:   /*
 40:   info=MatCheckVecs(mat,D1,D2,&flg);CHKERRQ(info);
 41:   if (flg==PETSC_FALSE){
 42:     SETERRQ(PETSC_ERR_SUP,"InCompatible matrix and vector for ADA^T matrix");
 43:   }
 44:   */
 45:   info = PetscNew(_p_TaoMatADACtx,&ctx);CHKERRQ(info);

 47:   ctx->A=mat;
 48:   ctx->D1=D1;
 49:   ctx->D2=D2;
 50:   if (D1){
 51:     info = VecDuplicate(D1,&ctx->W);CHKERRQ(info);
 52:     info =  PetscObjectReference((PetscObject)D1);CHKERRQ(info);
 53:   } else {
 54:     ctx->W=0;
 55:   }
 56:   if (D2){
 57:     info = VecDuplicate(D2,&ctx->W2);CHKERRQ(info);
 58:     info = VecDuplicate(D2,&ctx->ADADiag);CHKERRQ(info);
 59:     info =  PetscObjectReference((PetscObject)D2);CHKERRQ(info);
 60:   } else {
 61:     ctx->W2=0;
 62:     ctx->ADADiag=0;
 63:   }

 65:   ctx->GotDiag=0;
 66:   info =  PetscObjectReference((PetscObject)mat);CHKERRQ(info);

 68:   info=VecGetLocalSize(D2,&nloc);CHKERRQ(info);
 69:   info=VecGetSize(D2,&n);CHKERRQ(info);

 71:   info = MatCreateShell(comm,nloc,nloc,n,n,ctx,J);CHKERRQ(info);

 73:   info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_ADA);CHKERRQ(info);
 74:   info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_ADA);CHKERRQ(info);
 75:   info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_ADA);CHKERRQ(info);
 76:   info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_ADA);CHKERRQ(info);
 77:   info = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)())MatDiagonalSet_ADA);CHKERRQ(info);
 78:   info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_ADA);CHKERRQ(info);
 79:   info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_ADA);CHKERRQ(info);
 80:   info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_ADA);CHKERRQ(info);
 81:   info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_ADA);CHKERRQ(info);
 82:   info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_ADA);CHKERRQ(info);
 83:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_ADA);CHKERRQ(info);
 84:   info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_ADA);CHKERRQ(info);
 85:   info = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)())MatDuplicate_ADA);CHKERRQ(info);
 86:   info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_ADA);CHKERRQ(info);

 88:   info = PetscLogObjectParent(*J,ctx->W); CHKERRQ(info);
 89:   info = PetscLogObjectParent(mat,*J); CHKERRQ(info);

 91:   info = MatSetOption(*J,MAT_SYMMETRIC, PETSC_TRUE);CHKERRQ(info);

 93:   return(0);  
 94: }

 98: int MatMult_ADA(Mat mat,Vec a,Vec y)
 99: {
100:   TaoMatADACtx ctx;
101:   PetscScalar        one = 1.0;
102:   int           info;

105:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

107:   info = MatMult(ctx->A,a,ctx->W);CHKERRQ(info);
108:   if (ctx->D1){
109:     info = VecPointwiseMult(ctx->W,ctx->D1,ctx->W);CHKERRQ(info);
110:   }
111:   info = MatMultTranspose(ctx->A,ctx->W,y);CHKERRQ(info);
112:   if (ctx->D2){
113:     info = VecPointwiseMult(ctx->W2, ctx->D2, a);CHKERRQ(info);
114:     info = VecAXPY(y, one, ctx->W2);CHKERRQ(info);
115:   }
116:   return(0);
117: } 

121: int MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
122: {
123:   int info;

126:   info = MatMult_ADA(mat,a,y);CHKERRQ(info);
127:   return(0);
128: } 

132: int MatDiagonalSet_ADA(Vec D, Mat M)
133: {
134:   TaoMatADACtx ctx;
135:   PetscScalar        zero=0.0,one = 1.0;
136:   int           info;

139:   info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);

141:   if (ctx->D2==PETSC_NULL){
142:     info = VecDuplicate(D,&ctx->D2);CHKERRQ(info);
143:     info = VecSet(ctx->D2, zero);CHKERRQ(info);
144:   }
145:   info = VecAXPY(ctx->D2, one, D);CHKERRQ(info);

147:   return(0);
148: } 

152: int MatDestroy_ADA(Mat mat)
153: {
154:   int          info;
155:   TaoMatADACtx ctx;

158:   info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
159:   info=VecDestroy(ctx->W);CHKERRQ(info);
160:   info=VecDestroy(ctx->W2);CHKERRQ(info);
161:   info=VecDestroy(ctx->ADADiag);CHKERRQ(info);
162:   info=MatDestroy(ctx->A);CHKERRQ(info);
163:   info=VecDestroy(ctx->D1);CHKERRQ(info);
164:   info=VecDestroy(ctx->D2);CHKERRQ(info);
165:   info = PetscFree(ctx); CHKERRQ(info);
166:   return(0);
167: }

171: int MatView_ADA(Mat mat,PetscViewer viewer)
172: {

175:   /*
176:   info = ViewerGetFormat(viewer,&format);CHKERRQ(info);
177:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_LONG) {
178:     return(0);  / * do nothing for now * /
179:   }
180:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
181:   info = MatView(ctx->A,viewer);CHKERRQ(info);
182:   if (ctx->D1){
183:     info = VecView(ctx->D1,viewer);CHKERRQ(info);
184:   }
185:   if (ctx->D2){
186:     info = VecView(ctx->D2,viewer);CHKERRQ(info);
187:   }
188:   */
189:   return(0);
190: }

194: int MatShift_ADA(Mat Y, PetscScalar a)
195: {
196:   int          info;
197:   TaoMatADACtx ctx;

200:   info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
201:   info = VecShift(ctx->D2,a);CHKERRQ(info);
202:   return(0);
203: }

207: int MatDuplicate_ADA(Mat mat,MatDuplicateOption op,Mat *M)
208: {
209:   int          info;
210:   TaoMatADACtx ctx;
211:   Mat          A2;
212:   Vec          D1b=NULL,D2b;

215:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
216:   info = MatDuplicate(ctx->A,op,&A2);CHKERRQ(info);
217:   if (ctx->D1){
218:     info = VecDuplicate(ctx->D1,&D1b);CHKERRQ(info);
219:     info = VecCopy(ctx->D1,D1b);CHKERRQ(info);
220:   }
221:   info = VecDuplicate(ctx->D2,&D2b);CHKERRQ(info);
222:   info = VecCopy(ctx->D2,D2b);CHKERRQ(info);
223:   info = MatCreateADA(A2,D1b,D2b,M);CHKERRQ(info);
224:   if (ctx->D1){
225:   info = PetscObjectDereference((PetscObject)D1b);CHKERRQ(info);
226:   }
227:   info = PetscObjectDereference((PetscObject)D2b);CHKERRQ(info);
228:   info = PetscObjectDereference((PetscObject)A2);CHKERRQ(info);

230:   return(0);
231: }

235: int MatEqual_ADA(Mat A,Mat B,PetscTruth *flg)
236: {
237:   int          info;
238:   TaoMatADACtx  ctx1,ctx2;

241:   info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
242:   info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
243:   info = VecEqual(ctx1->D2,ctx2->D2,flg);CHKERRQ(info);
244:   if (*flg==PETSC_TRUE){
245:     info = VecEqual(ctx1->D1,ctx2->D1,flg);CHKERRQ(info);
246:   }
247:   if (*flg==PETSC_TRUE){
248:     info = MatEqual(ctx1->A,ctx2->A,flg);CHKERRQ(info);
249:   }
250:   return(0);
251: }

255: int MatScale_ADA(Mat mat, PetscScalar a)
256: {
257:   int          info;
258:   TaoMatADACtx ctx;

261:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
262:   info = VecScale(ctx->D1,a);CHKERRQ(info);
263:   if (ctx->D2){
264:     info = VecScale(ctx->D2,a);CHKERRQ(info);
265:   }
266:   return(0);
267: }

271: int MatTranspose_ADA(Mat mat,Mat *B)
272: {
273:   int          info;
274:   TaoMatADACtx ctx;

277:   if (*B){
278:     info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
279:     info = MatDuplicate(mat,MAT_COPY_VALUES,B);CHKERRQ(info);
280:   }
281:   return(0);
282: }

286: int MatADAComputeDiagonal(Mat mat)
287: {
288:   int          info;
289:   PetscInt i,m,n,low,high;
290:   PetscScalar       *dtemp,*dptr;
291:   TaoMatADACtx ctx;

294:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

296:   info = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(info);
297:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);
298:   
299:   info = PetscMalloc( n*sizeof(PetscScalar),&dtemp ); CHKERRQ(info);

301:   for (i=0; i<n; i++){
302:     info = MatGetColumnVector(ctx->A, ctx->W, i);CHKERRQ(info);
303:     info = VecPointwiseMult(ctx->W,ctx->W,ctx->W);CHKERRQ(info);
304:     info = VecDotBegin(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
305:   }
306:   for (i=0; i<n; i++){
307:     info = VecDotEnd(ctx->D1, ctx->W,dtemp+i);CHKERRQ(info);
308:   } 

310:   info = VecGetArray(ctx->ADADiag,&dptr);CHKERRQ(info);
311:   for (i=low; i<high; i++){
312:     dptr[i-low]= dtemp[i];
313:   }
314:   info = VecRestoreArray(ctx->ADADiag,&dptr);CHKERRQ(info);
315:   if (dtemp) {
316:     info = PetscFree(dtemp); CHKERRQ(info);
317:   }
318:   return(0);
319: }

323: int MatGetDiagonal_ADA(Mat mat,Vec v)
324: {
325:   int          info;
326:   PetscScalar       one=1.0;
327:   TaoMatADACtx ctx;

330:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
331:   info = MatADAComputeDiagonal(mat);
332:   info=VecCopy(ctx->ADADiag,v);CHKERRQ(info);
333:   if (ctx->D2){
334:     info=VecAXPY(v, one, ctx->D2);CHKERRQ(info);
335:   }

337:   return(0);
338: }

342: int MatGetSubMatrices_ADA(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
343: {
344:   int info,i;

347:   if (scall == MAT_INITIAL_MATRIX) {
348:     info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
349:   }

351:   for ( i=0; i<n; i++ ) {
352:     info = MatGetSubMatrix_ADA(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(info);
353:   }
354:   return(0);
355: }

359: int MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,MatReuse cll,
360:                         Mat *newmat)
361: {
362:   int          info;
363:   PetscInt       i,low,high;
364:   PetscInt       n,nlocal;
365:   const PetscInt       *iptr;
366:   PetscScalar       *dptr,*ddptr,zero=0.0;
367:   const VecType type_name;
368:   IS           ISrow;
369:   Vec          D1,D2;
370:   Mat          Atemp;
371:   TaoMatADACtx ctx;


375:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

377:   info = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(info);
378:   info = ISCreateStride(((PetscObject)mat)->comm,high-low,low,1,&ISrow);CHKERRQ(info);
379:   info = MatGetSubMatrix(ctx->A,ISrow,iscol,cll,&Atemp);CHKERRQ(info);
380:   info = ISDestroy(ISrow);CHKERRQ(info);

382:   if (ctx->D1){
383:     info=VecDuplicate(ctx->D1,&D1);CHKERRQ(info);
384:     info=VecCopy(ctx->D1,D1);CHKERRQ(info);
385:   } else {
386:     D1=PETSC_NULL;
387:   }

389:   if (ctx->D2){
390:     info=ISGetSize(isrow,&n);CHKERRQ(info);
391:     info=ISGetLocalSize(isrow,&nlocal);CHKERRQ(info);
392:     info=VecCreate(((PetscObject)ctx->D2)->comm,&D2);CHKERRQ(info);
393:     info=VecGetType(ctx->D2,&type_name);CHKERRQ(info);
394:     info=VecSetSizes(D2,nlocal,n);CHKERRQ(info);
395:     info=VecSetType(D2,type_name);CHKERRQ(info);
396:     info=VecSet(D2, zero);CHKERRQ(info);
397:     info=VecGetArray(ctx->D2, &dptr); CHKERRQ(info);
398:     info=VecGetArray(D2, &ddptr); CHKERRQ(info);
399:     info=ISGetIndices(isrow,&iptr); CHKERRQ(info);
400:     for (i=0;i<nlocal;i++){
401:       ddptr[i] = dptr[iptr[i]-low];
402:     }
403:     info=ISRestoreIndices(isrow,&iptr); CHKERRQ(info);
404:     info=VecRestoreArray(D2, &ddptr); CHKERRQ(info);
405:     info=VecRestoreArray(ctx->D2, &dptr); CHKERRQ(info);
406:    
407:   } else {
408:     D2=PETSC_NULL;
409:   }

411:   info = MatCreateADA(Atemp,D1,D2,newmat);CHKERRQ(info);
412:   info = MatShellGetContext(*newmat,(void **)&ctx);CHKERRQ(info);
413:   info = PetscObjectDereference((PetscObject)Atemp);CHKERRQ(info);
414:   if (ctx->D1){
415:     info = PetscObjectDereference((PetscObject)D1);CHKERRQ(info);
416:   }
417:   if (ctx->D2){
418:     info = PetscObjectDereference((PetscObject)D2);CHKERRQ(info);
419:   }
420:   return(0);
421: }

425: int MatGetRowADA(Mat mat,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals)
426: {
427:   int info;
428:   PetscInt m,n;

431:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);

433:   if (*ncols>0){
434:     info = PetscMalloc( (*ncols)*sizeof(int),cols );CHKERRQ(info);
435:     info = PetscMalloc( (*ncols)*sizeof(PetscScalar),vals );CHKERRQ(info);
436:   } else {
437:     *cols=PETSC_NULL;
438:     *vals=PETSC_NULL;
439:   }
440:   
441:   return(0);
442: }

446: int MatRestoreRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
447: {
448:   int info;
450:   if (*ncols>0){
451:     info = PetscFree(*cols);  CHKERRQ(info);
452:     info = PetscFree(*vals);  CHKERRQ(info);
453:   }
454:   *cols=PETSC_NULL;
455:   *vals=PETSC_NULL;
456:   return(0);
457: }

461: int MatGetColumnVector_ADA(Mat mat,Vec Y, PetscInt col)
462: {
463:   int    info;
464:   PetscInt low,high;
465:   PetscScalar zero=0.0,one=1.0;

468:   info=VecSet(Y, zero);CHKERRQ(info);
469:   info=VecGetOwnershipRange(Y,&low,&high);CHKERRQ(info);
470:   if (col>=low && col<high){
471:     info=VecSetValue(Y,col,one,INSERT_VALUES);CHKERRQ(info);
472:   }
473:   info=VecAssemblyBegin(Y);CHKERRQ(info);
474:   info=VecAssemblyEnd(Y);CHKERRQ(info);
475:   info=MatMult_ADA(mat,Y,Y);CHKERRQ(info);

477:   return(0);
478: }

480: int MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
481: {
482:   int info,size;
483:   TaoMatADACtx  ctx;
484:   PetscTruth sametype, issame, ismpidense, isseqdense;
486:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
487:   MPI_Comm_size(((PetscObject)mat)->comm,&size);

489:   info = PetscTypeCompare((PetscObject)mat,newtype,&sametype);CHKERRQ(info);
490:   info = PetscTypeCompare((PetscObject)mat,MATSAME,&issame); CHKERRQ(info);
491:   info = PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense); CHKERRQ(info);
492:   info = PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense); CHKERRQ(info);


495:   if (sametype || issame) {

497:     info=MatDuplicate(mat,MAT_COPY_VALUES,NewMat);CHKERRQ(info);

499:   } else if (ismpidense) {

501:     PetscInt i,j,low,high,m,n,M,N;
502:     PetscScalar *dptr;
503:     Vec X;

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

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

525:     PetscInt i,j,low,high,m,n,M,N;
526:     PetscScalar *dptr;
527:     Vec X;

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

547:   } else {
548:     SETERRQ(1,"No support to convert objects to that type");
549:   }
550:   return(0);
551: }

555: int MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
556: {
557:   int info;
558:   TaoMatADACtx  ctx;

561:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);

563:   if (type == NORM_FROBENIUS) {
564:     *norm = 1.0;
565:   } else if (type == NORM_1 || type == NORM_INFINITY) {
566:     *norm = 1.0;
567:   } else {
568:     SETERRQ(PETSC_ERR_SUP,"No two norm");
569:   }
570:   return(0);
571: }