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