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