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=mat->comm;
     34:   TaoMatADACtx ctx;
     35:   int          info,nloc,n;

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

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

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

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

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

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

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

     90:   info = MatSetOption(*J,MAT_SYMMETRIC);CHKERRQ(info);
     91:   return(0);
     92: }

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

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

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

    119: int MatMultTranspose_ADA(Mat mat,Vec a,Vec y)
    120: {
    121:   int info;

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

    130: int MatDiagonalShift_ADA(Vec D, Mat M)
    131: {
    132:   TaoMatADACtx ctx;
    133:   PetscScalar        zero=0.0,one = 1.0;
    134:   int           info;

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

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

    145:   return(0);
    146: }

    150: int MatDestroy_ADA(Mat mat)
    151: {
    152:   int          info;
    153:   TaoMatADACtx ctx;

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

    169: int MatView_ADA(Mat mat,PetscViewer viewer)
    170: {

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

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

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

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

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

    228:   return(0);
    229: }

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

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

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

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

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

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

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

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

    293:   info = MatGetOwnershipRange(mat, &low, &high);CHKERRQ(info);
    294:   info = MatGetSize(mat,&m,&n);CHKERRQ(info);
    295:
    296:   info = PetscMalloc( n*sizeof(PetscScalar),&dtemp ); CHKERRQ(info);

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

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

    318: int MatGetDiagonal_ADA(Mat mat,Vec v)
    319: {
    320:   int          info;
    321:   PetscScalar       one=1.0;
    322:   TaoMatADACtx ctx;

    325:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
    326:   info = MatADAComputeDiagonal(mat);
    327:   info=VecCopy(ctx->ADADiag,v);CHKERRQ(info);
    328:   if (ctx->D2){
    329:     info=VecAXPY(v, one, ctx->D2);CHKERRQ(info);
    330:   }

    332:   return(0);
    333: }

    337: int MatGetSubMatrices_ADA(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
    338: {
    339:   int info,i;

    342:   if (scall == MAT_INITIAL_MATRIX) {
    343:     info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
    344:   }

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

    354: int MatGetSubMatrix_ADA(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,
    355:                         Mat *newmat)
    356: {
    357:   int          info,low,high;
    358:   int          n,nlocal,i;
    359:   int          *iptr;
    360:   PetscScalar       *dptr,*ddptr,zero=0.0;
    361:   VecType type_name;
    362:   IS           ISrow;
    363:   Vec          D1,D2;
    364:   Mat          Atemp;
    365:   TaoMatADACtx ctx;

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

    371:   info = MatGetOwnershipRange(ctx->A,&low,&high);CHKERRQ(info);
    372:   info = ISCreateStride(mat->comm,high-low,low,1,&ISrow);CHKERRQ(info);
    373:   info = MatGetSubMatrix(ctx->A,ISrow,iscol,csize,cll,&Atemp);CHKERRQ(info);
    374:   info = ISDestroy(ISrow);CHKERRQ(info);

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

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

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

    419: int MatGetRowADA(Mat mat,int row,int *ncols,int **cols,PetscScalar **vals)
    420: {
    421:   int info,m,n;

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

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

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

    453: int MatGetColumnVector_ADA(Mat mat,Vec Y, int col)
    454: {
    455:   int    info,low,high;
    456:   PetscScalar zero=0.0,one=1.0;

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

    468:   return(0);
    469: }

    471: int MatConvert_ADA(Mat mat,MatType newtype,Mat *NewMat)
    472: {
    473:   int info,size;
    474:   TaoMatADACtx  ctx;

    477:   info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
    478:   MPI_Comm_size(mat->comm,&size);

    480:   if (newtype==MATSAME){

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

    484:   } else if (newtype==MATMPIDENSE){

    486:     int i,j,low,high,m,n,M,N;
    487:     PetscScalar *dptr;
    488:     Vec X;

    490:     info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
    491:     info=MatGetSize(mat,&M,&N);CHKERRQ(info);
    492:     info=MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
    493:     info = MatCreateMPIDense(mat->comm,m,m,N,N,PETSC_NULL,NewMat);
    494:     CHKERRQ(info);
    495:     info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
    496:     for (i=0;i<M;i++){
    497:       info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
    498:       info = VecGetArray(X,&dptr);CHKERRQ(info);
    499:       for (j=0; j<high-low; j++){
    500:         info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
    501:       }
    502:       info=VecRestoreArray(X,&dptr);CHKERRQ(info);
    503:     }
    504:     info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
    505:     info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
    506:     info = VecDestroy(X);CHKERRQ(info);

    508:   } else if (newtype==MATSEQDENSE && size==1){

    510:     int i,j,low,high,m,n,M,N;
    511:     PetscScalar *dptr;
    512:     Vec X;

    514:     info = VecDuplicate(ctx->D2,&X);CHKERRQ(info);
    515:     info = MatGetSize(mat,&M,&N);CHKERRQ(info);
    516:     info = MatGetLocalSize(mat,&m,&n);CHKERRQ(info);
    517:     info = MatCreateSeqDense(mat->comm,N,N,PETSC_NULL,NewMat);
    518:     CHKERRQ(info);
    519:     info = MatGetOwnershipRange(*NewMat,&low,&high);CHKERRQ(info);
    520:     for (i=0;i<M;i++){
    521:       info = MatGetColumnVector_ADA(mat,X,i);CHKERRQ(info);
    522:       info = VecGetArray(X,&dptr);CHKERRQ(info);
    523:       for (j=0; j<high-low; j++){
    524:         info = MatSetValue(*NewMat,low+j,i,dptr[j],INSERT_VALUES);CHKERRQ(info);
    525:       }
    526:       info=VecRestoreArray(X,&dptr);CHKERRQ(info);
    527:     }
    528:     info=MatAssemblyBegin(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
    529:     info=MatAssemblyEnd(*NewMat,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
    530:     info=VecDestroy(X);CHKERRQ(info);

    532:   } else {
    533:     SETERRQ(1,"No support to convert objects to that type");
    534:   }
    535:   return(0);
    536: }

    540: int MatNorm_ADA(Mat mat,NormType type,PetscReal *norm)
    541: {
    542:   int info;
    543:   TaoMatADACtx  ctx;

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

    548:   if (type == NORM_FROBENIUS) {
    549:     *norm = 1.0;
    550:   } else if (type == NORM_1 || type == NORM_INFINITY) {
    551:     *norm = 1.0;
    552:   } else {
    553:     SETERRQ(PETSC_ERR_SUP,"No two norm");
    554:   }
    555:   return(0);
    556: }