Actual source code: fcg.c
1: /*
2: This file implements the FCG (Flexible Conjugate Gradient) method
4: */
6: #include <../src/ksp/ksp/impls/fcg/fcgimpl.h>
7: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP,PetscReal*,PetscReal*);
8: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt*);
10: #define KSPFCG_DEFAULT_MMAX 30 /* maximum number of search directions to keep */
11: #define KSPFCG_DEFAULT_NPREALLOC 10 /* number of search directions to preallocate */
12: #define KSPFCG_DEFAULT_VECB 5 /* number of search directions to allocate each time new direction vectors are needed */
13: #define KSPFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
15: static PetscErrorCode KSPAllocateVectors_FCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
16: {
17: PetscInt i;
18: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
19: PetscInt nnewvecs, nvecsprev;
21: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
22: if (fcg->nvecs < PetscMin(fcg->mmax+1,nvecsneeded)) {
23: nvecsprev = fcg->nvecs;
24: nnewvecs = PetscMin(PetscMax(nvecsneeded-fcg->nvecs,chunksize),fcg->mmax+1-fcg->nvecs);
25: KSPCreateVecs(ksp,nnewvecs,&fcg->pCvecs[fcg->nchunks],0,NULL);
26: PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pCvecs[fcg->nchunks]);
27: KSPCreateVecs(ksp,nnewvecs,&fcg->pPvecs[fcg->nchunks],0,NULL);
28: PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pPvecs[fcg->nchunks]);
29: fcg->nvecs += nnewvecs;
30: for (i=0;i<nnewvecs;++i) {
31: fcg->Cvecs[nvecsprev + i] = fcg->pCvecs[fcg->nchunks][i];
32: fcg->Pvecs[nvecsprev + i] = fcg->pPvecs[fcg->nchunks][i];
33: }
34: fcg->chunksizes[fcg->nchunks] = nnewvecs;
35: ++fcg->nchunks;
36: }
37: return 0;
38: }
40: static PetscErrorCode KSPSetUp_FCG(KSP ksp)
41: {
42: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
43: PetscInt maxit = ksp->max_it;
44: const PetscInt nworkstd = 2;
47: /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
48: KSPSetWorkVecs(ksp,nworkstd);
50: /* Allocated space for pointers to additional work vectors
51: note that mmax is the number of previous directions, so we add 1 for the current direction,
52: and an extra 1 for the prealloc (which might be empty) */
53: PetscMalloc5(fcg->mmax+1,&fcg->Pvecs,fcg->mmax+1,&fcg->Cvecs,fcg->mmax+1,&fcg->pPvecs,fcg->mmax+1,&fcg->pCvecs,fcg->mmax+2,&fcg->chunksizes);
54: PetscLogObjectMemory((PetscObject)ksp,2*(fcg->mmax+1)*sizeof(Vec*) + 2*(fcg->mmax + 1)*sizeof(Vec**) + (fcg->mmax + 2)*sizeof(PetscInt));
56: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
57: if (fcg->nprealloc > fcg->mmax+1) {
58: PetscInfo(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",fcg->nprealloc, fcg->mmax+1);
59: }
61: /* Preallocate additional work vectors */
62: KSPAllocateVectors_FCG(ksp,fcg->nprealloc,fcg->nprealloc);
63: /*
64: If user requested computations of eigenvalues then allocate work
65: work space needed
66: */
67: if (ksp->calc_sings) {
68: /* get space to store tridiagonal matrix for Lanczos */
69: PetscMalloc4(maxit,&fcg->e,maxit,&fcg->d,maxit,&fcg->ee,maxit,&fcg->dd);
70: PetscLogObjectMemory((PetscObject)ksp,2*(maxit+1)*(sizeof(PetscScalar)+sizeof(PetscReal)));
72: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
73: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
74: }
75: return 0;
76: }
78: static PetscErrorCode KSPSolve_FCG(KSP ksp)
79: {
80: PetscInt i,k,idx,mi;
81: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
82: PetscScalar alpha=0.0,beta = 0.0,dpi,s;
83: PetscReal dp=0.0;
84: Vec B,R,Z,X,Pcurr,Ccurr;
85: Mat Amat,Pmat;
86: PetscInt eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
87: PetscInt stored_max_it = ksp->max_it;
88: PetscScalar alphaold = 0,betaold = 1.0,*e = NULL,*d = NULL;/* Variables for eigen estimation - FINISH */
91: #define VecXDot(x,y,a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))
92: #define VecXMDot(a,b,c,d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a,b,c,d) : VecMTDot(a,b,c,d))
94: X = ksp->vec_sol;
95: B = ksp->vec_rhs;
96: R = ksp->work[0];
97: Z = ksp->work[1];
99: PCGetOperators(ksp->pc,&Amat,&Pmat);
100: if (eigs) {e = fcg->e; d = fcg->d; e[0] = 0.0; }
101: /* Compute initial residual needed for convergence check*/
102: ksp->its = 0;
103: if (!ksp->guess_zero) {
104: KSP_MatMult(ksp,Amat,X,R);
105: VecAYPX(R,-1.0,B); /* r <- b - Ax */
106: } else {
107: VecCopy(B,R); /* r <- b (x is 0) */
108: }
109: switch (ksp->normtype) {
110: case KSP_NORM_PRECONDITIONED:
111: KSP_PCApply(ksp,R,Z); /* z <- Br */
112: VecNorm(Z,NORM_2,&dp); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
113: KSPCheckNorm(ksp,dp);
114: break;
115: case KSP_NORM_UNPRECONDITIONED:
116: VecNorm(R,NORM_2,&dp); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
117: KSPCheckNorm(ksp,dp);
118: break;
119: case KSP_NORM_NATURAL:
120: KSP_PCApply(ksp,R,Z); /* z <- Br */
121: VecXDot(R,Z,&s);
122: KSPCheckDot(ksp,s);
123: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
124: break;
125: case KSP_NORM_NONE:
126: dp = 0.0;
127: break;
128: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
129: }
131: /* Initial Convergence Check */
132: KSPLogResidualHistory(ksp,dp);
133: KSPMonitor(ksp,0,dp);
134: ksp->rnorm = dp;
135: if (ksp->normtype == KSP_NORM_NONE) {
136: KSPConvergedSkip(ksp,0,dp,&ksp->reason,ksp->cnvP);
137: } else {
138: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
139: }
140: if (ksp->reason) return 0;
142: /* Apply PC if not already done for convergence check */
143: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) {
144: KSP_PCApply(ksp,R,Z); /* z <- Br */
145: }
147: i = 0;
148: do {
149: ksp->its = i+1;
151: /* If needbe, allocate a new chunk of vectors in P and C */
152: KSPAllocateVectors_FCG(ksp,i+1,fcg->vecb);
154: /* Note that we wrap around and start clobbering old vectors */
155: idx = i % (fcg->mmax+1);
156: Pcurr = fcg->Pvecs[idx];
157: Ccurr = fcg->Cvecs[idx];
159: /* number of old directions to orthogonalize against */
160: switch(fcg->truncstrat) {
161: case KSP_FCD_TRUNC_TYPE_STANDARD:
162: mi = fcg->mmax;
163: break;
164: case KSP_FCD_TRUNC_TYPE_NOTAY:
165: mi = ((i-1) % fcg->mmax)+1;
166: break;
167: default:
168: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy");
169: }
171: /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
172: VecCopy(Z,Pcurr);
174: {
175: PetscInt l,ndots;
177: l = PetscMax(0,i-mi);
178: ndots = i-l;
179: if (ndots) {
180: PetscInt j;
181: Vec *Pold, *Cold;
182: PetscScalar *dots;
184: PetscMalloc3(ndots,&dots,ndots,&Cold,ndots,&Pold);
185: for (k=l,j=0;j<ndots;++k,++j) {
186: idx = k % (fcg->mmax+1);
187: Cold[j] = fcg->Cvecs[idx];
188: Pold[j] = fcg->Pvecs[idx];
189: }
190: VecXMDot(Z,ndots,Cold,dots);
191: for (k=0;k<ndots;++k) {
192: dots[k] = -dots[k];
193: }
194: VecMAXPY(Pcurr,ndots,dots,Pold);
195: PetscFree3(dots,Cold,Pold);
196: }
197: }
199: /* Update X and R */
200: betaold = beta;
201: VecXDot(Pcurr,R,&beta); /* beta <- pi'*r */
202: KSPCheckDot(ksp,beta);
203: KSP_MatMult(ksp,Amat,Pcurr,Ccurr); /* w <- A*pi (stored in ci) */
204: VecXDot(Pcurr,Ccurr,&dpi); /* dpi <- pi'*w */
205: alphaold = alpha;
206: alpha = beta / dpi; /* alpha <- beta/dpi */
207: VecAXPY(X,alpha,Pcurr); /* x <- x + alpha * pi */
208: VecAXPY(R,-alpha,Ccurr); /* r <- r - alpha * wi */
210: /* Compute norm for convergence check */
211: switch (ksp->normtype) {
212: case KSP_NORM_PRECONDITIONED:
213: KSP_PCApply(ksp,R,Z); /* z <- Br */
214: VecNorm(Z,NORM_2,&dp); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
215: KSPCheckNorm(ksp,dp);
216: break;
217: case KSP_NORM_UNPRECONDITIONED:
218: VecNorm(R,NORM_2,&dp); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
219: KSPCheckNorm(ksp,dp);
220: break;
221: case KSP_NORM_NATURAL:
222: KSP_PCApply(ksp,R,Z); /* z <- Br */
223: VecXDot(R,Z,&s);
224: KSPCheckDot(ksp,s);
225: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
226: break;
227: case KSP_NORM_NONE:
228: dp = 0.0;
229: break;
230: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
231: }
233: /* Check for convergence */
234: ksp->rnorm = dp;
235: KSPLogResidualHistory(ksp,dp);
236: KSPMonitor(ksp,i+1,dp);
237: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
238: if (ksp->reason) break;
240: /* Apply PC if not already done for convergence check */
241: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) {
242: KSP_PCApply(ksp,R,Z); /* z <- Br */
243: }
245: /* Compute current C (which is W/dpi) */
246: VecScale(Ccurr,1.0/dpi); /* w <- ci/dpi */
248: if (eigs) {
249: if (i > 0) {
251: e[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))/alphaold;
252: d[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))*e[i] + 1.0/alpha;
253: } else {
254: d[i] = PetscSqrtReal(PetscAbsScalar(beta))*e[i] + 1.0/alpha;
255: }
256: fcg->ned = ksp->its-1;
257: }
258: ++i;
259: } while (i<ksp->max_it);
260: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
261: return 0;
262: }
264: static PetscErrorCode KSPDestroy_FCG(KSP ksp)
265: {
266: PetscInt i;
267: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
270: /* Destroy "standard" work vecs */
271: VecDestroyVecs(ksp->nwork,&ksp->work);
273: /* Destroy P and C vectors and the arrays that manage pointers to them */
274: if (fcg->nvecs) {
275: for (i=0;i<fcg->nchunks;++i) {
276: VecDestroyVecs(fcg->chunksizes[i],&fcg->pPvecs[i]);
277: VecDestroyVecs(fcg->chunksizes[i],&fcg->pCvecs[i]);
278: }
279: }
280: PetscFree5(fcg->Pvecs,fcg->Cvecs,fcg->pPvecs,fcg->pCvecs,fcg->chunksizes);
281: /* free space used for singular value calculations */
282: if (ksp->calc_sings) {
283: PetscFree4(fcg->e,fcg->d,fcg->ee,fcg->dd);
284: }
285: KSPDestroyDefault(ksp);
286: return 0;
287: }
289: static PetscErrorCode KSPView_FCG(KSP ksp,PetscViewer viewer)
290: {
291: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
292: PetscBool iascii,isstring;
293: const char *truncstr;
295: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
296: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
298: if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy";
299: else if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy";
300: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCG truncation strategy");
302: if (iascii) {
303: PetscViewerASCIIPrintf(viewer," m_max=%D\n",fcg->mmax);
304: PetscViewerASCIIPrintf(viewer," preallocated %D directions\n",PetscMin(fcg->nprealloc,fcg->mmax+1));
305: PetscViewerASCIIPrintf(viewer," %s\n",truncstr);
306: } else if (isstring) {
307: PetscViewerStringSPrintf(viewer,"m_max %D nprealloc %D %s",fcg->mmax,fcg->nprealloc,truncstr);
308: }
309: return 0;
310: }
312: /*@
313: KSPFCGSetMmax - set the maximum number of previous directions FCG will store for orthogonalization
315: Note: mmax + 1 directions are stored (mmax previous ones along with a current one)
316: and whether all are used in each iteration also depends on the truncation strategy
317: (see KSPFCGSetTruncationType())
319: Logically Collective on ksp
321: Input Parameters:
322: + ksp - the Krylov space context
323: - mmax - the maximum number of previous directions to orthogonalize againt
325: Level: intermediate
327: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc()
328: @*/
329: PetscErrorCode KSPFCGSetMmax(KSP ksp,PetscInt mmax)
330: {
331: KSP_FCG *fcg = (KSP_FCG*)ksp->data;
335: fcg->mmax = mmax;
336: return 0;
337: }
339: /*@
340: KSPFCGGetMmax - get the maximum number of previous directions FCG will store
342: Note: FCG stores mmax+1 directions at most (mmax previous ones, and one current one)
344: Not Collective
346: Input Parameter:
347: . ksp - the Krylov space context
349: Output Parameter:
350: . mmax - the maximum number of previous directions allowed for orthogonalization
352: Level: intermediate
354: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc(), KSPFCGSetMmax()
355: @*/
357: PetscErrorCode KSPFCGGetMmax(KSP ksp,PetscInt *mmax)
358: {
359: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
362: *mmax = fcg->mmax;
363: return 0;
364: }
366: /*@
367: KSPFCGSetNprealloc - set the number of directions to preallocate with FCG
369: Logically Collective on ksp
371: Input Parameters:
372: + ksp - the Krylov space context
373: - nprealloc - the number of vectors to preallocate
375: Level: advanced
377: Options Database:
378: . -ksp_fcg_nprealloc <N> - number of directions to preallocate
380: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc()
381: @*/
382: PetscErrorCode KSPFCGSetNprealloc(KSP ksp,PetscInt nprealloc)
383: {
384: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
389: fcg->nprealloc = nprealloc;
390: return 0;
391: }
393: /*@
394: KSPFCGGetNprealloc - get the number of directions preallocate by FCG
396: Not Collective
398: Input Parameter:
399: . ksp - the Krylov space context
401: Output Parameter:
402: . nprealloc - the number of directions preallocated
404: Level: advanced
406: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGSetNprealloc()
407: @*/
408: PetscErrorCode KSPFCGGetNprealloc(KSP ksp,PetscInt *nprealloc)
409: {
410: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
413: *nprealloc = fcg->nprealloc;
414: return 0;
415: }
417: /*@
418: KSPFCGSetTruncationType - specify how many of its stored previous directions FCG uses during orthoganalization
420: Logically Collective on ksp
422: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
423: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
425: Input Parameters:
426: + ksp - the Krylov space context
427: - truncstrat - the choice of strategy
429: Level: intermediate
431: Options Database:
432: . -ksp_fcg_truncation_type <standard, notay> - specify how many of its stored previous directions FCG uses during orthoganalization
434: .seealso: KSPFCDTruncationType, KSPFCGGetTruncationType
435: @*/
436: PetscErrorCode KSPFCGSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat)
437: {
438: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
442: fcg->truncstrat=truncstrat;
443: return 0;
444: }
446: /*@
447: KSPFCGGetTruncationType - get the truncation strategy employed by FCG
449: Not Collective
451: Input Parameter:
452: . ksp - the Krylov space context
454: Output Parameter:
455: . truncstrat - the strategy type
457: Level: intermediate
459: .seealso: KSPFCG, KSPFCGSetTruncationType, KSPFCDTruncationType
460: @*/
461: PetscErrorCode KSPFCGGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat)
462: {
463: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
466: *truncstrat=fcg->truncstrat;
467: return 0;
468: }
470: static PetscErrorCode KSPSetFromOptions_FCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
471: {
472: KSP_FCG *fcg=(KSP_FCG*)ksp->data;
473: PetscInt mmax,nprealloc;
474: PetscBool flg;
476: PetscOptionsHead(PetscOptionsObject,"KSP FCG Options");
477: PetscOptionsInt("-ksp_fcg_mmax","Maximum number of search directions to store","KSPFCGSetMmax",fcg->mmax,&mmax,&flg);
478: if (flg) {
479: KSPFCGSetMmax(ksp,mmax);
480: }
481: PetscOptionsInt("-ksp_fcg_nprealloc","Number of directions to preallocate","KSPFCGSetNprealloc",fcg->nprealloc,&nprealloc,&flg);
482: if (flg) {
483: KSPFCGSetNprealloc(ksp,nprealloc);
484: }
485: PetscOptionsEnum("-ksp_fcg_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)fcg->truncstrat,(PetscEnum*)&fcg->truncstrat,NULL);
486: PetscOptionsTail();
487: return 0;
488: }
490: /*MC
491: KSPFCG - Implements the Flexible Conjugate Gradient method (FCG)
493: Options Database Keys:
494: + -ksp_fcg_mmax <N> - maximum number of search directions
495: . -ksp_fcg_nprealloc <N> - number of directions to preallocate
496: - -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions
498: Contributed by Patrick Sanan
500: Notes:
501: Supports left preconditioning only.
503: Level: beginner
505: References:
506: + * - Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, 2000
507: - * - Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable step Preconditioning",
508: SIAM J. Matrix Anal. Appl. 12:4, 1991
510: .seealso : KSPGCR, KSPFGMRES, KSPCG, KSPFCGSetMmax(), KSPFCGGetMmax(), KSPFCGSetNprealloc(), KSPFCGGetNprealloc(), KSPFCGSetTruncationType(), KSPFCGGetTruncationType()
512: M*/
513: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
514: {
515: KSP_FCG *fcg;
517: PetscNewLog(ksp,&fcg);
518: #if !defined(PETSC_USE_COMPLEX)
519: fcg->type = KSP_CG_SYMMETRIC;
520: #else
521: fcg->type = KSP_CG_HERMITIAN;
522: #endif
523: fcg->mmax = KSPFCG_DEFAULT_MMAX;
524: fcg->nprealloc = KSPFCG_DEFAULT_NPREALLOC;
525: fcg->nvecs = 0;
526: fcg->vecb = KSPFCG_DEFAULT_VECB;
527: fcg->nchunks = 0;
528: fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;
530: ksp->data = (void*)fcg;
532: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
533: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
534: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
535: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
537: ksp->ops->setup = KSPSetUp_FCG;
538: ksp->ops->solve = KSPSolve_FCG;
539: ksp->ops->destroy = KSPDestroy_FCG;
540: ksp->ops->view = KSPView_FCG;
541: ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
542: ksp->ops->buildsolution = KSPBuildSolutionDefault;
543: ksp->ops->buildresidual = KSPBuildResidualDefault;
544: return 0;
545: }