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: }