1 #include <../src/ksp/ksp/impls/cg/stcg/stcgimpl.h> /*I "petscksp.h" I*/
2
3 #define STCG_PRECONDITIONED_DIRECTION 0
4 #define STCG_UNPRECONDITIONED_DIRECTION 1
5 #define STCG_DIRECTION_TYPES 2
6
7 static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
8
KSPCGSolve_STCG(KSP ksp)9 static PetscErrorCode KSPCGSolve_STCG(KSP ksp)
10 {
11 #if defined(PETSC_USE_COMPLEX)
12 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "STCG is not available for complex systems");
13 #else
14 KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
15 Mat Qmat, Mmat;
16 Vec r, z, p, d;
17 PC pc;
18 PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
19 PetscReal alpha, beta, kappa, rz, rzm1;
20 PetscReal rr, r2, step;
21 PetscInt max_cg_its;
22 PetscBool diagonalscale;
23
24 /***************************************************************************/
25 /* Check the arguments and parameters. */
26 /***************************************************************************/
27
28 PetscFunctionBegin;
29 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
30 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
31 PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
32
33 /***************************************************************************/
34 /* Get the workspace vectors and initialize variables */
35 /***************************************************************************/
36
37 r2 = cg->radius * cg->radius;
38 r = ksp->work[0];
39 z = ksp->work[1];
40 p = ksp->work[2];
41 d = ksp->vec_sol;
42 pc = ksp->pc;
43
44 PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
45
46 PetscCall(VecGetSize(d, &max_cg_its));
47 max_cg_its = PetscMin(max_cg_its, ksp->max_it);
48 ksp->its = 0;
49
50 /***************************************************************************/
51 /* Initialize objective function and direction. */
52 /***************************************************************************/
53
54 cg->o_fcn = 0.0;
55
56 PetscCall(VecSet(d, 0.0)); /* d = 0 */
57 cg->norm_d = 0.0;
58
59 /***************************************************************************/
60 /* Begin the conjugate gradient method. Check the right-hand side for */
61 /* numerical problems. The check for not-a-number and infinite values */
62 /* need be performed only once. */
63 /***************************************************************************/
64
65 PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
66 PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
67 KSPCheckDot(ksp, rr);
68
69 /***************************************************************************/
70 /* Check the preconditioner for numerical problems and for positive */
71 /* definiteness. The check for not-a-number and infinite values need be */
72 /* performed only once. */
73 /***************************************************************************/
74
75 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
76 PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
77 if (PetscIsInfOrNanScalar(rz)) {
78 /*************************************************************************/
79 /* The preconditioner contains not-a-number or an infinite value. */
80 /* Return the gradient direction intersected with the trust region. */
81 /*************************************************************************/
82
83 ksp->reason = KSP_DIVERGED_NANORINF;
84 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad preconditioner: rz=%g\n", (double)rz));
85
86 if (cg->radius != 0) {
87 if (r2 >= rr) {
88 alpha = 1.0;
89 cg->norm_d = PetscSqrtReal(rr);
90 } else {
91 alpha = PetscSqrtReal(r2 / rr);
92 cg->norm_d = cg->radius;
93 }
94
95 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
96
97 /***********************************************************************/
98 /* Compute objective function. */
99 /***********************************************************************/
100
101 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
102 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
103 PetscCall(VecDot(d, z, &cg->o_fcn));
104 cg->o_fcn = -cg->o_fcn;
105 ++ksp->its;
106 }
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
110 if (rz < 0.0) {
111 /*************************************************************************/
112 /* The preconditioner is indefinite. Because this is the first */
113 /* and we do not have a direction yet, we use the gradient step. Note */
114 /* that we cannot use the preconditioned norm when computing the step */
115 /* because the matrix is indefinite. */
116 /*************************************************************************/
117
118 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
119 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: indefinite preconditioner: rz=%g\n", (double)rz));
120
121 if (cg->radius != 0.0) {
122 if (r2 >= rr) {
123 alpha = 1.0;
124 cg->norm_d = PetscSqrtReal(rr);
125 } else {
126 alpha = PetscSqrtReal(r2 / rr);
127 cg->norm_d = cg->radius;
128 }
129
130 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
131
132 /***********************************************************************/
133 /* Compute objective function. */
134 /***********************************************************************/
135
136 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
137 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
138 PetscCall(VecDot(d, z, &cg->o_fcn));
139 cg->o_fcn = -cg->o_fcn;
140 ++ksp->its;
141 }
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
145 /***************************************************************************/
146 /* As far as we know, the preconditioner is positive semidefinite. */
147 /* Compute and log the residual. Check convergence because this */
148 /* initializes things, but do not terminate until at least one conjugate */
149 /* gradient iteration has been performed. */
150 /***************************************************************************/
151
152 switch (ksp->normtype) {
153 case KSP_NORM_PRECONDITIONED:
154 PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
155 break;
156
157 case KSP_NORM_UNPRECONDITIONED:
158 norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
159 break;
160
161 case KSP_NORM_NATURAL:
162 norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
163 break;
164
165 default:
166 norm_r = 0.0;
167 break;
168 }
169
170 PetscCall(KSPLogResidualHistory(ksp, norm_r));
171 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
172 ksp->rnorm = norm_r;
173
174 PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
175
176 /***************************************************************************/
177 /* Compute the first direction and update the iteration. */
178 /***************************************************************************/
179
180 PetscCall(VecCopy(z, p)); /* p = z */
181 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
182 ++ksp->its;
183
184 /***************************************************************************/
185 /* Check the matrix for numerical problems. */
186 /***************************************************************************/
187
188 PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
189 if (PetscIsInfOrNanScalar(kappa)) {
190 /*************************************************************************/
191 /* The matrix produced not-a-number or an infinite value. In this case, */
192 /* we must stop and use the gradient direction. This condition need */
193 /* only be checked once. */
194 /*************************************************************************/
195
196 ksp->reason = KSP_DIVERGED_NANORINF;
197 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: bad matrix: kappa=%g\n", (double)kappa));
198
199 if (cg->radius) {
200 if (r2 >= rr) {
201 alpha = 1.0;
202 cg->norm_d = PetscSqrtReal(rr);
203 } else {
204 alpha = PetscSqrtReal(r2 / rr);
205 cg->norm_d = cg->radius;
206 }
207
208 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
209
210 /***********************************************************************/
211 /* Compute objective function. */
212 /***********************************************************************/
213
214 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
215 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
216 PetscCall(VecDot(d, z, &cg->o_fcn));
217 cg->o_fcn = -cg->o_fcn;
218 ++ksp->its;
219 }
220 PetscFunctionReturn(PETSC_SUCCESS);
221 }
222
223 /***************************************************************************/
224 /* Initialize variables for calculating the norm of the direction. */
225 /***************************************************************************/
226
227 dMp = 0.0;
228 norm_d = 0.0;
229 switch (cg->dtype) {
230 case STCG_PRECONDITIONED_DIRECTION:
231 norm_p = rz;
232 break;
233
234 default:
235 PetscCall(VecDot(p, p, &norm_p));
236 break;
237 }
238
239 /***************************************************************************/
240 /* Check for negative curvature. */
241 /***************************************************************************/
242
243 if (kappa <= 0.0) {
244 /*************************************************************************/
245 /* In this case, the matrix is indefinite and we have encountered a */
246 /* direction of negative curvature. Because negative curvature occurs */
247 /* during the first step, we must follow a direction. */
248 /*************************************************************************/
249
250 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
251 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
252
253 if (cg->radius != 0.0 && norm_p > 0.0) {
254 /***********************************************************************/
255 /* Follow direction of negative curvature to the boundary of the */
256 /* trust region. */
257 /***********************************************************************/
258
259 step = PetscSqrtReal(r2 / norm_p);
260 cg->norm_d = cg->radius;
261
262 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
263
264 /***********************************************************************/
265 /* Update objective function. */
266 /***********************************************************************/
267
268 cg->o_fcn += step * (0.5 * step * kappa - rz);
269 } else if (cg->radius != 0.0) {
270 /***********************************************************************/
271 /* The norm of the preconditioned direction is zero; use the gradient */
272 /* step. */
273 /***********************************************************************/
274
275 if (r2 >= rr) {
276 alpha = 1.0;
277 cg->norm_d = PetscSqrtReal(rr);
278 } else {
279 alpha = PetscSqrtReal(r2 / rr);
280 cg->norm_d = cg->radius;
281 }
282
283 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
284
285 /***********************************************************************/
286 /* Compute objective function. */
287 /***********************************************************************/
288
289 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
290 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
291 PetscCall(VecDot(d, z, &cg->o_fcn));
292
293 cg->o_fcn = -cg->o_fcn;
294 ++ksp->its;
295 }
296 PetscFunctionReturn(PETSC_SUCCESS);
297 }
298
299 /***************************************************************************/
300 /* Run the conjugate gradient method until either the problem is solved, */
301 /* we encounter the boundary of the trust region, or the conjugate */
302 /* gradient method breaks down. */
303 /***************************************************************************/
304
305 while (1) {
306 /*************************************************************************/
307 /* Know that kappa is nonzero, because we have not broken down, so we */
308 /* can compute the steplength. */
309 /*************************************************************************/
310
311 alpha = rz / kappa;
312
313 /*************************************************************************/
314 /* Compute the steplength and check for intersection with the trust */
315 /* region. */
316 /*************************************************************************/
317
318 norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
319 if (cg->radius != 0.0 && norm_dp1 >= r2) {
320 /***********************************************************************/
321 /* In this case, the matrix is positive definite as far as we know. */
322 /* However, the full step goes beyond the trust region. */
323 /***********************************************************************/
324
325 ksp->reason = KSP_CONVERGED_STEP_LENGTH;
326 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: constrained step: radius=%g\n", (double)cg->radius));
327
328 if (norm_p > 0.0) {
329 /*********************************************************************/
330 /* Follow the direction to the boundary of the trust region. */
331 /* Final residual norm is never computed. */
332 /*********************************************************************/
333
334 step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
335 cg->norm_d = cg->radius;
336
337 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
338
339 /*********************************************************************/
340 /* Update objective function. */
341 /*********************************************************************/
342
343 cg->o_fcn += step * (0.5 * step * kappa - rz);
344 } else {
345 /*********************************************************************/
346 /* The norm of the direction is zero; there is nothing to follow. */
347 /*********************************************************************/
348 }
349 break;
350 }
351
352 /*************************************************************************/
353 /* Now we can update the direction and residual. */
354 /*************************************************************************/
355
356 PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
357 PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
358 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
359
360 switch (cg->dtype) {
361 case STCG_PRECONDITIONED_DIRECTION:
362 norm_d = norm_dp1;
363 break;
364
365 default:
366 PetscCall(VecDot(d, d, &norm_d));
367 break;
368 }
369 cg->norm_d = PetscSqrtReal(norm_d);
370
371 /*************************************************************************/
372 /* Update objective function. */
373 /*************************************************************************/
374
375 cg->o_fcn -= 0.5 * alpha * rz;
376
377 /*************************************************************************/
378 /* Check that the preconditioner appears positive semidefinite. */
379 /*************************************************************************/
380
381 rzm1 = rz;
382 PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
383 if (rz < 0.0) {
384 /***********************************************************************/
385 /* The preconditioner is indefinite. */
386 /***********************************************************************/
387
388 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
389 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: cg indefinite preconditioner: rz=%g\n", (double)rz));
390 break;
391 }
392
393 /*************************************************************************/
394 /* As far as we know, the preconditioner is positive semidefinite. */
395 /* Compute the residual and check for convergence. */
396 /*************************************************************************/
397
398 switch (ksp->normtype) {
399 case KSP_NORM_PRECONDITIONED:
400 PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
401 break;
402
403 case KSP_NORM_UNPRECONDITIONED:
404 PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
405 break;
406
407 case KSP_NORM_NATURAL:
408 norm_r = PetscSqrtReal(rz); /* norm_r = |r|_M */
409 break;
410
411 default:
412 norm_r = 0.0;
413 break;
414 }
415
416 PetscCall(KSPLogResidualHistory(ksp, norm_r));
417 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
418 ksp->rnorm = norm_r;
419
420 PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
421 if (ksp->reason) {
422 /***********************************************************************/
423 /* The method has converged. */
424 /***********************************************************************/
425
426 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
427 break;
428 }
429
430 /*************************************************************************/
431 /* We have not converged yet. Check for breakdown. */
432 /*************************************************************************/
433
434 beta = rz / rzm1;
435 if (PetscAbsScalar(beta) <= 0.0) {
436 /***********************************************************************/
437 /* Conjugate gradients has broken down. */
438 /***********************************************************************/
439
440 ksp->reason = KSP_DIVERGED_BREAKDOWN;
441 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: breakdown: beta=%g\n", (double)beta));
442 break;
443 }
444
445 /*************************************************************************/
446 /* Check iteration limit. */
447 /*************************************************************************/
448
449 if (ksp->its >= max_cg_its) {
450 ksp->reason = KSP_DIVERGED_ITS;
451 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
452 break;
453 }
454
455 /*************************************************************************/
456 /* Update p and the norms. */
457 /*************************************************************************/
458
459 PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
460
461 switch (cg->dtype) {
462 case STCG_PRECONDITIONED_DIRECTION:
463 dMp = beta * (dMp + alpha * norm_p);
464 norm_p = beta * (rzm1 + beta * norm_p);
465 break;
466
467 default:
468 PetscCall(VecDot(d, p, &dMp));
469 PetscCall(VecDot(p, p, &norm_p));
470 break;
471 }
472
473 /*************************************************************************/
474 /* Compute the new direction and update the iteration. */
475 /*************************************************************************/
476
477 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
478 PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
479 ++ksp->its;
480
481 /*************************************************************************/
482 /* Check for negative curvature. */
483 /*************************************************************************/
484
485 if (kappa <= 0.0) {
486 /***********************************************************************/
487 /* In this case, the matrix is indefinite and we have encountered */
488 /* a direction of negative curvature. Follow the direction to the */
489 /* boundary of the trust region. */
490 /***********************************************************************/
491
492 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
493 PetscCall(PetscInfo(ksp, "KSPCGSolve_STCG: negative curvature: kappa=%g\n", (double)kappa));
494
495 if (cg->radius != 0.0 && norm_p > 0.0) {
496 /*********************************************************************/
497 /* Follow direction of negative curvature to boundary. */
498 /* Final residual norm is never computed. */
499 /*********************************************************************/
500
501 step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
502 cg->norm_d = cg->radius;
503
504 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
505
506 /*********************************************************************/
507 /* Update objective function. */
508 /*********************************************************************/
509
510 cg->o_fcn += step * (0.5 * step * kappa - rz);
511 } else if (cg->radius != 0.0) {
512 /*********************************************************************/
513 /* The norm of the direction is zero; there is nothing to follow. */
514 /*********************************************************************/
515 }
516 break;
517 }
518 }
519 PetscFunctionReturn(PETSC_SUCCESS);
520 #endif
521 }
522
KSPCGSetUp_STCG(KSP ksp)523 static PetscErrorCode KSPCGSetUp_STCG(KSP ksp)
524 {
525 PetscFunctionBegin;
526 /***************************************************************************/
527 /* Set work vectors needed by conjugate gradient method and allocate */
528 /***************************************************************************/
529
530 PetscCall(KSPSetWorkVecs(ksp, 3));
531 PetscFunctionReturn(PETSC_SUCCESS);
532 }
533
KSPCGDestroy_STCG(KSP ksp)534 static PetscErrorCode KSPCGDestroy_STCG(KSP ksp)
535 {
536 PetscFunctionBegin;
537 /***************************************************************************/
538 /* Clear composed functions */
539 /***************************************************************************/
540
541 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
542 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
543 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
544
545 /***************************************************************************/
546 /* Destroy KSP object. */
547 /***************************************************************************/
548
549 PetscCall(KSPDestroyDefault(ksp));
550 PetscFunctionReturn(PETSC_SUCCESS);
551 }
552
KSPCGSetRadius_STCG(KSP ksp,PetscReal radius)553 static PetscErrorCode KSPCGSetRadius_STCG(KSP ksp, PetscReal radius)
554 {
555 KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
556
557 PetscFunctionBegin;
558 cg->radius = radius;
559 PetscFunctionReturn(PETSC_SUCCESS);
560 }
561
KSPCGGetNormD_STCG(KSP ksp,PetscReal * norm_d)562 static PetscErrorCode KSPCGGetNormD_STCG(KSP ksp, PetscReal *norm_d)
563 {
564 KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
565
566 PetscFunctionBegin;
567 *norm_d = cg->norm_d;
568 PetscFunctionReturn(PETSC_SUCCESS);
569 }
570
KSPCGGetObjFcn_STCG(KSP ksp,PetscReal * o_fcn)571 static PetscErrorCode KSPCGGetObjFcn_STCG(KSP ksp, PetscReal *o_fcn)
572 {
573 KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
574
575 PetscFunctionBegin;
576 *o_fcn = cg->o_fcn;
577 PetscFunctionReturn(PETSC_SUCCESS);
578 }
579
KSPCGSetFromOptions_STCG(KSP ksp,PetscOptionItems PetscOptionsObject)580 static PetscErrorCode KSPCGSetFromOptions_STCG(KSP ksp, PetscOptionItems PetscOptionsObject)
581 {
582 KSPCG_STCG *cg = (KSPCG_STCG *)ksp->data;
583
584 PetscFunctionBegin;
585 PetscOptionsHeadBegin(PetscOptionsObject, "KSPCG STCG options");
586 PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
587 PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, STCG_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
588 PetscOptionsHeadEnd();
589 PetscFunctionReturn(PETSC_SUCCESS);
590 }
591
592 /*MC
593 KSPSTCG - Code to run conjugate gradient method subject to a constraint on the solution norm for use in a trust region method
594 {cite}`steihaug:83`, {cite}`toint1981towards`
595
596 Options Database Key:
597 . -ksp_cg_radius <r> - Trust Region Radius
598
599 Level: developer
600
601 Notes:
602 This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
603
604 Use preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
605
606 $$
607 q(s) = g^T s + \frac{1}{2} s^T H s
608 $$
609
610 subject to the trust region constraint
611
612 $$
613 || s || \le \delta,
614 $$
615
616 where
617 .vb
618 delta is the trust region radius,
619 g is the gradient vector,
620 H is the Hessian approximation, and
621 .ve
622
623 `KSPConvergedReason` may include
624 + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
625 - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step,
626
627 The preconditioner supplied should be symmetric and positive definite.
628
629 .seealso: [](ch_ksp), `KSPCreate()`, `KSPCGSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPNASH`, `KSPGLTR`, `KSPQCG`
630 M*/
631
KSPCreate_STCG(KSP ksp)632 PETSC_EXTERN PetscErrorCode KSPCreate_STCG(KSP ksp)
633 {
634 KSPCG_STCG *cg;
635
636 PetscFunctionBegin;
637 PetscCall(PetscNew(&cg));
638
639 cg->radius = 0.0;
640 cg->dtype = STCG_UNPRECONDITIONED_DIRECTION;
641
642 ksp->data = (void *)cg;
643 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
644 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
645 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
646 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
647 PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
648
649 /***************************************************************************/
650 /* Sets the functions that are associated with this data structure */
651 /* (in C++ this is the same as defining virtual functions). */
652 /***************************************************************************/
653
654 ksp->ops->setup = KSPCGSetUp_STCG;
655 ksp->ops->solve = KSPCGSolve_STCG;
656 ksp->ops->destroy = KSPCGDestroy_STCG;
657 ksp->ops->setfromoptions = KSPCGSetFromOptions_STCG;
658 ksp->ops->buildsolution = KSPBuildSolutionDefault;
659 ksp->ops->buildresidual = KSPBuildResidualDefault;
660 ksp->ops->view = NULL;
661
662 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_STCG));
663 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_STCG));
664 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_STCG));
665 PetscFunctionReturn(PETSC_SUCCESS);
666 }
667