Actual source code: pchara.cu

petsc-master 2020-10-29
Report Typos and Errors
  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: }