1
2 #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h> /*I "petscksp.h" I*/
3 #include <petscblaslapack.h>
4
5 #define GLTR_PRECONDITIONED_DIRECTION 0
6 #define GLTR_UNPRECONDITIONED_DIRECTION 1
7 #define GLTR_DIRECTION_TYPES 2
8
9 static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
10
11 /*@
12 KSPGLTRGetMinEig - Get minimum eigenvalue.
13
14 Collective on ksp
15
16 Input Parameters:
17 + ksp - the iterative context
18 - e_min - the minimum eigenvalue
19
20 Level: advanced
21
22 @*/
KSPGLTRGetMinEig(KSP ksp,PetscReal * e_min)23 PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
24 {
25 PetscErrorCode ierr;
26
27 PetscFunctionBegin;
28 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
29 ierr = PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));CHKERRQ(ierr);
30 PetscFunctionReturn(0);
31 }
32
33 /*@
34 KSPGLTRGetLambda - Get multiplier on trust-region constraint.
35
36 Not Collective
37
38 Input Parameters:
39 + ksp - the iterative context
40 - lambda - the multiplier
41
42 Level: advanced
43
44 @*/
KSPGLTRGetLambda(KSP ksp,PetscReal * lambda)45 PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
46 {
47 PetscErrorCode ierr;
48
49 PetscFunctionBegin;
50 PetscValidHeaderSpecific(ksp, KSP_CLASSID, 1);
51 ierr = PetscUseMethod(ksp,"KSPGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));CHKERRQ(ierr);
52 PetscFunctionReturn(0);
53 }
54
KSPCGSolve_GLTR(KSP ksp)55 static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
56 {
57 #if defined(PETSC_USE_COMPLEX)
58 SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "GLTR is not available for complex systems");
59 #else
60 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
61 PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
62 PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
63
64 PetscErrorCode ierr;
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 /***************************************************************************/
85 /* Check the arguments and parameters. */
86 /***************************************************************************/
87
88 ierr = PCGetDiagonalScale(ksp->pc, &diagonalscale);CHKERRQ(ierr);
89 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
90 if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
91
92 /***************************************************************************/
93 /* Get the workspace vectors and initialize variables */
94 /***************************************************************************/
95
96 r2 = cg->radius * cg->radius;
97 r = ksp->work[0];
98 z = ksp->work[1];
99 p = ksp->work[2];
100 d = ksp->vec_sol;
101 pc = ksp->pc;
102
103 ierr = PCGetOperators(pc, &Qmat, &Mmat);CHKERRQ(ierr);
104
105 ierr = VecGetSize(d, &max_cg_its);CHKERRQ(ierr);
106 max_cg_its = PetscMin(max_cg_its, ksp->max_it);
107 max_lanczos_its = cg->max_lanczos_its;
108 max_newton_its = cg->max_newton_its;
109 ksp->its = 0;
110
111 /***************************************************************************/
112 /* Initialize objective function direction, and minimum eigenvalue. */
113 /***************************************************************************/
114
115 cg->o_fcn = 0.0;
116
117 ierr = VecSet(d, 0.0);CHKERRQ(ierr); /* d = 0 */
118 cg->norm_d = 0.0;
119
120 cg->e_min = 0.0;
121 cg->lambda = 0.0;
122
123 /***************************************************************************/
124 /* The first phase of GLTR performs a standard conjugate gradient method, */
125 /* but stores the values required for the Lanczos matrix. We switch to */
126 /* the Lanczos when the conjugate gradient method breaks down. Check the */
127 /* right-hand side for numerical problems. The check for not-a-number and */
128 /* infinite values need be performed only once. */
129 /***************************************************************************/
130
131 ierr = VecCopy(ksp->vec_rhs, r);CHKERRQ(ierr); /* r = -grad */
132 ierr = VecDot(r, r, &rr);CHKERRQ(ierr); /* rr = r^T r */
133 KSPCheckDot(ksp,rr);
134
135 /***************************************************************************/
136 /* Check the preconditioner for numerical problems and for positive */
137 /* definiteness. The check for not-a-number and infinite values need be */
138 /* performed only once. */
139 /***************************************************************************/
140
141 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
142 ierr = VecDot(r, z, &rz);CHKERRQ(ierr); /* rz = r^T inv(M) r */
143 if (PetscIsInfOrNanScalar(rz)) {
144 /*************************************************************************/
145 /* The preconditioner contains not-a-number or an infinite value. */
146 /* Return the gradient direction intersected with the trust region. */
147 /*************************************************************************/
148
149 ksp->reason = KSP_DIVERGED_NANORINF;
150 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);CHKERRQ(ierr);
151
152 if (cg->radius) {
153 if (r2 >= rr) {
154 alpha = 1.0;
155 cg->norm_d = PetscSqrtReal(rr);
156 } else {
157 alpha = PetscSqrtReal(r2 / rr);
158 cg->norm_d = cg->radius;
159 }
160
161 ierr = VecAXPY(d, alpha, r);CHKERRQ(ierr); /* d = d + alpha r */
162
163 /***********************************************************************/
164 /* Compute objective function. */
165 /***********************************************************************/
166
167 ierr = KSP_MatMult(ksp, Qmat, d, z);CHKERRQ(ierr);
168 ierr = VecAYPX(z, -0.5, ksp->vec_rhs);CHKERRQ(ierr);
169 ierr = VecDot(d, z, &cg->o_fcn);CHKERRQ(ierr);
170 cg->o_fcn = -cg->o_fcn;
171 ++ksp->its;
172 }
173 PetscFunctionReturn(0);
174 }
175
176 if (rz < 0.0) {
177 /*************************************************************************/
178 /* The preconditioner is indefinite. Because this is the first */
179 /* and we do not have a direction yet, we use the gradient step. Note */
180 /* that we cannot use the preconditioned norm when computing the step */
181 /* because the matrix is indefinite. */
182 /*************************************************************************/
183
184 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
185 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);CHKERRQ(ierr);
186
187 if (cg->radius) {
188 if (r2 >= rr) {
189 alpha = 1.0;
190 cg->norm_d = PetscSqrtReal(rr);
191 } else {
192 alpha = PetscSqrtReal(r2 / rr);
193 cg->norm_d = cg->radius;
194 }
195
196 ierr = VecAXPY(d, alpha, r);CHKERRQ(ierr); /* d = d + alpha r */
197
198 /***********************************************************************/
199 /* Compute objective function. */
200 /***********************************************************************/
201
202 ierr = KSP_MatMult(ksp, Qmat, d, z);CHKERRQ(ierr);
203 ierr = VecAYPX(z, -0.5, ksp->vec_rhs);CHKERRQ(ierr);
204 ierr = VecDot(d, z, &cg->o_fcn);CHKERRQ(ierr);
205 cg->o_fcn = -cg->o_fcn;
206 ++ksp->its;
207 }
208 PetscFunctionReturn(0);
209 }
210
211 /***************************************************************************/
212 /* As far as we know, the preconditioner is positive semidefinite. */
213 /* Compute and log the residual. Check convergence because this */
214 /* initializes things, but do not terminate until at least one conjugate */
215 /* gradient iteration has been performed. */
216 /***************************************************************************/
217
218 cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
219
220 switch (ksp->normtype) {
221 case KSP_NORM_PRECONDITIONED:
222 ierr = VecNorm(z, NORM_2, &norm_r);CHKERRQ(ierr); /* norm_r = |z| */
223 break;
224
225 case KSP_NORM_UNPRECONDITIONED:
226 norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
227 break;
228
229 case KSP_NORM_NATURAL:
230 norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
231 break;
232
233 default:
234 norm_r = 0.0;
235 break;
236 }
237
238 ierr = KSPLogResidualHistory(ksp, norm_r);CHKERRQ(ierr);
239 ierr = KSPMonitor(ksp, ksp->its, norm_r);CHKERRQ(ierr);
240 ksp->rnorm = norm_r;
241
242 ierr = (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);CHKERRQ(ierr);
243
244 /***************************************************************************/
245 /* Compute the first direction and update the iteration. */
246 /***************************************************************************/
247
248 ierr = VecCopy(z, p);CHKERRQ(ierr); /* p = z */
249 ierr = KSP_MatMult(ksp, Qmat, p, z);CHKERRQ(ierr); /* z = Q * p */
250 ++ksp->its;
251
252 /***************************************************************************/
253 /* Check the matrix for numerical problems. */
254 /***************************************************************************/
255
256 ierr = VecDot(p, z, &kappa);CHKERRQ(ierr); /* kappa = p^T Q p */
257 if (PetscIsInfOrNanScalar(kappa)) {
258 /*************************************************************************/
259 /* The matrix produced not-a-number or an infinite value. In this case, */
260 /* we must stop and use the gradient direction. This condition need */
261 /* only be checked once. */
262 /*************************************************************************/
263
264 ksp->reason = KSP_DIVERGED_NANORINF;
265 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);CHKERRQ(ierr);
266
267 if (cg->radius) {
268 if (r2 >= rr) {
269 alpha = 1.0;
270 cg->norm_d = PetscSqrtReal(rr);
271 } else {
272 alpha = PetscSqrtReal(r2 / rr);
273 cg->norm_d = cg->radius;
274 }
275
276 ierr = VecAXPY(d, alpha, r);CHKERRQ(ierr); /* d = d + alpha r */
277
278 /***********************************************************************/
279 /* Compute objective function. */
280 /***********************************************************************/
281
282 ierr = KSP_MatMult(ksp, Qmat, d, z);CHKERRQ(ierr);
283 ierr = VecAYPX(z, -0.5, ksp->vec_rhs);CHKERRQ(ierr);
284 ierr = VecDot(d, z, &cg->o_fcn);CHKERRQ(ierr);
285 cg->o_fcn = -cg->o_fcn;
286 ++ksp->its;
287 }
288 PetscFunctionReturn(0);
289 }
290
291 /***************************************************************************/
292 /* Initialize variables for calculating the norm of the direction and for */
293 /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
294 /* of the Cholesky factorization of the Lanczos matrix in order to */
295 /* determine when negative curvature is encountered. */
296 /***************************************************************************/
297
298 dMp = 0.0;
299 norm_d = 0.0;
300 switch (cg->dtype) {
301 case GLTR_PRECONDITIONED_DIRECTION:
302 norm_p = rz;
303 break;
304
305 default:
306 ierr = VecDot(p, p, &norm_p);CHKERRQ(ierr);
307 break;
308 }
309
310 cg->diag[t_size] = kappa / rz;
311 cg->offd[t_size] = 0.0;
312 ++t_size;
313
314 piv = 1.0;
315
316 /***************************************************************************/
317 /* Check for breakdown of the conjugate gradient method; this occurs when */
318 /* kappa is zero. */
319 /***************************************************************************/
320
321 if (PetscAbsReal(kappa) <= 0.0) {
322 /*************************************************************************/
323 /* The curvature is zero. In this case, we must stop and use follow */
324 /* the direction of negative curvature since the Lanczos matrix is zero. */
325 /*************************************************************************/
326
327 ksp->reason = KSP_DIVERGED_BREAKDOWN;
328 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);CHKERRQ(ierr);
329
330 if (cg->radius && norm_p > 0.0) {
331 /***********************************************************************/
332 /* Follow direction of negative curvature to the boundary of the */
333 /* trust region. */
334 /***********************************************************************/
335
336 step = PetscSqrtReal(r2 / norm_p);
337 cg->norm_d = cg->radius;
338
339 ierr = VecAXPY(d, step, p);CHKERRQ(ierr); /* d = d + step p */
340
341 /***********************************************************************/
342 /* Update objective function. */
343 /***********************************************************************/
344
345 cg->o_fcn += step * (0.5 * step * kappa - rz);
346 } else if (cg->radius) {
347 /***********************************************************************/
348 /* The norm of the preconditioned direction is zero; use the gradient */
349 /* step. */
350 /***********************************************************************/
351
352 if (r2 >= rr) {
353 alpha = 1.0;
354 cg->norm_d = PetscSqrtReal(rr);
355 } else {
356 alpha = PetscSqrtReal(r2 / rr);
357 cg->norm_d = cg->radius;
358 }
359
360 ierr = VecAXPY(d, alpha, r);CHKERRQ(ierr); /* d = d + alpha r */
361
362 /***********************************************************************/
363 /* Compute objective function. */
364 /***********************************************************************/
365
366 ierr = KSP_MatMult(ksp, Qmat, d, z);CHKERRQ(ierr);
367 ierr = VecAYPX(z, -0.5, ksp->vec_rhs);CHKERRQ(ierr);
368 ierr = VecDot(d, z, &cg->o_fcn);CHKERRQ(ierr);
369 cg->o_fcn = -cg->o_fcn;
370 ++ksp->its;
371 }
372 PetscFunctionReturn(0);
373 }
374
375 /***************************************************************************/
376 /* Begin the first part of the GLTR algorithm which runs the conjugate */
377 /* gradient method until either the problem is solved, we encounter the */
378 /* boundary of the trust region, or the conjugate gradient method breaks */
379 /* down. */
380 /***************************************************************************/
381
382 while (1) {
383 /*************************************************************************/
384 /* Know that kappa is nonzero, because we have not broken down, so we */
385 /* can compute the steplength. */
386 /*************************************************************************/
387
388 alpha = rz / kappa;
389 cg->alpha[l_size] = alpha;
390
391 /*************************************************************************/
392 /* Compute the diagonal value of the Cholesky factorization of the */
393 /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
394 /* This indicates a direction of negative curvature. */
395 /*************************************************************************/
396
397 piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
398 if (piv <= 0.0) {
399 /***********************************************************************/
400 /* In this case, the matrix is indefinite and we have encountered */
401 /* a direction of negative curvature. Follow the direction to the */
402 /* boundary of the trust region. */
403 /***********************************************************************/
404
405 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
406 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);CHKERRQ(ierr);
407
408 if (cg->radius && norm_p > 0.0) {
409 /*********************************************************************/
410 /* Follow direction of negative curvature to boundary. */
411 /*********************************************************************/
412
413 step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
414 cg->norm_d = cg->radius;
415
416 ierr = VecAXPY(d, step, p);CHKERRQ(ierr); /* d = d + step p */
417
418 /*********************************************************************/
419 /* Update objective function. */
420 /*********************************************************************/
421
422 cg->o_fcn += step * (0.5 * step * kappa - rz);
423 } else if (cg->radius) {
424 /*********************************************************************/
425 /* The norm of the direction is zero; there is nothing to follow. */
426 /*********************************************************************/
427 }
428 break;
429 }
430
431 /*************************************************************************/
432 /* Compute the steplength and check for intersection with the trust */
433 /* region. */
434 /*************************************************************************/
435
436 norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
437 if (cg->radius && norm_dp1 >= r2) {
438 /***********************************************************************/
439 /* In this case, the matrix is positive definite as far as we know. */
440 /* However, the full step goes beyond the trust region. */
441 /***********************************************************************/
442
443 ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
444 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);CHKERRQ(ierr);
445
446 if (norm_p > 0.0) {
447 /*********************************************************************/
448 /* Follow the direction to the boundary of the trust region. */
449 /*********************************************************************/
450
451 step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
452 cg->norm_d = cg->radius;
453
454 ierr = VecAXPY(d, step, p);CHKERRQ(ierr); /* d = d + step p */
455
456 /*********************************************************************/
457 /* Update objective function. */
458 /*********************************************************************/
459
460 cg->o_fcn += step * (0.5 * step * kappa - rz);
461 } else {
462 /*********************************************************************/
463 /* The norm of the direction is zero; there is nothing to follow. */
464 /*********************************************************************/
465 }
466 break;
467 }
468
469 /*************************************************************************/
470 /* Now we can update the direction and residual. */
471 /*************************************************************************/
472
473 ierr = VecAXPY(d, alpha, p);CHKERRQ(ierr); /* d = d + alpha p */
474 ierr = VecAXPY(r, -alpha, z);CHKERRQ(ierr); /* r = r - alpha Q p */
475 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
476
477 switch (cg->dtype) {
478 case GLTR_PRECONDITIONED_DIRECTION:
479 norm_d = norm_dp1;
480 break;
481
482 default:
483 ierr = VecDot(d, d, &norm_d);CHKERRQ(ierr);
484 break;
485 }
486 cg->norm_d = PetscSqrtReal(norm_d);
487
488 /*************************************************************************/
489 /* Update objective function. */
490 /*************************************************************************/
491
492 cg->o_fcn -= 0.5 * alpha * rz;
493
494 /*************************************************************************/
495 /* Check that the preconditioner appears positive semidefinite. */
496 /*************************************************************************/
497
498 rzm1 = rz;
499 ierr = VecDot(r, z, &rz);CHKERRQ(ierr); /* rz = r^T z */
500 if (rz < 0.0) {
501 /***********************************************************************/
502 /* The preconditioner is indefinite. */
503 /***********************************************************************/
504
505 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
506 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);CHKERRQ(ierr);
507 break;
508 }
509
510 /*************************************************************************/
511 /* As far as we know, the preconditioner is positive semidefinite. */
512 /* Compute the residual and check for convergence. */
513 /*************************************************************************/
514
515 cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
516
517 switch (ksp->normtype) {
518 case KSP_NORM_PRECONDITIONED:
519 ierr = VecNorm(z, NORM_2, &norm_r);CHKERRQ(ierr); /* norm_r = |z| */
520 break;
521
522 case KSP_NORM_UNPRECONDITIONED:
523 ierr = VecNorm(r, NORM_2, &norm_r);CHKERRQ(ierr); /* norm_r = |r| */
524 break;
525
526 case KSP_NORM_NATURAL:
527 norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
528 break;
529
530 default:
531 norm_r = 0.0;
532 break;
533 }
534
535 ierr = KSPLogResidualHistory(ksp, norm_r);CHKERRQ(ierr);
536 ierr = KSPMonitor(ksp, ksp->its, norm_r);CHKERRQ(ierr);
537 ksp->rnorm = norm_r;
538
539 ierr = (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);CHKERRQ(ierr);
540 if (ksp->reason) {
541 /***********************************************************************/
542 /* The method has converged. */
543 /***********************************************************************/
544
545 ierr = PetscInfo2(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);CHKERRQ(ierr);
546 break;
547 }
548
549 /*************************************************************************/
550 /* We have not converged yet. Check for breakdown. */
551 /*************************************************************************/
552
553 beta = rz / rzm1;
554 if (PetscAbsReal(beta) <= 0.0) {
555 /***********************************************************************/
556 /* Conjugate gradients has broken down. */
557 /***********************************************************************/
558
559 ksp->reason = KSP_DIVERGED_BREAKDOWN;
560 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);CHKERRQ(ierr);
561 break;
562 }
563
564 /*************************************************************************/
565 /* Check iteration limit. */
566 /*************************************************************************/
567
568 if (ksp->its >= max_cg_its) {
569 ksp->reason = KSP_DIVERGED_ITS;
570 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);CHKERRQ(ierr);
571 break;
572 }
573
574 /*************************************************************************/
575 /* Update p and the norms. */
576 /*************************************************************************/
577
578 cg->beta[l_size] = beta;
579 ierr = VecAYPX(p, beta, z);CHKERRQ(ierr); /* p = z + beta p */
580
581 switch (cg->dtype) {
582 case GLTR_PRECONDITIONED_DIRECTION:
583 dMp = beta*(dMp + alpha*norm_p);
584 norm_p = beta*(rzm1 + beta*norm_p);
585 break;
586
587 default:
588 ierr = VecDot(d, p, &dMp);CHKERRQ(ierr);
589 ierr = VecDot(p, p, &norm_p);CHKERRQ(ierr);
590 break;
591 }
592
593 /*************************************************************************/
594 /* Compute the new direction and update the iteration. */
595 /*************************************************************************/
596
597 ierr = KSP_MatMult(ksp, Qmat, p, z);CHKERRQ(ierr); /* z = Q * p */
598 ierr = VecDot(p, z, &kappa);CHKERRQ(ierr); /* kappa = p^T Q p */
599 ++ksp->its;
600
601 /*************************************************************************/
602 /* Update the Lanczos tridiagonal matrix. */
603 /*************************************************************************/
604
605 ++l_size;
606 cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
607 cg->diag[t_size] = kappa / rz + beta / alpha;
608 ++t_size;
609
610 /*************************************************************************/
611 /* Check for breakdown of the conjugate gradient method; this occurs */
612 /* when kappa is zero. */
613 /*************************************************************************/
614
615 if (PetscAbsReal(kappa) <= 0.0) {
616 /***********************************************************************/
617 /* The method breaks down; move along the direction as if the matrix */
618 /* were indefinite. */
619 /***********************************************************************/
620
621 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
622 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);CHKERRQ(ierr);
623
624 if (cg->radius && norm_p > 0.0) {
625 /*********************************************************************/
626 /* Follow direction to boundary. */
627 /*********************************************************************/
628
629 step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
630 cg->norm_d = cg->radius;
631
632 ierr = VecAXPY(d, step, p);CHKERRQ(ierr); /* d = d + step p */
633
634 /*********************************************************************/
635 /* Update objective function. */
636 /*********************************************************************/
637
638 cg->o_fcn += step * (0.5 * step * kappa - rz);
639 } else if (cg->radius) {
640 /*********************************************************************/
641 /* The norm of the direction is zero; there is nothing to follow. */
642 /*********************************************************************/
643 }
644 break;
645 }
646 }
647
648 /***************************************************************************/
649 /* Check to see if we need to continue with the Lanczos method. */
650 /***************************************************************************/
651
652 if (!cg->radius) {
653 /*************************************************************************/
654 /* There is no radius. Therefore, we cannot move along the boundary. */
655 /*************************************************************************/
656 PetscFunctionReturn(0);
657 }
658
659 if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
660 /*************************************************************************/
661 /* The method either converged to an interior point, hit the boundary of */
662 /* the trust-region without encountering a direction of negative */
663 /* curvature or the iteration limit was reached. */
664 /*************************************************************************/
665 PetscFunctionReturn(0);
666 }
667
668 /***************************************************************************/
669 /* Switch to contructing the Lanczos basis by way of the conjugate */
670 /* directions. */
671 /***************************************************************************/
672
673 for (i = 0; i < max_lanczos_its; ++i) {
674 /*************************************************************************/
675 /* Check for breakdown of the conjugate gradient method; this occurs */
676 /* when kappa is zero. */
677 /*************************************************************************/
678
679 if (PetscAbsReal(kappa) <= 0.0) {
680 ksp->reason = KSP_DIVERGED_BREAKDOWN;
681 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);CHKERRQ(ierr);
682 break;
683 }
684
685 /*************************************************************************/
686 /* Update the direction and residual. */
687 /*************************************************************************/
688
689 alpha = rz / kappa;
690 cg->alpha[l_size] = alpha;
691
692 ierr = VecAXPY(r, -alpha, z);CHKERRQ(ierr); /* r = r - alpha Q p */
693 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
694
695 /*************************************************************************/
696 /* Check that the preconditioner appears positive semidefinite. */
697 /*************************************************************************/
698
699 rzm1 = rz;
700 ierr = VecDot(r, z, &rz);CHKERRQ(ierr); /* rz = r^T z */
701 if (rz < 0.0) {
702 /***********************************************************************/
703 /* The preconditioner is indefinite. */
704 /***********************************************************************/
705
706 ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
707 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);CHKERRQ(ierr);
708 break;
709 }
710
711 /*************************************************************************/
712 /* As far as we know, the preconditioner is positive definite. Compute */
713 /* the residual. Do NOT check for convergence. */
714 /*************************************************************************/
715
716 cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
717
718 switch (ksp->normtype) {
719 case KSP_NORM_PRECONDITIONED:
720 ierr = VecNorm(z, NORM_2, &norm_r);CHKERRQ(ierr); /* norm_r = |z| */
721 break;
722
723 case KSP_NORM_UNPRECONDITIONED:
724 ierr = VecNorm(r, NORM_2, &norm_r);CHKERRQ(ierr); /* norm_r = |r| */
725 break;
726
727 case KSP_NORM_NATURAL:
728 norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
729 break;
730
731 default:
732 norm_r = 0.0;
733 break;
734 }
735
736 ierr = KSPLogResidualHistory(ksp, norm_r);CHKERRQ(ierr);
737 ierr = KSPMonitor(ksp, ksp->its, norm_r);CHKERRQ(ierr);
738 ksp->rnorm = norm_r;
739
740 /*************************************************************************/
741 /* Check for breakdown. */
742 /*************************************************************************/
743
744 beta = rz / rzm1;
745 if (PetscAbsReal(beta) <= 0.0) {
746 /***********************************************************************/
747 /* Conjugate gradients has broken down. */
748 /***********************************************************************/
749
750 ksp->reason = KSP_DIVERGED_BREAKDOWN;
751 ierr = PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n",(double) beta);CHKERRQ(ierr);
752 break;
753 }
754
755 /*************************************************************************/
756 /* Update p and the norms. */
757 /*************************************************************************/
758
759 cg->beta[l_size] = beta;
760 ierr = VecAYPX(p, beta, z);CHKERRQ(ierr); /* p = z + beta p */
761
762 /*************************************************************************/
763 /* Compute the new direction and update the iteration. */
764 /*************************************************************************/
765
766 ierr = KSP_MatMult(ksp, Qmat, p, z);CHKERRQ(ierr); /* z = Q * p */
767 ierr = VecDot(p, z, &kappa);CHKERRQ(ierr); /* kappa = p^T Q p */
768 ++ksp->its;
769
770 /*************************************************************************/
771 /* Update the Lanczos tridiagonal matrix. */
772 /*************************************************************************/
773
774 ++l_size;
775 cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
776 cg->diag[t_size] = kappa / rz + beta / alpha;
777 ++t_size;
778 }
779
780 /***************************************************************************/
781 /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
782 /* to obtain the Lanczos direction. We know that the solution lies on */
783 /* the boundary of the trust region. We start by checking that the */
784 /* workspace allocated is large enough. */
785 /***************************************************************************/
786 /* Note that the current version only computes the solution by using the */
787 /* preconditioned direction. Need to think about how to do the */
788 /* unpreconditioned direction calculation. */
789 /***************************************************************************/
790
791 if (t_size > cg->alloced) {
792 if (cg->alloced) {
793 ierr = PetscFree(cg->rwork);CHKERRQ(ierr);
794 ierr = PetscFree(cg->iwork);CHKERRQ(ierr);
795 cg->alloced += cg->init_alloc;
796 } else {
797 cg->alloced = cg->init_alloc;
798 }
799
800 while (t_size > cg->alloced) {
801 cg->alloced += cg->init_alloc;
802 }
803
804 cg->alloced = PetscMin(cg->alloced, t_size);
805 ierr = PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);CHKERRQ(ierr);
806 }
807
808 /***************************************************************************/
809 /* Set up the required vectors. */
810 /***************************************************************************/
811
812 t_soln = cg->rwork + 0*t_size; /* Solution */
813 t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
814 t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
815 e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
816 e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
817 e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
818
819 e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
820 e_splt = cg->iwork + 1*t_size; /* Eigen splits */
821 e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
822
823 /***************************************************************************/
824 /* Compute the minimum eigenvalue of T. */
825 /***************************************************************************/
826
827 vl = 0.0;
828 vu = 0.0;
829 il = 1;
830 iu = 1;
831
832 PetscStackCallBLAS("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));
833
834 if ((0 != info) || (1 != e_valus)) {
835 /*************************************************************************/
836 /* Calculation of the minimum eigenvalue failed. Return the */
837 /* Steihaug-Toint direction. */
838 /*************************************************************************/
839
840 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");CHKERRQ(ierr);
841 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
842 PetscFunctionReturn(0);
843 }
844
845 cg->e_min = e_valu[0];
846
847 /***************************************************************************/
848 /* Compute the initial value of lambda to make (T + lamba I) positive */
849 /* definite. */
850 /***************************************************************************/
851
852 pert = cg->init_pert;
853 if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
854
855 while (1) {
856 for (i = 0; i < t_size; ++i) {
857 t_diag[i] = cg->diag[i] + cg->lambda;
858 t_offd[i] = cg->offd[i];
859 }
860
861 PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
862
863 if (0 == info) break;
864
865 pert += pert;
866 cg->lambda = cg->lambda * (1.0 + pert) + pert;
867 }
868
869 /***************************************************************************/
870 /* Compute the initial step and its norm. */
871 /***************************************************************************/
872
873 nrhs = 1;
874 nldb = t_size;
875
876 t_soln[0] = -cg->norm_r[0];
877 for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
878
879 PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
880
881 if (0 != info) {
882 /*************************************************************************/
883 /* Calculation of the initial step failed; return the Steihaug-Toint */
884 /* direction. */
885 /*************************************************************************/
886
887 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");CHKERRQ(ierr);
888 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
889 PetscFunctionReturn(0);
890 }
891
892 norm_t = 0.;
893 for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
894 norm_t = PetscSqrtReal(norm_t);
895
896 /***************************************************************************/
897 /* Determine the case we are in. */
898 /***************************************************************************/
899
900 if (norm_t <= cg->radius) {
901 /*************************************************************************/
902 /* The step is within the trust region; check if we are in the hard case */
903 /* and need to move to the boundary by following a direction of negative */
904 /* curvature. */
905 /*************************************************************************/
906
907 if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
908 /***********************************************************************/
909 /* This is the hard case; compute the eigenvector associated with the */
910 /* minimum eigenvalue and move along this direction to the boundary. */
911 /***********************************************************************/
912
913 PetscStackCallBLAS("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));
914
915 if (0 != info) {
916 /*********************************************************************/
917 /* Calculation of the minimum eigenvalue failed. Return the */
918 /* Steihaug-Toint direction. */
919 /*********************************************************************/
920
921 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");CHKERRQ(ierr);
922 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
923 PetscFunctionReturn(0);
924 }
925
926 coef1 = 0.0;
927 coef2 = 0.0;
928 coef3 = -cg->radius * cg->radius;
929 for (i = 0; i < t_size; ++i) {
930 coef1 += e_vect[i] * e_vect[i];
931 coef2 += e_vect[i] * t_soln[i];
932 coef3 += t_soln[i] * t_soln[i];
933 }
934
935 coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
936 root1 = (-coef2 + coef3) / coef1;
937 root2 = (-coef2 - coef3) / coef1;
938
939 /***********************************************************************/
940 /* Compute objective value for (t_soln + root1 * e_vect) */
941 /***********************************************************************/
942
943 for (i = 0; i < t_size; ++i) {
944 e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
945 }
946
947 obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
948 cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
949 for (i = 1; i < t_size - 1; ++i) {
950 obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
951 cg->diag[i]*e_rwrk[i]+
952 cg->offd[i+1]*e_rwrk[i+1]);
953 }
954 obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
955 cg->diag[i]*e_rwrk[i]);
956
957 /***********************************************************************/
958 /* Compute objective value for (t_soln + root2 * e_vect) */
959 /***********************************************************************/
960
961 for (i = 0; i < t_size; ++i) {
962 e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
963 }
964
965 obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
966 cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
967 for (i = 1; i < t_size - 1; ++i) {
968 obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
969 cg->diag[i]*e_rwrk[i]+
970 cg->offd[i+1]*e_rwrk[i+1]);
971 }
972 obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
973 cg->diag[i]*e_rwrk[i]);
974
975 /***********************************************************************/
976 /* Choose the point with the best objective function value. */
977 /***********************************************************************/
978
979 if (obj1 <= obj2) {
980 for (i = 0; i < t_size; ++i) {
981 t_soln[i] += root1 * e_vect[i];
982 }
983 }
984 else {
985 for (i = 0; i < t_size; ++i) {
986 t_soln[i] += root2 * e_vect[i];
987 }
988 }
989 } else {
990 /***********************************************************************/
991 /* The matrix is positive definite or there was no room to move; the */
992 /* solution is already contained in t_soln. */
993 /***********************************************************************/
994 }
995 } else {
996 /*************************************************************************/
997 /* The step is outside the trust-region. Compute the correct value for */
998 /* lambda by performing Newton's method. */
999 /*************************************************************************/
1000
1001 for (i = 0; i < max_newton_its; ++i) {
1002 /***********************************************************************/
1003 /* Check for convergence. */
1004 /***********************************************************************/
1005
1006 if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1007
1008 /***********************************************************************/
1009 /* Compute the update. */
1010 /***********************************************************************/
1011
1012 ierr = PetscArraycpy(e_rwrk, t_soln, t_size);CHKERRQ(ierr);
1013
1014 PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1015
1016 if (0 != info) {
1017 /*********************************************************************/
1018 /* Calculation of the step failed; return the Steihaug-Toint */
1019 /* direction. */
1020 /*********************************************************************/
1021
1022 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");CHKERRQ(ierr);
1023 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1024 PetscFunctionReturn(0);
1025 }
1026
1027 /***********************************************************************/
1028 /* Modify lambda. */
1029 /***********************************************************************/
1030
1031 norm_w = 0.;
1032 for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1033
1034 cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1035
1036 /***********************************************************************/
1037 /* Factor T + lambda I */
1038 /***********************************************************************/
1039
1040 for (j = 0; j < t_size; ++j) {
1041 t_diag[j] = cg->diag[j] + cg->lambda;
1042 t_offd[j] = cg->offd[j];
1043 }
1044
1045 PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1046
1047 if (0 != info) {
1048 /*********************************************************************/
1049 /* Calculation of factorization failed; return the Steihaug-Toint */
1050 /* direction. */
1051 /*********************************************************************/
1052
1053 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");CHKERRQ(ierr);
1054 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1055 PetscFunctionReturn(0);
1056 }
1057
1058 /***********************************************************************/
1059 /* Compute the new step and its norm. */
1060 /***********************************************************************/
1061
1062 t_soln[0] = -cg->norm_r[0];
1063 for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1064
1065 PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1066
1067 if (0 != info) {
1068 /*********************************************************************/
1069 /* Calculation of the step failed; return the Steihaug-Toint */
1070 /* direction. */
1071 /*********************************************************************/
1072
1073 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");CHKERRQ(ierr);
1074 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1075 PetscFunctionReturn(0);
1076 }
1077
1078 norm_t = 0.;
1079 for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1080 norm_t = PetscSqrtReal(norm_t);
1081 }
1082
1083 /*************************************************************************/
1084 /* Check for convergence. */
1085 /*************************************************************************/
1086
1087 if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1088 /***********************************************************************/
1089 /* Newton method failed to converge in iteration limit. */
1090 /***********************************************************************/
1091
1092 ierr = PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");CHKERRQ(ierr);
1093 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1094 PetscFunctionReturn(0);
1095 }
1096 }
1097
1098 /***************************************************************************/
1099 /* Recover the norm of the direction and objective function value. */
1100 /***************************************************************************/
1101
1102 cg->norm_d = norm_t;
1103
1104 cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1105 for (i = 1; i < t_size - 1; ++i) {
1106 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]);
1107 }
1108 cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1109
1110 /***************************************************************************/
1111 /* Recover the direction. */
1112 /***************************************************************************/
1113
1114 sigma = -1;
1115
1116 /***************************************************************************/
1117 /* Start conjugate gradient method from the beginning */
1118 /***************************************************************************/
1119
1120 ierr = VecCopy(ksp->vec_rhs, r);CHKERRQ(ierr); /* r = -grad */
1121 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
1122
1123 /***************************************************************************/
1124 /* Accumulate Q * s */
1125 /***************************************************************************/
1126
1127 ierr = VecCopy(z, d);CHKERRQ(ierr);
1128 ierr = VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);CHKERRQ(ierr);
1129
1130 /***************************************************************************/
1131 /* Compute the first direction. */
1132 /***************************************************************************/
1133
1134 ierr = VecCopy(z, p);CHKERRQ(ierr); /* p = z */
1135 ierr = KSP_MatMult(ksp, Qmat, p, z);CHKERRQ(ierr); /* z = Q * p */
1136 ++ksp->its;
1137
1138 for (i = 0; i < l_size - 1; ++i) {
1139 /*************************************************************************/
1140 /* Update the residual and direction. */
1141 /*************************************************************************/
1142
1143 alpha = cg->alpha[i];
1144 if (alpha >= 0.0) sigma = -sigma;
1145
1146 ierr = VecAXPY(r, -alpha, z);CHKERRQ(ierr); /* r = r - alpha Q p */
1147 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
1148
1149 /*************************************************************************/
1150 /* Accumulate Q * s */
1151 /*************************************************************************/
1152
1153 ierr = VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);CHKERRQ(ierr);
1154
1155 /*************************************************************************/
1156 /* Update p. */
1157 /*************************************************************************/
1158
1159 beta = cg->beta[i];
1160 ierr = VecAYPX(p, beta, z);CHKERRQ(ierr); /* p = z + beta p */
1161 ierr = KSP_MatMult(ksp, Qmat, p, z);CHKERRQ(ierr); /* z = Q * p */
1162 ++ksp->its;
1163 }
1164
1165 /***************************************************************************/
1166 /* Update the residual and direction. */
1167 /***************************************************************************/
1168
1169 alpha = cg->alpha[i];
1170 if (alpha >= 0.0) sigma = -sigma;
1171
1172 ierr = VecAXPY(r, -alpha, z);CHKERRQ(ierr); /* r = r - alpha Q p */
1173 ierr = KSP_PCApply(ksp, r, z);CHKERRQ(ierr); /* z = inv(M) r */
1174
1175 /***************************************************************************/
1176 /* Accumulate Q * s */
1177 /***************************************************************************/
1178
1179 ierr = VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);CHKERRQ(ierr);
1180
1181 /***************************************************************************/
1182 /* Set the termination reason. */
1183 /***************************************************************************/
1184
1185 ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1186 PetscFunctionReturn(0);
1187 #endif
1188 }
1189
KSPCGSetUp_GLTR(KSP ksp)1190 static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1191 {
1192 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1193 PetscInt max_its;
1194 PetscErrorCode ierr;
1195
1196 PetscFunctionBegin;
1197 /***************************************************************************/
1198 /* Determine the total maximum number of iterations. */
1199 /***************************************************************************/
1200
1201 max_its = ksp->max_it + cg->max_lanczos_its + 1;
1202
1203 /***************************************************************************/
1204 /* Set work vectors needed by conjugate gradient method and allocate */
1205 /* workspace for Lanczos matrix. */
1206 /***************************************************************************/
1207
1208 ierr = KSPSetWorkVecs(ksp, 3);CHKERRQ(ierr);
1209 if (cg->diag) {
1210 ierr = PetscArrayzero(cg->diag, max_its);CHKERRQ(ierr);
1211 ierr = PetscArrayzero(cg->offd, max_its);CHKERRQ(ierr);
1212 ierr = PetscArrayzero(cg->alpha, max_its);CHKERRQ(ierr);
1213 ierr = PetscArrayzero(cg->beta, max_its);CHKERRQ(ierr);
1214 ierr = PetscArrayzero(cg->norm_r, max_its);CHKERRQ(ierr);
1215 } else {
1216 ierr = PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);CHKERRQ(ierr);
1217 ierr = PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));CHKERRQ(ierr);
1218 }
1219 PetscFunctionReturn(0);
1220 }
1221
KSPCGDestroy_GLTR(KSP ksp)1222 static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1223 {
1224 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1225 PetscErrorCode ierr;
1226
1227 PetscFunctionBegin;
1228 /***************************************************************************/
1229 /* Free memory allocated for the data. */
1230 /***************************************************************************/
1231
1232 ierr = PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);CHKERRQ(ierr);
1233 if (cg->alloced) {
1234 ierr = PetscFree2(cg->rwork,cg->iwork);CHKERRQ(ierr);
1235 }
1236
1237 /***************************************************************************/
1238 /* Clear composed functions */
1239 /***************************************************************************/
1240
1241 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);CHKERRQ(ierr);
1242 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);CHKERRQ(ierr);
1243 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);CHKERRQ(ierr);
1244 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);CHKERRQ(ierr);
1245 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);CHKERRQ(ierr);
1246
1247 /***************************************************************************/
1248 /* Destroy KSP object. */
1249 /***************************************************************************/
1250 ierr = KSPDestroyDefault(ksp);CHKERRQ(ierr);
1251 PetscFunctionReturn(0);
1252 }
1253
KSPCGSetRadius_GLTR(KSP ksp,PetscReal radius)1254 static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1255 {
1256 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1257
1258 PetscFunctionBegin;
1259 cg->radius = radius;
1260 PetscFunctionReturn(0);
1261 }
1262
KSPCGGetNormD_GLTR(KSP ksp,PetscReal * norm_d)1263 static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1264 {
1265 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1266
1267 PetscFunctionBegin;
1268 *norm_d = cg->norm_d;
1269 PetscFunctionReturn(0);
1270 }
1271
KSPCGGetObjFcn_GLTR(KSP ksp,PetscReal * o_fcn)1272 static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1273 {
1274 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1275
1276 PetscFunctionBegin;
1277 *o_fcn = cg->o_fcn;
1278 PetscFunctionReturn(0);
1279 }
1280
KSPGLTRGetMinEig_GLTR(KSP ksp,PetscReal * e_min)1281 static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1282 {
1283 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1284
1285 PetscFunctionBegin;
1286 *e_min = cg->e_min;
1287 PetscFunctionReturn(0);
1288 }
1289
KSPGLTRGetLambda_GLTR(KSP ksp,PetscReal * lambda)1290 static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1291 {
1292 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1293
1294 PetscFunctionBegin;
1295 *lambda = cg->lambda;
1296 PetscFunctionReturn(0);
1297 }
1298
KSPCGSetFromOptions_GLTR(PetscOptionItems * PetscOptionsObject,KSP ksp)1299 static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1300 {
1301 PetscErrorCode ierr;
1302 KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1303
1304 PetscFunctionBegin;
1305 ierr = PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");CHKERRQ(ierr);
1306
1307 ierr = PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);CHKERRQ(ierr);
1308
1309 ierr = PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);CHKERRQ(ierr);
1310
1311 ierr = PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);CHKERRQ(ierr);
1312 ierr = PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);CHKERRQ(ierr);
1313 ierr = PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);CHKERRQ(ierr);
1314
1315 ierr = PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);CHKERRQ(ierr);
1316 ierr = PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);CHKERRQ(ierr);
1317
1318 ierr = PetscOptionsTail();CHKERRQ(ierr);
1319 PetscFunctionReturn(0);
1320 }
1321
1322 /*MC
1323 KSPGLTR - Code to run conjugate gradient method subject to a constraint
1324 on the solution norm. This is used in Trust Region methods for
1325 nonlinear equations, SNESNEWTONTR
1326
1327 Options Database Keys:
1328 . -ksp_cg_radius <r> - Trust Region Radius
1329
1330 Notes:
1331 This is rarely used directly
1332
1333 Use preconditioned conjugate gradient to compute
1334 an approximate minimizer of the quadratic function
1335
1336 q(s) = g^T * s + .5 * s^T * H * s
1337
1338 subject to the trust region constraint
1339
1340 || s || <= delta,
1341
1342 where
1343
1344 delta is the trust region radius,
1345 g is the gradient vector,
1346 H is the Hessian approximation,
1347 M is the positive definite preconditioner matrix.
1348
1349 KSPConvergedReason may be
1350 $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1351 $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1352 $ other KSP converged/diverged reasons
1353
1354 Notes:
1355 The preconditioner supplied should be symmetric and positive definite.
1356
1357 Reference:
1358 Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1359 SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525
1360
1361 Level: developer
1362
1363 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1364 M*/
1365
KSPCreate_GLTR(KSP ksp)1366 PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1367 {
1368 PetscErrorCode ierr;
1369 KSPCG_GLTR *cg;
1370
1371 PetscFunctionBegin;
1372 ierr = PetscNewLog(ksp,&cg);CHKERRQ(ierr);
1373 cg->radius = 0.0;
1374 cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1375
1376 cg->init_pert = 1.0e-8;
1377 cg->eigen_tol = 1.0e-10;
1378 cg->newton_tol = 1.0e-6;
1379
1380 cg->alloced = 0;
1381 cg->init_alloc = 1024;
1382
1383 cg->max_lanczos_its = 20;
1384 cg->max_newton_its = 10;
1385
1386 ksp->data = (void*) cg;
1387 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
1388 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
1389 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);CHKERRQ(ierr);
1390 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
1391
1392 /***************************************************************************/
1393 /* Sets the functions that are associated with this data structure */
1394 /* (in C++ this is the same as defining virtual functions). */
1395 /***************************************************************************/
1396
1397 ksp->ops->setup = KSPCGSetUp_GLTR;
1398 ksp->ops->solve = KSPCGSolve_GLTR;
1399 ksp->ops->destroy = KSPCGDestroy_GLTR;
1400 ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1401 ksp->ops->buildsolution = KSPBuildSolutionDefault;
1402 ksp->ops->buildresidual = KSPBuildResidualDefault;
1403 ksp->ops->view = NULL;
1404
1405 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);CHKERRQ(ierr);
1406 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);CHKERRQ(ierr);
1407 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);CHKERRQ(ierr);
1408 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);CHKERRQ(ierr);
1409 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);CHKERRQ(ierr);
1410 PetscFunctionReturn(0);
1411 }
1412