Actual source code: pchara.cu

      1: #include <petsc/private/pcimpl.h>
      2: #include <petsc/private/matimpl.h>

      4: typedef struct {
      5:   Mat         M;
      6:   PetscScalar s0;

      8:   /* sampler for Newton-Schultz */
      9:   Mat      S;
     10:   PetscInt hyperorder;
     11:   Vec      wns[4];
     12:   Mat      wnsmat[4];

     14:   /* convergence testing */
     15:   Mat T;
     16:   Vec w;

     18:   /* Support for PCSetCoordinates */
     19:   PetscInt  sdim;
     20:   PetscInt  nlocc;
     21:   PetscReal *coords;

     23:   /* Newton-Schultz customization */
     24:   PetscInt  maxits;
     25:   PetscReal rtol,atol;
     26:   PetscBool monitor;
     27:   PetscBool useapproximatenorms;
     28: } PC_HARA;

     30: static PetscErrorCode PCReset_HARA(PC pc)
     31: {
     32:   PC_HARA        *pchara = (PC_HARA*)pc->data;

     36:   pchara->sdim  = 0;
     37:   pchara->nlocc = 0;
     38:   PetscFree(pchara->coords);
     39:   MatDestroy(&pchara->M);
     40:   MatDestroy(&pchara->T);
     41:   VecDestroy(&pchara->w);
     42:   MatDestroy(&pchara->S);
     43:   VecDestroy(&pchara->wns[0]);
     44:   VecDestroy(&pchara->wns[1]);
     45:   VecDestroy(&pchara->wns[2]);
     46:   VecDestroy(&pchara->wns[3]);
     47:   MatDestroy(&pchara->wnsmat[0]);
     48:   MatDestroy(&pchara->wnsmat[1]);
     49:   MatDestroy(&pchara->wnsmat[2]);
     50:   MatDestroy(&pchara->wnsmat[3]);
     51:   return(0);
     52: }

     54: static PetscErrorCode PCSetCoordinates_HARA(PC pc, PetscInt sdim, PetscInt nlocc, PetscReal *coords)
     55: {
     56:   PC_HARA        *pchara = (PC_HARA*)pc->data;
     57:   PetscBool      reset = PETSC_TRUE;

     61:   if (pchara->sdim && sdim == pchara->sdim && nlocc == pchara->nlocc) {
     62:     PetscArraycmp(pchara->coords,coords,sdim*nlocc,&reset);
     63:     reset = (PetscBool)!reset;
     64:   }
     65:   MPIU_Allreduce(MPI_IN_PLACE,&reset,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)pc));
     66:   if (reset) {
     67:     PCReset_HARA(pc);
     68:     PetscMalloc1(sdim*nlocc,&pchara->coords);
     69:     PetscArraycpy(pchara->coords,coords,sdim*nlocc);
     70:     pchara->sdim  = sdim;
     71:     pchara->nlocc = nlocc;
     72:   }
     73:   return(0);
     74: }

     76: static PetscErrorCode PCDestroy_HARA(PC pc)
     77: {

     81:   PCReset_HARA(pc);
     82:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
     83:   PetscFree(pc->data);
     84:   return(0);
     85: }

     87: static PetscErrorCode PCSetFromOptions_HARA(PetscOptionItems *PetscOptionsObject,PC pc)
     88: {
     89:   PC_HARA       *pchara = (PC_HARA*)pc->data;

     93:   PetscOptionsHead(PetscOptionsObject,"HARA options");
     94:   PetscOptionsInt("-pc_hara_maxits","Maximum number of iterations for Newton-Schultz",NULL,pchara->maxits,&pchara->maxits,NULL);
     95:   PetscOptionsBool("-pc_hara_monitor","Monitor Newton-Schultz convergence",NULL,pchara->monitor,&pchara->monitor,NULL);
     96:   PetscOptionsReal("-pc_hara_atol","Absolute tolerance",NULL,pchara->atol,&pchara->atol,NULL);
     97:   PetscOptionsReal("-pc_hara_rtol","Relative tolerance",NULL,pchara->rtol,&pchara->rtol,NULL);
     98:   PetscOptionsInt("-pc_hara_hyperorder","Hyper power order of sampling",NULL,pchara->hyperorder,&pchara->hyperorder,NULL);
     99:   PetscOptionsTail();
    100:   return(0);
    101: }

    103: static PetscErrorCode PCApplyKernel_HARA(PC pc, Vec x, Vec y, PetscBool t)
    104: {
    105:   PC_HARA        *pchara = (PC_HARA*)pc->data;
    107:   PetscBool      flg = PETSC_FALSE;

    110:   MatAssembled(pchara->M,&flg);
    111:   if (flg) {
    112:     if (t) {
    113:       MatMultTranspose(pchara->M,x,y);
    114:     } else {
    115:       MatMult(pchara->M,x,y);
    116:     }
    117:   } else { /* Not assembled, initial approximation */
    118:     Mat A = pc->useAmat ? pc->mat : pc->pmat;

    120:     if (pchara->s0 < 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong scaling");
    121:     /* X_0 = s0 * A^T */
    122:     if (t) {
    123:       MatMult(A,x,y);
    124:     } else {
    125:       MatMultTranspose(A,x,y);
    126:     }
    127:     VecScale(y,pchara->s0);
    128:   }
    129:   return(0);
    130: }

    132: static PetscErrorCode PCApplyMatKernel_HARA(PC pc, Mat X, Mat Y, PetscBool t)
    133: {
    134:   PC_HARA        *pchara = (PC_HARA*)pc->data;
    136:   PetscBool      flg = PETSC_FALSE;

    139:   MatAssembled(pchara->M,&flg);
    140:   if (flg) {
    141:     if (t) {
    142:       MatTransposeMatMult(pchara->M,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Y);
    143:     } else {
    144:       MatMatMult(pchara->M,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Y);
    145:     }
    146:   } else { /* Not assembled, initial approximation */
    147:     Mat A = pc->useAmat ? pc->mat : pc->pmat;

    149:     if (pchara->s0 < 0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong scaling");
    150:     /* X_0 = s0 * A^T */
    151:     if (t) {
    152:       MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Y);
    153:     } else {
    154:       MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Y);
    155:     }
    156:     MatScale(Y,pchara->s0);
    157:   }
    158:   return(0);
    159: }

    161: static PetscErrorCode PCApplyMat_HARA(PC pc, Mat X, Mat Y)
    162: {

    166:   PCApplyMatKernel_HARA(pc,X,Y,PETSC_FALSE);
    167:   return(0);
    168: }

    170: static PetscErrorCode PCApplyTransposeMat_HARA(PC pc, Mat X, Mat Y)
    171: {

    175:   PCApplyMatKernel_HARA(pc,X,Y,PETSC_TRUE);
    176:   return(0);
    177: }

    179: static PetscErrorCode PCApply_HARA(PC pc, Vec x, Vec y)
    180: {

    184:   PCApplyKernel_HARA(pc,x,y,PETSC_FALSE);
    185:   return(0);
    186: }

    188: static PetscErrorCode PCApplyTranspose_HARA(PC pc, Vec x, Vec y)
    189: {

    193:   PCApplyKernel_HARA(pc,x,y,PETSC_TRUE);
    194:   return(0);
    195: }

    197: /* used to test norm of (M^-1 A - I) */
    198: static PetscErrorCode MatMultKernel_MAmI(Mat M, Vec x, Vec y, PetscBool t)
    199: {
    200:   PC             pc;
    201:   Mat            A;
    202:   PC_HARA        *pchara;

    206:   MatShellGetContext(M,(void*)&pc);
    207:   pchara = (PC_HARA*)pc->data;
    208:   if (!pchara->w) {
    209:     MatCreateVecs(pchara->M,&pchara->w,NULL);
    210:   }
    211:   A = pc->useAmat ? pc->mat : pc->pmat;
    212:   if (t) {
    213:     PCApplyTranspose_HARA(pc,x,pchara->w);
    214:     MatMultTranspose(A,pchara->w,y);
    215:   } else {
    216:     MatMult(A,x,pchara->w);
    217:     PCApply_HARA(pc,pchara->w,y);
    218:   }
    219:   VecAXPY(y,-1.0,x);
    220:   return(0);
    221: }

    223: static PetscErrorCode MatMult_MAmI(Mat A, Vec x, Vec y)
    224: {

    228:   MatMultKernel_MAmI(A,x,y,PETSC_FALSE);
    229:   return(0);
    230: }

    232: static PetscErrorCode MatMultTranspose_MAmI(Mat A, Vec x, Vec y)
    233: {

    237:   MatMultKernel_MAmI(A,x,y,PETSC_TRUE);
    238:   return(0);
    239: }

    241: /* HyperPower kernel:
    242: Y = R = x
    243: for i = 1 . . . l - 1 do
    244:   R = (I - AXk)R
    245:   Y = Y + R
    246: Y = XkY
    247: */
    248: static PetscErrorCode MatMultKernel_Hyper(Mat M, Vec x, Vec y, PetscBool t)
    249: {
    250:   PC             pc;
    251:   Mat            A;
    252:   PC_HARA        *pchara;
    253:   PetscInt       i;

    257:   MatShellGetContext(M,(void*)&pc);
    258:   A = pc->useAmat ? pc->mat : pc->pmat;
    259:   pchara = (PC_HARA*)pc->data;
    260:   MatCreateVecs(pchara->M,pchara->wns[0] ? NULL : &pchara->wns[0],pchara->wns[1] ? NULL : &pchara->wns[1]);
    261:   MatCreateVecs(pchara->M,pchara->wns[2] ? NULL : &pchara->wns[2],pchara->wns[3] ? NULL : &pchara->wns[3]);
    262:   VecCopy(x,pchara->wns[0]);
    263:   VecCopy(x,pchara->wns[3]);
    264:   if (t) {
    265:     for (i=0;i<pchara->hyperorder-1;i++) {
    266:       MatMultTranspose(A,pchara->wns[0],pchara->wns[1]);
    267:       PCApplyTranspose_HARA(pc,pchara->wns[1],pchara->wns[2]);
    268:       VecAXPBYPCZ(pchara->wns[3],-1.,1.,1.,pchara->wns[2],pchara->wns[0]);
    269:     }
    270:     PCApplyTranspose_HARA(pc,pchara->wns[3],y);
    271:   } else {
    272:     for (i=0;i<pchara->hyperorder-1;i++) {
    273:       PCApply_HARA(pc,pchara->wns[0],pchara->wns[1]);
    274:       MatMult(A,pchara->wns[1],pchara->wns[2]);
    275:       VecAXPBYPCZ(pchara->wns[3],-1.,1.,1.,pchara->wns[2],pchara->wns[0]);
    276:     }
    277:     PCApply_HARA(pc,pchara->wns[3],y);
    278:   }
    279:   return(0);
    280: }

    282: static PetscErrorCode MatMult_Hyper(Mat M, Vec x, Vec y)
    283: {

    287:   MatMultKernel_Hyper(M,x,y,PETSC_FALSE);
    288:   return(0);
    289: }

    291: static PetscErrorCode MatMultTranspose_Hyper(Mat M, Vec x, Vec y)
    292: {

    296:   MatMultKernel_Hyper(M,x,y,PETSC_TRUE);
    297:   return(0);
    298: }

    300: /* Hyper power kernel, MatMat version */
    301: static PetscErrorCode MatMatMultKernel_Hyper(Mat M, Mat X, Mat Y, PetscBool t)
    302: {
    303:   PC             pc;
    304:   Mat            A;
    305:   PC_HARA        *pchara;
    306:   PetscInt       i;

    310:   MatShellGetContext(M,(void*)&pc);
    311:   A = pc->useAmat ? pc->mat : pc->pmat;
    312:   pchara = (PC_HARA*)pc->data;
    313:   MatDuplicate(X,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[0]);
    314:   MatDuplicate(Y,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[1]);
    315:   MatDuplicate(X,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[2]);
    316:   MatDuplicate(Y,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[3]);
    317:   MatCopy(X,pchara->wnsmat[0],SAME_NONZERO_PATTERN);
    318:   MatCopy(X,pchara->wnsmat[3],SAME_NONZERO_PATTERN);
    319:   if (t) {
    320:     for (i=0;i<pchara->hyperorder-1;i++) {
    321:       MatTransposeMatMult(A,pchara->wnsmat[0],MAT_REUSE_MATRIX,PETSC_DEFAULT,&pchara->wnsmat[1]);
    322:       PCApplyTransposeMat_HARA(pc,pchara->wnsmat[1],pchara->wnsmat[2]);
    323:       MatAXPY(pchara->wnsmat[0],-1.,pchara->wnsmat[2],SAME_NONZERO_PATTERN);
    324:       MatAXPY(pchara->wnsmat[3],1.,pchara->wnsmat[0],SAME_NONZERO_PATTERN);
    325:     }
    326:     PCApplyTransposeMat_HARA(pc,pchara->wnsmat[3],Y);
    327:   } else {
    328:     for (i=0;i<pchara->hyperorder-1;i++) {
    329:       PCApplyMat_HARA(pc,pchara->wnsmat[0],pchara->wnsmat[1]);
    330:       MatMatMult(A,pchara->wnsmat[1],MAT_REUSE_MATRIX,PETSC_DEFAULT,&pchara->wnsmat[2]);
    331:       MatAXPY(pchara->wnsmat[0],-1.,pchara->wnsmat[2],SAME_NONZERO_PATTERN);
    332:       MatAXPY(pchara->wnsmat[3],1.,pchara->wnsmat[0],SAME_NONZERO_PATTERN);
    333:     }
    334:     PCApplyMat_HARA(pc,pchara->wnsmat[3],Y);
    335:   }
    336:   MatDestroy(&pchara->wnsmat[0]);
    337:   MatDestroy(&pchara->wnsmat[1]);
    338:   MatDestroy(&pchara->wnsmat[2]);
    339:   MatDestroy(&pchara->wnsmat[3]);
    340:   return(0);
    341: }

    343: static PetscErrorCode MatMatMultNumeric_Hyper(Mat M, Mat X, Mat Y,void *ctx)
    344: {

    348:   MatMatMultKernel_Hyper(M,X,Y,PETSC_FALSE);
    349:   return(0);
    350: }

    352: /* Basic Newton-Schultz sampler: (2 * I - M * A) * M */
    353: static PetscErrorCode MatMultKernel_NS(Mat M, Vec x, Vec y, PetscBool t)
    354: {
    355:   PC             pc;
    356:   Mat            A;
    357:   PC_HARA        *pchara;

    361:   MatShellGetContext(M,(void*)&pc);
    362:   A = pc->useAmat ? pc->mat : pc->pmat;
    363:   pchara = (PC_HARA*)pc->data;
    364:   MatCreateVecs(pchara->M,pchara->wns[0] ? NULL : &pchara->wns[0],pchara->wns[1] ? NULL : &pchara->wns[1]);
    365:   if (t) {
    366:     PCApplyTranspose_HARA(pc,x,y);
    367:     MatMultTranspose(A,y,pchara->wns[1]);
    368:     PCApplyTranspose_HARA(pc,pchara->wns[1],pchara->wns[0]);
    369:     VecAXPBY(y,-1.,2.,pchara->wns[0]);
    370:   } else {
    371:     PCApply_HARA(pc,x,y);
    372:     MatMult(A,y,pchara->wns[0]);
    373:     PCApply_HARA(pc,pchara->wns[0],pchara->wns[1]);
    374:     VecAXPBY(y,-1.,2.,pchara->wns[1]);
    375:   }
    376:   return(0);
    377: }

    379: static PetscErrorCode MatMult_NS(Mat M, Vec x, Vec y)
    380: {

    384:   MatMultKernel_NS(M,x,y,PETSC_FALSE);
    385:   return(0);
    386: }

    388: static PetscErrorCode MatMultTranspose_NS(Mat M, Vec x, Vec y)
    389: {

    393:   MatMultKernel_NS(M,x,y,PETSC_TRUE);
    394:   return(0);
    395: }

    397: /* (2 * I - M * A) * M, MatMat version */
    398: static PetscErrorCode MatMatMultKernel_NS(Mat M, Mat X, Mat Y, PetscBool t)
    399: {
    400:   PC             pc;
    401:   Mat            A;
    402:   PC_HARA        *pchara;

    406:   MatShellGetContext(M,(void*)&pc);
    407:   A = pc->useAmat ? pc->mat : pc->pmat;
    408:   pchara = (PC_HARA*)pc->data;
    409:   MatDuplicate(X,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[0]);
    410:   MatDuplicate(Y,MAT_SHARE_NONZERO_PATTERN,&pchara->wnsmat[1]);
    411:   if (t) {
    412:     PCApplyTransposeMat_HARA(pc,X,Y);
    413:     MatTransposeMatMult(A,Y,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pchara->wnsmat[1]);
    414:     PCApplyTransposeMat_HARA(pc,pchara->wnsmat[1],pchara->wnsmat[0]);
    415:     MatScale(Y,2.);
    416:     MatAXPY(Y,-1.,pchara->wnsmat[0],SAME_NONZERO_PATTERN);
    417:   } else {
    418:     PCApplyMat_HARA(pc,X,Y);
    419:     MatMatMult(A,Y,MAT_REUSE_MATRIX,PETSC_DEFAULT,&pchara->wnsmat[0]);
    420:     PCApplyMat_HARA(pc,pchara->wnsmat[0],pchara->wnsmat[1]);
    421:     MatScale(Y,2.);
    422:     MatAXPY(Y,-1.,pchara->wnsmat[1],SAME_NONZERO_PATTERN);
    423:   }
    424:   MatDestroy(&pchara->wnsmat[0]);
    425:   MatDestroy(&pchara->wnsmat[1]);
    426:   return(0);
    427: }

    429: static PetscErrorCode MatMatMultNumeric_NS(Mat M, Mat X, Mat Y, void *)
    430: {

    434:   MatMatMultKernel_NS(M,X,Y,PETSC_FALSE);
    435:   return(0);
    436: }

    438: PETSC_EXTERN PetscErrorCode MatNorm_HARA(Mat,NormType,PetscReal*);

    440: static PetscErrorCode PCHaraSetUpSampler_Private(PC pc)
    441: {
    442:   PC_HARA        *pchara = (PC_HARA*)pc->data;
    443:   Mat            A = pc->useAmat ? pc->mat : pc->pmat;

    447:   if (!pchara->S) {
    448:     PetscInt M,N,m,n;

    450:     MatGetSize(A,&M,&N);
    451:     MatGetLocalSize(A,&m,&n);
    452:     MatCreateShell(PetscObjectComm((PetscObject)A),m,n,M,N,pc,&pchara->S);
    453:     MatSetBlockSizesFromMats(pchara->S,A,A);
    454: #if defined(PETSC_HAVE_CUDA)
    455:     MatShellSetVecType(pchara->S,VECCUDA);
    456: #endif
    457:   }
    458:   if (pchara->hyperorder >= 2) {
    459:     MatShellSetOperation(pchara->S,MATOP_MULT,(void (*)(void))MatMult_Hyper);
    460:     MatShellSetOperation(pchara->S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_Hyper);
    461:     MatShellSetMatProductOperation(pchara->S,MATPRODUCT_AB,NULL,MatMatMultNumeric_Hyper,NULL,MATDENSE,MATDENSE);
    462:     MatShellSetMatProductOperation(pchara->S,MATPRODUCT_AB,NULL,MatMatMultNumeric_Hyper,NULL,MATDENSECUDA,MATDENSECUDA);
    463:   } else {
    464:     MatShellSetOperation(pchara->S,MATOP_MULT,(void (*)(void))MatMult_NS);
    465:     MatShellSetOperation(pchara->S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_NS);
    466:     MatShellSetMatProductOperation(pchara->S,MATPRODUCT_AB,NULL,MatMatMultNumeric_NS,NULL,MATDENSE,MATDENSE);
    467:     MatShellSetMatProductOperation(pchara->S,MATPRODUCT_AB,NULL,MatMatMultNumeric_NS,NULL,MATDENSECUDA,MATDENSECUDA);
    468:   }
    469:   MatPropagateSymmetryOptions(A,pchara->S);
    470:   return(0);
    471: }

    473: /* NS */
    474: static PetscErrorCode PCSetUp_HARA(PC pc)
    475: {
    476:   PC_HARA        *pchara = (PC_HARA*)pc->data;
    477:   Mat            A = pc->useAmat ? pc->mat : pc->pmat;
    479:   NormType       norm = NORM_2;
    480:   PetscReal      initerr,err;

    483:   if (!pchara->T) {
    484:     PetscInt M,N,m,n;

    486:     MatGetSize(pc->pmat,&M,&N);
    487:     MatGetLocalSize(pc->pmat,&m,&n);
    488:     MatCreateShell(PetscObjectComm((PetscObject)pc->pmat),m,n,M,N,pc,&pchara->T);
    489:     MatSetBlockSizesFromMats(pchara->T,pc->pmat,pc->pmat);
    490:     MatShellSetOperation(pchara->T,MATOP_MULT,(void (*)(void))MatMult_MAmI);
    491:     MatShellSetOperation(pchara->T,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_MAmI);
    492:     MatShellSetOperation(pchara->T,MATOP_NORM,(void (*)(void))MatNorm_HARA);
    493: #if defined(PETSC_HAVE_CUDA)
    494:     MatShellSetVecType(pchara->T,VECCUDA);
    495: #endif
    496:     PetscLogObjectParent((PetscObject)pc,(PetscObject)pchara->T);
    497:   }
    498:   if (!pchara->M) {
    499:     Mat       Ain = pc->pmat;
    500:     PetscBool ishara,flg;
    501:     PetscReal onenormA,infnormA;
    502:     void      (*normfunc)(void);

    504:     PetscObjectTypeCompare((PetscObject)Ain,MATHARA,&ishara);
    505:     if (!ishara) {
    506:       Ain  = pc->mat;
    507:       PetscObjectTypeCompare((PetscObject)Ain,MATHARA,&ishara);
    508:     }
    509:     if (!ishara) {
    510:       MatCreateHaraFromMat(A,pchara->sdim,pchara->coords,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&pchara->M);
    511:     } else {
    512:       MatDuplicate(Ain,MAT_SHARE_NONZERO_PATTERN,&pchara->M);
    513:     }

    515:     MatGetOperation(A,MATOP_NORM,&normfunc);
    516:     if (!normfunc || pchara->useapproximatenorms) {
    517:       MatSetOperation(A,MATOP_NORM,(void (*)(void))MatNorm_HARA);
    518:     }
    519:     MatNorm(A,NORM_1,&onenormA);
    520:     MatGetOption(A,MAT_SYMMETRIC,&flg);
    521:     if (!flg) {
    522:       MatNorm(A,NORM_INFINITY,&infnormA);
    523:     } else infnormA = onenormA;
    524:     MatSetOperation(A,MATOP_NORM,normfunc);
    525:     pchara->s0 = 1./(infnormA*onenormA);
    526:   }
    527:   MatNorm(pchara->T,norm,&initerr);
    528:   if (initerr > pchara->atol) {
    529:     PetscInt i;

    531:     PCHaraSetUpSampler_Private(pc);
    532:     err  = initerr;
    533:     if (pchara->monitor) { PetscPrintf(PetscObjectComm((PetscObject)pc),"%D: %g %g\n",0,(double)err,(double)(err/initerr)); }
    534:     for (i = 0; i < pchara->maxits; i++) {
    535:       Mat         M;
    536:       const char* prefix;

    538:       MatDuplicate(pchara->M,MAT_SHARE_NONZERO_PATTERN,&M);
    539:       MatGetOptionsPrefix(M,&prefix);
    540:       if (!prefix) {
    541:         PCGetOptionsPrefix(pc,&prefix);
    542:         MatSetOptionsPrefix(M,prefix);
    543:         MatAppendOptionsPrefix(M,"pc_hara_inv_");
    544:       }
    545: #if 0
    546:   {
    547:      Mat Sd1,Sd2,Id;
    548:      PetscReal err;
    549:      MatComputeOperator(pchara->S,MATDENSE,&Sd1);
    550:      MatDuplicate(Sd1,MAT_COPY_VALUES,&Id);
    551:      MatZeroEntries(Id);
    552:      MatShift(Id,1.);
    553:      MatMatMult(pchara->S,Id,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Sd2);
    554:      MatAXPY(Sd2,-1.,Sd1,SAME_NONZERO_PATTERN);
    555:      MatNorm(Sd2,NORM_FROBENIUS,&err);
    556:      PetscPrintf(PetscObjectComm((PetscObject)Sd2),"ERR %g\n",err);
    557:      MatViewFromOptions(Sd2,NULL,"-Sd_view");
    558:      MatDestroy(&Sd1);
    559:      MatDestroy(&Sd2);
    560:      MatDestroy(&Id);
    561:   }
    562: #endif
    563:       MatHaraSetSamplingMat(M,pchara->S,1,PETSC_DECIDE);
    564:       if (pc->setfromoptionscalled) {
    565:         MatSetFromOptions(M);
    566:       }
    567:       MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
    568:       MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
    569: #if 0
    570:       {
    571:          Mat Md;
    572:          MatComputeOperator(M,MATDENSE,&Md);
    573:          MatViewFromOptions(Md,NULL,"-Md_view");
    574:          MatDestroy(&Md);
    575:          MatComputeOperator(pchara->S,MATDENSE,&Md);
    576:          MatViewFromOptions(Md,NULL,"-Md_view");
    577:          MatDestroy(&Md);
    578:       }
    579: #endif
    580:       MatDestroy(&pchara->M);
    581:       pchara->M = M;
    582:       MatNorm(pchara->T,norm,&err);
    583:       if (pchara->monitor) { PetscPrintf(PetscObjectComm((PetscObject)pc),"%D: %g %g\n",i+1,(double)err,(double)(err/initerr)); }
    584:       if (err < pchara->atol || err < pchara->rtol*initerr) break;
    585:     }
    586:   }
    587:   return(0);
    588: }

    590: PETSC_EXTERN PetscErrorCode PCCreate_HARA(PC pc)
    591: {
    593:   PC_HARA        *pchara;

    596:   PetscNewLog(pc,&pchara);
    597:   pc->data = (void*)pchara;

    599:   pchara->atol       = 1.e-2;
    600:   pchara->rtol       = 1.e-6;
    601:   pchara->maxits     = 50;
    602:   pchara->hyperorder = 1; /* default to basic NewtonSchultz */

    604:   pc->ops->destroy        = PCDestroy_HARA;
    605:   pc->ops->setup          = PCSetUp_HARA;
    606:   pc->ops->apply          = PCApply_HARA;
    607:   pc->ops->matapply       = PCApplyMat_HARA;
    608:   pc->ops->applytranspose = PCApplyTranspose_HARA;
    609:   pc->ops->reset          = PCReset_HARA;
    610:   pc->ops->setfromoptions = PCSetFromOptions_HARA;

    612:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HARA);
    613:   return(0);
    614: }