1 #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h> /*I "petscksp.h" I*/
2 #include <petscblaslapack.h>
3
4 #define GLTR_PRECONDITIONED_DIRECTION 0
5 #define GLTR_UNPRECONDITIONED_DIRECTION 1
6 #define GLTR_DIRECTION_TYPES 2
7
8 static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
9
10 /*@
11 KSPGLTRGetMinEig - Get minimum eigenvalue computed by `KSPGLTR`
12
13 Collective
14
15 Input Parameter:
16 . ksp - the iterative context
17
18 Output Parameter:
19 . e_min - the minimum eigenvalue
20
21 Level: advanced
22
23 .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetLambda()`
24 @*/
KSPGLTRGetMinEig(KSP ksp,PetscReal * e_min)25 PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
26 {
27 PetscFunctionBegin;
28 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
29 PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
30 PetscFunctionReturn(PETSC_SUCCESS);
31 }
32
33 /*@
34 KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`
35
36 Not Collective
37
38 Input Parameter:
39 . ksp - the iterative context
40
41 Output Parameter:
42 . lambda - the multiplier
43
44 Level: advanced
45
46 .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetMinEig()`
47 @*/
KSPGLTRGetLambda(KSP ksp,PetscReal * lambda)48 PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
49 {
50 PetscFunctionBegin;
51 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
52 PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
53 PetscFunctionReturn(PETSC_SUCCESS);
54 }
55
KSPCGSolve_GLTR(KSP ksp)56 static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
57 {
58 #if defined(PETSC_USE_COMPLEX)
59 SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "GLTR is not available for complex systems");
60 #else
61 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
62 PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
63 PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
64
65 Mat Qmat, Mmat;
66 Vec r, z, p, d;
67 PC pc;
68
69 PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
70 PetscReal alpha, beta, kappa, rz, rzm1;
71 PetscReal rr, r2, piv, step;
72 PetscReal vl, vu;
73 PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
74 PetscReal norm_t, norm_w, pert;
75
76 PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
77 PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
78 PetscBLASInt nrhs, nldb;
79
80 PetscBLASInt e_valus = 0, e_splts;
81 PetscBool diagonalscale;
82
83 PetscFunctionBegin;
84 /* Check the arguments and parameters. */
85 PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
86 PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
87 PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
88
89 /* Get the workspace vectors and initialize variables */
90 r2 = cg->radius * cg->radius;
91 r = ksp->work[0];
92 z = ksp->work[1];
93 p = ksp->work[2];
94 d = ksp->vec_sol;
95 pc = ksp->pc;
96
97 PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
98
99 PetscCall(VecGetSize(d, &max_cg_its));
100 max_cg_its = PetscMin(max_cg_its, ksp->max_it);
101 max_lanczos_its = cg->max_lanczos_its;
102 max_newton_its = cg->max_newton_its;
103 ksp->its = 0;
104
105 /* Initialize objective function direction, and minimum eigenvalue. */
106 cg->o_fcn = 0.0;
107
108 PetscCall(VecSet(d, 0.0)); /* d = 0 */
109 cg->norm_d = 0.0;
110
111 cg->e_min = 0.0;
112 cg->lambda = 0.0;
113
114 /*
115 The first phase of GLTR performs a standard conjugate gradient method,
116 but stores the values required for the Lanczos matrix. We switch to
117 the Lanczos when the conjugate gradient method breaks down. Check the
118 right-hand side for numerical problems. The check for not-a-number and
119 infinite values need be performed only once.
120 */
121 PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
122 PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
123 KSPCheckDot(ksp, rr);
124
125 /*
126 Check the preconditioner for numerical problems and for positive
127 definiteness. The check for not-a-number and infinite values need be
128 performed only once.
129 */
130 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
131 PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
132 if (PetscIsInfOrNanScalar(rz)) {
133 /*
134 The preconditioner contains not-a-number or an infinite value.
135 Return the gradient direction intersected with the trust region.
136 */
137 ksp->reason = KSP_DIVERGED_NANORINF;
138 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));
139
140 if (cg->radius) {
141 if (r2 >= rr) {
142 alpha = 1.0;
143 cg->norm_d = PetscSqrtReal(rr);
144 } else {
145 alpha = PetscSqrtReal(r2 / rr);
146 cg->norm_d = cg->radius;
147 }
148
149 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
150
151 /* Compute objective function. */
152 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
153 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
154 PetscCall(VecDot(d, z, &cg->o_fcn));
155 cg->o_fcn = -cg->o_fcn;
156 ++ksp->its;
157 }
158 PetscFunctionReturn(PETSC_SUCCESS);
159 }
160
161 if (rz < 0.0) {
162 /*
163 The preconditioner is indefinite. Because this is the first
164 and we do not have a direction yet, we use the gradient step. Note
165 that we cannot use the preconditioned norm when computing the step
166 because the matrix is indefinite.
167 */
168 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
169 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));
170
171 if (cg->radius) {
172 if (r2 >= rr) {
173 alpha = 1.0;
174 cg->norm_d = PetscSqrtReal(rr);
175 } else {
176 alpha = PetscSqrtReal(r2 / rr);
177 cg->norm_d = cg->radius;
178 }
179
180 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
181
182 /* Compute objective function. */
183 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
184 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
185 PetscCall(VecDot(d, z, &cg->o_fcn));
186 cg->o_fcn = -cg->o_fcn;
187 ++ksp->its;
188 }
189 PetscFunctionReturn(PETSC_SUCCESS);
190 }
191
192 /*
193 As far as we know, the preconditioner is positive semidefinite.
194 Compute and log the residual. Check convergence because this
195 initializes things, but do not terminate until at least one conjugate
196 gradient iteration has been performed.
197 */
198 cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
199
200 switch (ksp->normtype) {
201 case KSP_NORM_PRECONDITIONED:
202 PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
203 break;
204
205 case KSP_NORM_UNPRECONDITIONED:
206 norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
207 break;
208
209 case KSP_NORM_NATURAL:
210 norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
211 break;
212
213 default:
214 norm_r = 0.0;
215 break;
216 }
217
218 PetscCall(KSPLogResidualHistory(ksp, norm_r));
219 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
220 ksp->rnorm = norm_r;
221
222 PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
223
224 /* Compute the first direction and update the iteration. */
225 PetscCall(VecCopy(z, p)); /* p = z */
226 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
227 ++ksp->its;
228
229 /* Check the matrix for numerical problems. */
230 PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
231 if (PetscIsInfOrNanScalar(kappa)) {
232 /*
233 The matrix produced not-a-number or an infinite value. In this case
234 we must stop and use the gradient direction. This condition need
235 only be checked once.
236 */
237 ksp->reason = KSP_DIVERGED_NANORINF;
238 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));
239
240 if (cg->radius) {
241 if (r2 >= rr) {
242 alpha = 1.0;
243 cg->norm_d = PetscSqrtReal(rr);
244 } else {
245 alpha = PetscSqrtReal(r2 / rr);
246 cg->norm_d = cg->radius;
247 }
248
249 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
250
251 /* Compute objective function. */
252 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
253 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
254 PetscCall(VecDot(d, z, &cg->o_fcn));
255 cg->o_fcn = -cg->o_fcn;
256 ++ksp->its;
257 }
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
261 /*
262 Initialize variables for calculating the norm of the direction and for
263 the Lanczos tridiagonal matrix. Note that we track the diagonal value
264 of the Cholesky factorization of the Lanczos matrix in order to
265 determine when negative curvature is encountered.
266 */
267 dMp = 0.0;
268 norm_d = 0.0;
269 switch (cg->dtype) {
270 case GLTR_PRECONDITIONED_DIRECTION:
271 norm_p = rz;
272 break;
273
274 default:
275 PetscCall(VecDot(p, p, &norm_p));
276 break;
277 }
278
279 cg->diag[t_size] = kappa / rz;
280 cg->offd[t_size] = 0.0;
281 ++t_size;
282
283 piv = 1.0;
284
285 /*
286 Check for breakdown of the conjugate gradient method; this occurs when
287 kappa is zero.
288 */
289 if (PetscAbsReal(kappa) <= 0.0) {
290 /* The curvature is zero. In this case, we must stop and use follow
291 the direction of negative curvature since the Lanczos matrix is zero. */
292 ksp->reason = KSP_DIVERGED_BREAKDOWN;
293 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));
294
295 if (cg->radius && norm_p > 0.0) {
296 /* Follow direction of negative curvature to the boundary of the
297 trust region. */
298 step = PetscSqrtReal(r2 / norm_p);
299 cg->norm_d = cg->radius;
300
301 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
302
303 /* Update objective function. */
304 cg->o_fcn += step * (0.5 * step * kappa - rz);
305 } else if (cg->radius) {
306 /* The norm of the preconditioned direction is zero; use the gradient
307 step. */
308 if (r2 >= rr) {
309 alpha = 1.0;
310 cg->norm_d = PetscSqrtReal(rr);
311 } else {
312 alpha = PetscSqrtReal(r2 / rr);
313 cg->norm_d = cg->radius;
314 }
315
316 PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
317
318 /* Compute objective function. */
319 PetscCall(KSP_MatMult(ksp, Qmat, d, z));
320 PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
321 PetscCall(VecDot(d, z, &cg->o_fcn));
322 cg->o_fcn = -cg->o_fcn;
323 ++ksp->its;
324 }
325 PetscFunctionReturn(PETSC_SUCCESS);
326 }
327
328 /*
329 Begin the first part of the GLTR algorithm which runs the conjugate
330 gradient method until either the problem is solved, we encounter the
331 boundary of the trust region, or the conjugate gradient method breaks
332 down.
333 */
334 while (1) {
335 /* Know that kappa is nonzero, because we have not broken down, so we */
336 /* can compute the steplength. */
337 alpha = rz / kappa;
338 cg->alpha[l_size] = alpha;
339
340 /* Compute the diagonal value of the Cholesky factorization of the */
341 /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
342 /* This indicates a direction of negative curvature. */
343 piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
344 if (piv <= 0.0) {
345 /* In this case, the matrix is indefinite and we have encountered */
346 /* a direction of negative curvature. Follow the direction to the */
347 /* boundary of the trust region. */
348 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
349 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));
350
351 if (cg->radius && norm_p > 0.0) {
352 /* Follow direction of negative curvature to boundary. */
353 step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
354 cg->norm_d = cg->radius;
355
356 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
357
358 /* Update objective function. */
359 cg->o_fcn += step * (0.5 * step * kappa - rz);
360 } else if (cg->radius) {
361 /* The norm of the direction is zero; there is nothing to follow. */
362 }
363 break;
364 }
365
366 /* Compute the steplength and check for intersection with the trust */
367 /* region. */
368 norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
369 if (cg->radius && norm_dp1 >= r2) {
370 /* In this case, the matrix is positive definite as far as we know. */
371 /* However, the full step goes beyond the trust region. */
372 ksp->reason = KSP_CONVERGED_STEP_LENGTH;
373 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));
374
375 if (norm_p > 0.0) {
376 /* Follow the direction to the boundary of the trust region. */
377
378 step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
379 cg->norm_d = cg->radius;
380
381 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
382
383 /* Update objective function. */
384 cg->o_fcn += step * (0.5 * step * kappa - rz);
385 } else {
386 /* The norm of the direction is zero; there is nothing to follow. */
387 }
388 break;
389 }
390
391 /* Now we can update the direction and residual. */
392 PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
393 PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
394 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
395
396 switch (cg->dtype) {
397 case GLTR_PRECONDITIONED_DIRECTION:
398 norm_d = norm_dp1;
399 break;
400
401 default:
402 PetscCall(VecDot(d, d, &norm_d));
403 break;
404 }
405 cg->norm_d = PetscSqrtReal(norm_d);
406
407 /* Update objective function. */
408 cg->o_fcn -= 0.5 * alpha * rz;
409
410 /* Check that the preconditioner appears positive semidefinite. */
411 rzm1 = rz;
412 PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
413 if (rz < 0.0) {
414 /* The preconditioner is indefinite. */
415 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
416 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
417 break;
418 }
419
420 /* As far as we know, the preconditioner is positive semidefinite. */
421 /* Compute the residual and check for convergence. */
422 cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
423
424 switch (ksp->normtype) {
425 case KSP_NORM_PRECONDITIONED:
426 PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
427 break;
428
429 case KSP_NORM_UNPRECONDITIONED:
430 PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
431 break;
432
433 case KSP_NORM_NATURAL:
434 norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
435 break;
436
437 default:
438 norm_r = 0.0;
439 break;
440 }
441
442 PetscCall(KSPLogResidualHistory(ksp, norm_r));
443 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
444 ksp->rnorm = norm_r;
445
446 PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
447 if (ksp->reason) {
448 /* The method has converged. */
449 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
450 break;
451 }
452
453 /* We have not converged yet. Check for breakdown. */
454 beta = rz / rzm1;
455 if (PetscAbsReal(beta) <= 0.0) {
456 /* Conjugate gradients has broken down. */
457 ksp->reason = KSP_DIVERGED_BREAKDOWN;
458 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
459 break;
460 }
461
462 /* Check iteration limit. */
463 if (ksp->its >= max_cg_its) {
464 ksp->reason = KSP_DIVERGED_ITS;
465 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
466 break;
467 }
468
469 /* Update p and the norms. */
470 cg->beta[l_size] = beta;
471 PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
472
473 switch (cg->dtype) {
474 case GLTR_PRECONDITIONED_DIRECTION:
475 dMp = beta * (dMp + alpha * norm_p);
476 norm_p = beta * (rzm1 + beta * norm_p);
477 break;
478
479 default:
480 PetscCall(VecDot(d, p, &dMp));
481 PetscCall(VecDot(p, p, &norm_p));
482 break;
483 }
484
485 /* Compute the new direction and update the iteration. */
486 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
487 PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
488 ++ksp->its;
489
490 /* Update the Lanczos tridiagonal matrix. */
491 ++l_size;
492 cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
493 cg->diag[t_size] = kappa / rz + beta / alpha;
494 ++t_size;
495
496 /* Check for breakdown of the conjugate gradient method; this occurs */
497 /* when kappa is zero. */
498 if (PetscAbsReal(kappa) <= 0.0) {
499 /* The method breaks down; move along the direction as if the matrix */
500 /* were indefinite. */
501 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
502 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));
503
504 if (cg->radius && norm_p > 0.0) {
505 /* Follow direction to boundary. */
506 step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
507 cg->norm_d = cg->radius;
508
509 PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
510
511 /* Update objective function. */
512 cg->o_fcn += step * (0.5 * step * kappa - rz);
513 } else if (cg->radius) {
514 /* The norm of the direction is zero; there is nothing to follow. */
515 }
516 break;
517 }
518 }
519
520 /* Check to see if we need to continue with the Lanczos method. */
521 if (!cg->radius) {
522 /* There is no radius. Therefore, we cannot move along the boundary. */
523 PetscFunctionReturn(PETSC_SUCCESS);
524 }
525
526 if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
527 /* The method either converged to an interior point, hit the boundary of */
528 /* the trust-region without encountering a direction of negative */
529 /* curvature or the iteration limit was reached. */
530 PetscFunctionReturn(PETSC_SUCCESS);
531 }
532
533 /* Switch to constructing the Lanczos basis by way of the conjugate */
534 /* directions. */
535 for (i = 0; i < max_lanczos_its; ++i) {
536 /* Check for breakdown of the conjugate gradient method; this occurs */
537 /* when kappa is zero. */
538 if (PetscAbsReal(kappa) <= 0.0) {
539 ksp->reason = KSP_DIVERGED_BREAKDOWN;
540 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
541 break;
542 }
543
544 /* Update the direction and residual. */
545 alpha = rz / kappa;
546 cg->alpha[l_size] = alpha;
547
548 PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
549 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
550
551 /* Check that the preconditioner appears positive semidefinite. */
552 rzm1 = rz;
553 PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
554 if (rz < 0.0) {
555 /* The preconditioner is indefinite. */
556 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
557 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
558 break;
559 }
560
561 /* As far as we know, the preconditioner is positive definite. Compute */
562 /* the residual. Do NOT check for convergence. */
563 cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
564
565 switch (ksp->normtype) {
566 case KSP_NORM_PRECONDITIONED:
567 PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
568 break;
569
570 case KSP_NORM_UNPRECONDITIONED:
571 PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
572 break;
573
574 case KSP_NORM_NATURAL:
575 norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
576 break;
577
578 default:
579 norm_r = 0.0;
580 break;
581 }
582
583 PetscCall(KSPLogResidualHistory(ksp, norm_r));
584 PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
585 ksp->rnorm = norm_r;
586
587 /* Check for breakdown. */
588 beta = rz / rzm1;
589 if (PetscAbsReal(beta) <= 0.0) {
590 /* Conjugate gradients has broken down. */
591 ksp->reason = KSP_DIVERGED_BREAKDOWN;
592 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
593 break;
594 }
595
596 /* Update p and the norms. */
597 cg->beta[l_size] = beta;
598 PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
599
600 /* Compute the new direction and update the iteration. */
601 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
602 PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
603 ++ksp->its;
604
605 /* Update the Lanczos tridiagonal matrix. */
606 ++l_size;
607 cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
608 cg->diag[t_size] = kappa / rz + beta / alpha;
609 ++t_size;
610 }
611
612 /*
613 We have the Lanczos basis, solve the tridiagonal trust-region problem
614 to obtain the Lanczos direction. We know that the solution lies on
615 the boundary of the trust region. We start by checking that the
616 workspace allocated is large enough.
617
618 Note that the current version only computes the solution by using the
619 preconditioned direction. Need to think about how to do the
620 unpreconditioned direction calculation.
621 */
622
623 if (t_size > cg->alloced) {
624 if (cg->alloced) {
625 PetscCall(PetscFree(cg->rwork));
626 PetscCall(PetscFree(cg->iwork));
627 cg->alloced += cg->init_alloc;
628 } else {
629 cg->alloced = cg->init_alloc;
630 }
631
632 while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
633
634 cg->alloced = PetscMin(cg->alloced, t_size);
635 PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
636 }
637
638 /* Set up the required vectors. */
639 t_soln = cg->rwork + 0 * t_size; /* Solution */
640 t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
641 t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
642 e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
643 e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
644 e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
645
646 e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
647 e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
648 e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
649
650 /* Compute the minimum eigenvalue of T. */
651 vl = 0.0;
652 vu = 0.0;
653 il = 1;
654 iu = 1;
655
656 PetscCallBLAS("LAPACKstebz", LAPACKstebz_("I", "E", &t_size, &vl, &vu, &il, &iu, &cg->eigen_tol, cg->diag, cg->offd + 1, &e_valus, &e_splts, e_valu, e_iblk, e_splt, e_rwrk, e_iwrk, &info));
657
658 if ((0 != info) || (1 != e_valus)) {
659 /* Calculation of the minimum eigenvalue failed. Return the */
660 /* Steihaug-Toint direction. */
661 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
662 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
663 PetscFunctionReturn(PETSC_SUCCESS);
664 }
665
666 cg->e_min = e_valu[0];
667
668 /* Compute the initial value of lambda to make (T + lambda I) positive */
669 /* definite. */
670 pert = cg->init_pert;
671 if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
672
673 while (1) {
674 for (i = 0; i < t_size; ++i) {
675 t_diag[i] = cg->diag[i] + cg->lambda;
676 t_offd[i] = cg->offd[i];
677 }
678
679 PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
680
681 if (0 == info) break;
682
683 pert += pert;
684 cg->lambda = cg->lambda * (1.0 + pert) + pert;
685 }
686
687 /* Compute the initial step and its norm. */
688 nrhs = 1;
689 nldb = t_size;
690
691 t_soln[0] = -cg->norm_r[0];
692 for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
693
694 PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
695
696 if (0 != info) {
697 /* Calculation of the initial step failed; return the Steihaug-Toint */
698 /* direction. */
699 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
700 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
701 PetscFunctionReturn(PETSC_SUCCESS);
702 }
703
704 norm_t = 0.;
705 for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
706 norm_t = PetscSqrtReal(norm_t);
707
708 /* Determine the case we are in. */
709 if (norm_t <= cg->radius) {
710 /* The step is within the trust region; check if we are in the hard case */
711 /* and need to move to the boundary by following a direction of negative */
712 /* curvature. */
713 if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
714 /* This is the hard case; compute the eigenvector associated with the */
715 /* minimum eigenvalue and move along this direction to the boundary. */
716 PetscCallBLAS("LAPACKstein", LAPACKstein_(&t_size, cg->diag, cg->offd + 1, &e_valus, e_valu, e_iblk, e_splt, e_vect, &nldb, e_rwrk, e_iwrk, e_iwrk + t_size, &info));
717
718 if (0 != info) {
719 /* Calculation of the minimum eigenvalue failed. Return the */
720 /* Steihaug-Toint direction. */
721 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
722 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
723 PetscFunctionReturn(PETSC_SUCCESS);
724 }
725
726 coef1 = 0.0;
727 coef2 = 0.0;
728 coef3 = -cg->radius * cg->radius;
729 for (i = 0; i < t_size; ++i) {
730 coef1 += e_vect[i] * e_vect[i];
731 coef2 += e_vect[i] * t_soln[i];
732 coef3 += t_soln[i] * t_soln[i];
733 }
734
735 coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
736 root1 = (-coef2 + coef3) / coef1;
737 root2 = (-coef2 - coef3) / coef1;
738
739 /* Compute objective value for (t_soln + root1 * e_vect) */
740 for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
741
742 obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
743 for (i = 1; i < t_size - 1; ++i) obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
744 obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
745
746 /* Compute objective value for (t_soln + root2 * e_vect) */
747 for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
748
749 obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
750 for (i = 1; i < t_size - 1; ++i) obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
751 obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
752
753 /* Choose the point with the best objective function value. */
754 if (obj1 <= obj2) {
755 for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
756 } else {
757 for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
758 }
759 } else {
760 /* The matrix is positive definite or there was no room to move; the */
761 /* solution is already contained in t_soln. */
762 }
763 } else {
764 /* The step is outside the trust-region. Compute the correct value for */
765 /* lambda by performing Newton's method. */
766
767 for (i = 0; i < max_newton_its; ++i) {
768 /* Check for convergence. */
769 if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
770
771 /* Compute the update. */
772 PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
773
774 PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
775
776 if (0 != info) {
777 /* Calculation of the step failed; return the Steihaug-Toint */
778 /* direction. */
779 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
780 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
781 PetscFunctionReturn(PETSC_SUCCESS);
782 }
783
784 /* Modify lambda. */
785 norm_w = 0.;
786 for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
787
788 cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
789
790 /* Factor T + lambda I */
791 for (j = 0; j < t_size; ++j) {
792 t_diag[j] = cg->diag[j] + cg->lambda;
793 t_offd[j] = cg->offd[j];
794 }
795
796 PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
797
798 if (0 != info) {
799 /* Calculation of factorization failed; return the Steihaug-Toint */
800 /* direction. */
801 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
802 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
803 PetscFunctionReturn(PETSC_SUCCESS);
804 }
805
806 /* Compute the new step and its norm. */
807 t_soln[0] = -cg->norm_r[0];
808 for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
809
810 PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
811
812 if (0 != info) {
813 /* Calculation of the step failed; return the Steihaug-Toint */
814 /* direction. */
815 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
816 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
817 PetscFunctionReturn(PETSC_SUCCESS);
818 }
819
820 norm_t = 0.;
821 for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
822 norm_t = PetscSqrtReal(norm_t);
823 }
824
825 /* Check for convergence. */
826 if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
827 /* Newton method failed to converge in iteration limit. */
828 PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
829 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
830 PetscFunctionReturn(PETSC_SUCCESS);
831 }
832 }
833
834 /* Recover the norm of the direction and objective function value. */
835 cg->norm_d = norm_t;
836
837 cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
838 for (i = 1; i < t_size - 1; ++i) cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i] + cg->offd[i + 1] * t_soln[i + 1]);
839 cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
840
841 /* Recover the direction. */
842 sigma = -1;
843
844 /* Start conjugate gradient method from the beginning */
845 PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
846 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
847
848 /* Accumulate Q * s */
849 PetscCall(VecCopy(z, d));
850 PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
851
852 /* Compute the first direction. */
853 PetscCall(VecCopy(z, p)); /* p = z */
854 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
855 ++ksp->its;
856
857 for (i = 0; i < l_size - 1; ++i) {
858 /* Update the residual and direction. */
859 alpha = cg->alpha[i];
860 if (alpha >= 0.0) sigma = -sigma;
861
862 PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
863 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
864
865 /* Accumulate Q * s */
866 PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
867
868 /* Update p. */
869 beta = cg->beta[i];
870 PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
871 PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
872 ++ksp->its;
873 }
874
875 /* Update the residual and direction. */
876 alpha = cg->alpha[i];
877 if (alpha >= 0.0) sigma = -sigma;
878
879 PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
880 PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
881
882 /* Accumulate Q * s */
883 PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
884
885 /* Set the termination reason. */
886 ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
887 PetscFunctionReturn(PETSC_SUCCESS);
888 #endif
889 }
890
KSPCGSetUp_GLTR(KSP ksp)891 static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
892 {
893 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
894 PetscInt max_its;
895
896 PetscFunctionBegin;
897 /* Determine the total maximum number of iterations. */
898 max_its = ksp->max_it + cg->max_lanczos_its + 1;
899
900 /* Set work vectors needed by conjugate gradient method and allocate */
901 /* workspace for Lanczos matrix. */
902 PetscCall(KSPSetWorkVecs(ksp, 3));
903 if (cg->diag) {
904 PetscCall(PetscArrayzero(cg->diag, max_its));
905 PetscCall(PetscArrayzero(cg->offd, max_its));
906 PetscCall(PetscArrayzero(cg->alpha, max_its));
907 PetscCall(PetscArrayzero(cg->beta, max_its));
908 PetscCall(PetscArrayzero(cg->norm_r, max_its));
909 } else {
910 PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
911 }
912 PetscFunctionReturn(PETSC_SUCCESS);
913 }
914
KSPCGDestroy_GLTR(KSP ksp)915 static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
916 {
917 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
918
919 PetscFunctionBegin;
920 PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
921 if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
922 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
923 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
924 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
925 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
926 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
927 PetscCall(KSPDestroyDefault(ksp));
928 PetscFunctionReturn(PETSC_SUCCESS);
929 }
930
KSPCGSetRadius_GLTR(KSP ksp,PetscReal radius)931 static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
932 {
933 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
934
935 PetscFunctionBegin;
936 cg->radius = radius;
937 PetscFunctionReturn(PETSC_SUCCESS);
938 }
939
KSPCGGetNormD_GLTR(KSP ksp,PetscReal * norm_d)940 static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
941 {
942 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
943
944 PetscFunctionBegin;
945 *norm_d = cg->norm_d;
946 PetscFunctionReturn(PETSC_SUCCESS);
947 }
948
KSPCGGetObjFcn_GLTR(KSP ksp,PetscReal * o_fcn)949 static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
950 {
951 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
952
953 PetscFunctionBegin;
954 *o_fcn = cg->o_fcn;
955 PetscFunctionReturn(PETSC_SUCCESS);
956 }
957
KSPGLTRGetMinEig_GLTR(KSP ksp,PetscReal * e_min)958 static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
959 {
960 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
961
962 PetscFunctionBegin;
963 *e_min = cg->e_min;
964 PetscFunctionReturn(PETSC_SUCCESS);
965 }
966
KSPGLTRGetLambda_GLTR(KSP ksp,PetscReal * lambda)967 static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
968 {
969 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
970
971 PetscFunctionBegin;
972 *lambda = cg->lambda;
973 PetscFunctionReturn(PETSC_SUCCESS);
974 }
975
KSPCGSetFromOptions_GLTR(KSP ksp,PetscOptionItems PetscOptionsObject)976 static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems PetscOptionsObject)
977 {
978 KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
979
980 PetscFunctionBegin;
981 PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
982
983 PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
984
985 PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
986
987 PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
988 PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
989 PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
990
991 PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
992 PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
993
994 PetscOptionsHeadEnd();
995 PetscFunctionReturn(PETSC_SUCCESS);
996 }
997
998 /*MC
999 KSPGLTR - Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`
1000
1001 Options Database Key:
1002 . -ksp_cg_radius <r> - Trust Region Radius
1003
1004 Level: developer
1005
1006 Notes:
1007 Uses preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
1008
1009 $$
1010 q(s) = g^T * s + .5 * s^T * H * s
1011 $$
1012
1013 subject to the trust region constraint
1014
1015 $$
1016 || s || \le delta,
1017 $$
1018
1019 where
1020 .vb
1021 delta is the trust region radius,
1022 g is the gradient vector,
1023 H is the Hessian approximation,
1024 .ve
1025
1026 `KSPConvergedReason` may have the additional values
1027 + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
1028 - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.
1029
1030 The operator and the preconditioner supplied must be symmetric and positive definite.
1031
1032 This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1033
1034 .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1035 M*/
1036
KSPCreate_GLTR(KSP ksp)1037 PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1038 {
1039 KSPCG_GLTR *cg;
1040
1041 PetscFunctionBegin;
1042 PetscCall(PetscNew(&cg));
1043 cg->radius = 0.0;
1044 cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1045
1046 cg->init_pert = 1.0e-8;
1047 cg->eigen_tol = 1.0e-10;
1048 cg->newton_tol = 1.0e-6;
1049
1050 cg->alloced = 0;
1051 cg->init_alloc = 1024;
1052
1053 cg->max_lanczos_its = 20;
1054 cg->max_newton_its = 10;
1055
1056 ksp->data = (void *)cg;
1057 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1058 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1059 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1060 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1061 PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1062
1063 /* Sets the functions that are associated with this data structure */
1064 /* (in C++ this is the same as defining virtual functions). */
1065
1066 ksp->ops->setup = KSPCGSetUp_GLTR;
1067 ksp->ops->solve = KSPCGSolve_GLTR;
1068 ksp->ops->destroy = KSPCGDestroy_GLTR;
1069 ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1070 ksp->ops->buildsolution = KSPBuildSolutionDefault;
1071 ksp->ops->buildresidual = KSPBuildResidualDefault;
1072 ksp->ops->view = NULL;
1073
1074 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1075 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1076 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1077 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1078 PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1079 PetscFunctionReturn(PETSC_SUCCESS);
1080 }
1081