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