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