Actual source code: ntl.c
1: #include "src/matrix/lmvmmat.h"
2: #include "ntl.h"
4: #include "petscksp.h"
5: #include "petscpc.h"
6: #include "private/kspimpl.h"
7: #include "private/pcimpl.h"
9: #define NTL_KSP_NASH 0
10: #define NTL_KSP_STCG 1
11: #define NTL_KSP_GLTR 2
12: #define NTL_KSP_TYPES 3
14: #define NTL_PC_NONE 0
15: #define NTL_PC_AHESS 1
16: #define NTL_PC_BFGS 2
17: #define NTL_PC_PETSC 3
18: #define NTL_PC_TYPES 4
20: #define BFGS_SCALE_AHESS 0
21: #define BFGS_SCALE_BFGS 1
22: #define BFGS_SCALE_TYPES 2
24: #define NTL_INIT_CONSTANT 0
25: #define NTL_INIT_DIRECTION 1
26: #define NTL_INIT_INTERPOLATION 2
27: #define NTL_INIT_TYPES 3
29: #define NTL_UPDATE_REDUCTION 0
30: #define NTL_UPDATE_INTERPOLATION 1
31: #define NTL_UPDATE_TYPES 2
33: static const char *NTL_KSP[64] = {
34: "nash", "stcg", "gltr"
35: };
37: static const char *NTL_PC[64] = {
38: "none", "ahess", "bfgs", "petsc"
39: };
41: static const char *BFGS_SCALE[64] = {
42: "ahess", "bfgs"
43: };
45: static const char *NTL_INIT[64] = {
46: "constant", "direction", "interpolation"
47: };
49: static const char *NTL_UPDATE[64] = {
50: "reduction", "interpolation"
51: };
53: /* Routine for BFGS preconditioner */
57: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
58: {
60: Mat M;
66: PCShellGetContext(pc,(void**)&M);
67: MatLMVMSolve(M, b, x);
68: return(0);
69: }
71: /* Implements Newton's Method with a trust-region, line-search approach for
72: solving unconstrained minimization problems. A More'-Thuente line search
73: is used to guarantee that the bfgs preconditioner remains positive
74: definite. */
76: #define NTL_NEWTON 0
77: #define NTL_BFGS 1
78: #define NTL_SCALED_GRADIENT 2
79: #define NTL_GRADIENT 3
83: static PetscErrorCode TaoSolve_NTL(TaoSolver tao)
84: {
85: TAO_NTL *tl = (TAO_NTL *)tao->data;
87: PC pc;
88: KSPConvergedReason ksp_reason;
89: TaoSolverTerminationReason reason;
90: TaoLineSearchTerminationReason ls_reason;
92: PetscReal fmin, ftrial, prered, actred, kappa, sigma;
93: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
94: PetscReal f, fold, gdx, gnorm;
95: PetscReal step = 1.0;
97: PetscReal delta;
98: PetscReal norm_d = 0.0;
99: MatStructure matflag;
101: PetscInt stepType;
102: PetscInt iter = 0,its;
104: PetscInt bfgsUpdates = 0;
105: PetscInt needH;
107: PetscInt i_max = 5;
108: PetscInt j_max = 1;
109: PetscInt i, j, n, N;
111: PetscInt tr_reject;
115: if (tao->XL || tao->XU || tao->ops->computebounds) {
116: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n");
117: }
119: /* Initialize trust-region radius */
120: tao->trust = tao->trust0;
122: /* Modify the radius if it is too large or small */
123: tao->trust = PetscMax(tao->trust, tl->min_radius);
124: tao->trust = PetscMin(tao->trust, tl->max_radius);
126: if (NTL_PC_BFGS == tl->pc_type && !tl->M) {
127: VecGetLocalSize(tao->solution,&n);
128: VecGetSize(tao->solution,&N);
129: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tl->M);
130: MatLMVMAllocateVectors(tl->M,tao->solution);
131: }
133: /* Check convergence criteria */
134: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
135: VecNorm(tao->gradient, NORM_2, &gnorm);
136: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
137: SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
138: }
139: needH = 1;
141: TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);
142: if (reason != TAO_CONTINUE_ITERATING) {
143: return(0);
144: }
146: /* Create vectors for the limited memory preconditioner */
147: if ((NTL_PC_BFGS == tl->pc_type) &&
148: (BFGS_SCALE_BFGS != tl->bfgs_scale_type)) {
149: if (!tl->Diag) {
150: VecDuplicate(tao->solution, &tl->Diag);
151: }
152: }
154: /* Modify the linear solver to a conjugate gradient method */
155: switch(tl->ksp_type) {
156: case NTL_KSP_NASH:
157: KSPSetType(tao->ksp, KSPNASH);
158: if (tao->ksp->ops->setfromoptions) {
159: (*tao->ksp->ops->setfromoptions)(tao->ksp);
160: }
161: break;
163: case NTL_KSP_STCG:
164: KSPSetType(tao->ksp, KSPSTCG);
165: if (tao->ksp->ops->setfromoptions) {
166: (*tao->ksp->ops->setfromoptions)(tao->ksp);
167: }
168: break;
170: default:
171: KSPSetType(tao->ksp, KSPGLTR);
172: if (tao->ksp->ops->setfromoptions) {
173: (*tao->ksp->ops->setfromoptions)(tao->ksp);
174: }
175: break;
176: }
178: /* Modify the preconditioner to use the bfgs approximation */
179: KSPGetPC(tao->ksp, &pc);
180: switch(tl->pc_type) {
181: case NTL_PC_NONE:
182: PCSetType(pc, PCNONE);
183: if (pc->ops->setfromoptions) {
184: (*pc->ops->setfromoptions)(pc);
185: }
186: break;
188: case NTL_PC_AHESS:
189: PCSetType(pc, PCJACOBI);
190: if (pc->ops->setfromoptions) {
191: (*pc->ops->setfromoptions)(pc);
192: }
193: PCJacobiSetUseAbs(pc);
194: break;
196: case NTL_PC_BFGS:
197: PCSetType(pc, PCSHELL);
198: if (pc->ops->setfromoptions) {
199: (*pc->ops->setfromoptions)(pc);
200: }
201: PCShellSetName(pc, "bfgs");
202: PCShellSetContext(pc, tl->M);
203: PCShellSetApply(pc, MatLMVMSolveShell);
204: break;
206: default:
207: /* Use the pc method set by pc_type */
208: break;
209: }
211: /* Initialize trust-region radius. The initialization is only performed
212: when we are using Steihaug-Toint or the Generalized Lanczos method. */
213: switch(tl->init_type) {
214: case NTL_INIT_CONSTANT:
215: /* Use the initial radius specified */
216: break;
218: case NTL_INIT_INTERPOLATION:
219: /* Use the initial radius specified */
220: max_radius = 0.0;
221:
222: for (j = 0; j < j_max; ++j) {
223: fmin = f;
224: sigma = 0.0;
225:
226: if (needH) {
227: TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);
228: needH = 0;
229: }
230:
231: for (i = 0; i < i_max; ++i) {
232: VecCopy(tao->solution, tl->W);
233: VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient);
235: TaoComputeObjective(tao, tl->W, &ftrial);
236: if (PetscIsInfOrNanReal(ftrial)) {
237: tau = tl->gamma1_i;
238: }
239: else {
240: if (ftrial < fmin) {
241: fmin = ftrial;
242: sigma = -tao->trust / gnorm;
243: }
245: MatMult(tao->hessian, tao->gradient, tao->stepdirection);
246: VecDot(tao->gradient, tao->stepdirection, &prered);
248: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
249: actred = f - ftrial;
250: if ((PetscAbsScalar(actred) <= tl->epsilon) &&
251: (PetscAbsScalar(prered) <= tl->epsilon)) {
252: kappa = 1.0;
253: }
254: else {
255: kappa = actred / prered;
256: }
258: tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
259: tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
260: tau_min = PetscMin(tau_1, tau_2);
261: tau_max = PetscMax(tau_1, tau_2);
263: if (PetscAbsScalar(kappa - 1.0) <= tl->mu1_i) {
264: /* Great agreement */
265: max_radius = PetscMax(max_radius, tao->trust);
267: if (tau_max < 1.0) {
268: tau = tl->gamma3_i;
269: }
270: else if (tau_max > tl->gamma4_i) {
271: tau = tl->gamma4_i;
272: }
273: else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
274: tau = tau_1;
275: }
276: else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
277: tau = tau_2;
278: }
279: else {
280: tau = tau_max;
281: }
282: }
283: else if (PetscAbsScalar(kappa - 1.0) <= tl->mu2_i) {
284: /* Good agreement */
285: max_radius = PetscMax(max_radius, tao->trust);
287: if (tau_max < tl->gamma2_i) {
288: tau = tl->gamma2_i;
289: }
290: else if (tau_max > tl->gamma3_i) {
291: tau = tl->gamma3_i;
292: }
293: else {
294: tau = tau_max;
295: }
296: }
297: else {
298: /* Not good agreement */
299: if (tau_min > 1.0) {
300: tau = tl->gamma2_i;
301: }
302: else if (tau_max < tl->gamma1_i) {
303: tau = tl->gamma1_i;
304: }
305: else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
306: tau = tl->gamma1_i;
307: }
308: else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) &&
309: ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
310: tau = tau_1;
311: }
312: else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) &&
313: ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
314: tau = tau_2;
315: }
316: else {
317: tau = tau_max;
318: }
319: }
320: }
321: tao->trust = tau * tao->trust;
322: }
323:
324: if (fmin < f) {
325: f = fmin;
326: VecAXPY(tao->solution, sigma, tao->gradient);
327: TaoComputeGradient(tao, tao->solution, tao->gradient);
329: VecNorm(tao->gradient, NORM_2, &gnorm);
330: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
331: SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
332: }
333: needH = 1;
334:
335: TaoMonitor(tao, iter, f, gnorm, 0.0, 1.0, &reason);
336: if (reason != TAO_CONTINUE_ITERATING) {
337: return(0);
338: }
339: }
340: }
341: tao->trust = PetscMax(tao->trust, max_radius);
343: /* Modify the radius if it is too large or small */
344: tao->trust = PetscMax(tao->trust, tl->min_radius);
345: tao->trust = PetscMin(tao->trust, tl->max_radius);
346: break;
348: default:
349: /* Norm of the first direction will initialize radius */
350: tao->trust = 0.0;
351: break;
352: }
354: /* Set initial scaling for the BFGS preconditioner
355: This step is done after computing the initial trust-region radius
356: since the function value may have decreased */
357: if (NTL_PC_BFGS == tl->pc_type) {
358: if (f != 0.0) {
359: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
360: }
361: else {
362: delta = 2.0 / (gnorm*gnorm);
363: }
364: MatLMVMSetDelta(tl->M, delta);
365: }
367: /* Set counter for gradient/reset steps */
368: tl->ntrust = 0;
369: tl->newt = 0;
370: tl->bfgs = 0;
371: tl->sgrad = 0;
372: tl->grad = 0;
374: /* Have not converged; continue with Newton method */
375: while (reason == TAO_CONTINUE_ITERATING) {
376: ++iter;
378: /* Compute the Hessian */
379: if (needH) {
380: TaoComputeHessian(tao, tao->solution, &tao->hessian, &tao->hessian_pre, &matflag);
381: needH = 0;
382: }
384: if (NTL_PC_BFGS == tl->pc_type) {
385: if (BFGS_SCALE_AHESS == tl->bfgs_scale_type) {
386: /* Obtain diagonal for the bfgs preconditioner */
387: MatGetDiagonal(tao->hessian, tl->Diag);
388: VecAbs(tl->Diag);
389: VecReciprocal(tl->Diag);
390: MatLMVMSetScale(tl->M, tl->Diag);
391: }
393: /* Update the limited memory preconditioner */
394: MatLMVMUpdate(tl->M,tao->solution, tao->gradient);
395: ++bfgsUpdates;
396: }
397: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre, matflag);
398: /* Solve the Newton system of equations */
399: if (NTL_KSP_NASH == tl->ksp_type) {
400: KSPNASHSetRadius(tao->ksp,tl->max_radius);
401: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
402: KSPGetIterationNumber(tao->ksp,&its);
403: tao->ksp_its+=its;
404: KSPNASHGetNormD(tao->ksp, &norm_d);
405: } else if (NTL_KSP_STCG == tl->ksp_type) {
406: KSPSTCGSetRadius(tao->ksp,tl->max_radius);
407: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
408: KSPGetIterationNumber(tao->ksp,&its);
409: tao->ksp_its+=its;
410: KSPSTCGGetNormD(tao->ksp, &norm_d);
411: } else { /* NTL_KSP_GLTR */
412: KSPGLTRSetRadius(tao->ksp,tl->max_radius);
413: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
414: KSPGetIterationNumber(tao->ksp,&its);
415: tao->ksp_its+=its;
416: KSPGLTRGetNormD(tao->ksp, &norm_d);
417: }
419: if (0.0 == tao->trust) {
420: /* Radius was uninitialized; use the norm of the direction */
421: if (norm_d > 0.0) {
422: tao->trust = norm_d;
424: /* Modify the radius if it is too large or small */
425: tao->trust = PetscMax(tao->trust, tl->min_radius);
426: tao->trust = PetscMin(tao->trust, tl->max_radius);
427: }
428: else {
429: /* The direction was bad; set radius to default value and re-solve
430: the trust-region subproblem to get a direction */
431: tao->trust = tao->trust0;
433: /* Modify the radius if it is too large or small */
434: tao->trust = PetscMax(tao->trust, tl->min_radius);
435: tao->trust = PetscMin(tao->trust, tl->max_radius);
437: if (NTL_KSP_NASH == tl->ksp_type) {
438: KSPNASHSetRadius(tao->ksp,tl->max_radius);
439: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
440: KSPGetIterationNumber(tao->ksp,&its);
441: tao->ksp_its+=its;
442: KSPNASHGetNormD(tao->ksp, &norm_d);
443: } else if (NTL_KSP_STCG == tl->ksp_type) {
444: KSPSTCGSetRadius(tao->ksp,tl->max_radius);
445: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
446: KSPGetIterationNumber(tao->ksp,&its);
447: tao->ksp_its+=its;
448: KSPSTCGGetNormD(tao->ksp, &norm_d);
449: } else { /* NTL_KSP_GLTR */
450: KSPGLTRSetRadius(tao->ksp,tl->max_radius);
451: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
452: KSPGetIterationNumber(tao->ksp,&its);
453: tao->ksp_its+=its;
454: KSPGLTRGetNormD(tao->ksp, &norm_d);
455: }
458: if (norm_d == 0.0) {
459: SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
460: }
461: }
462: }
464: VecScale(tao->stepdirection, -1.0);
465: KSPGetConvergedReason(tao->ksp, &ksp_reason);
466: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) &&
467: (NTL_PC_BFGS == tl->pc_type) && (bfgsUpdates > 1)) {
468: /* Preconditioner is numerically indefinite; reset the
469: approximate if using BFGS preconditioning. */
471: if (f != 0.0) {
472: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
473: }
474: else {
475: delta = 2.0 / (gnorm*gnorm);
476: }
477: MatLMVMSetDelta(tl->M, delta);
478: MatLMVMReset(tl->M);
479: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
480: bfgsUpdates = 1;
481: }
483: /* Check trust-region reduction conditions */
484: tr_reject = 0;
485: if (NTL_UPDATE_REDUCTION == tl->update_type) {
486: /* Get predicted reduction */
487: if (NTL_KSP_NASH == tl->ksp_type) {
488: KSPNASHGetObjFcn(tao->ksp,&prered);
489: } else if (NTL_KSP_STCG == tl->ksp_type) {
490: KSPSTCGGetObjFcn(tao->ksp,&prered);
491: } else { /* gltr */
492: KSPGLTRGetObjFcn(tao->ksp,&prered);
493: }
495: if (prered >= 0.0) {
496: /* The predicted reduction has the wrong sign. This cannot
497: happen in infinite precision arithmetic. Step should
498: be rejected! */
499: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
500: tr_reject = 1;
501: }
502: else {
503: /* Compute trial step and function value */
504: VecCopy(tao->solution, tl->W);
505: VecAXPY(tl->W, 1.0, tao->stepdirection);
506: TaoComputeObjective(tao, tl->W, &ftrial);
508: if (PetscIsInfOrNanReal(ftrial)) {
509: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
510: tr_reject = 1;
511: }
512: else {
513: /* Compute and actual reduction */
514: actred = f - ftrial;
515: prered = -prered;
516: if ((PetscAbsScalar(actred) <= tl->epsilon) &&
517: (PetscAbsScalar(prered) <= tl->epsilon)) {
518: kappa = 1.0;
519: }
520: else {
521: kappa = actred / prered;
522: }
524: /* Accept of reject the step and update radius */
525: if (kappa < tl->eta1) {
526: /* Reject the step */
527: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
528: tr_reject = 1;
529: }
530: else {
531: /* Accept the step */
532: if (kappa < tl->eta2) {
533: /* Marginal bad step */
534: tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
535: }
536: else if (kappa < tl->eta3) {
537: /* Reasonable step */
538: tao->trust = tl->alpha3 * tao->trust;
539: }
540: else if (kappa < tl->eta4) {
541: /* Good step */
542: tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
543: }
544: else {
545: /* Very good step */
546: tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
547: }
548: }
549: }
550: }
551: }
552: else {
553: /* Get predicted reduction */
554: if (NTL_KSP_NASH == tl->ksp_type) {
555: KSPNASHGetObjFcn(tao->ksp,&prered);
556: } else if (NTL_KSP_STCG == tl->ksp_type) {
557: KSPSTCGGetObjFcn(tao->ksp,&prered);
558: } else { /* gltr */
559: KSPGLTRGetObjFcn(tao->ksp,&prered);
560: }
562: if (prered >= 0.0) {
563: /* The predicted reduction has the wrong sign. This cannot
564: happen in infinite precision arithmetic. Step should
565: be rejected! */
566: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
567: tr_reject = 1;
568: }
569: else {
570: VecCopy(tao->solution, tl->W);
571: VecAXPY(tl->W, 1.0, tao->stepdirection);
572: TaoComputeObjective(tao, tl->W, &ftrial);
573: if (PetscIsInfOrNanReal(ftrial)) {
574: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
575: tr_reject = 1;
576: }
577: else {
578: VecDot(tao->gradient, tao->stepdirection, &gdx);
580: actred = f - ftrial;
581: prered = -prered;
582: if ((PetscAbsScalar(actred) <= tl->epsilon) &&
583: (PetscAbsScalar(prered) <= tl->epsilon)) {
584: kappa = 1.0;
585: }
586: else {
587: kappa = actred / prered;
588: }
590: tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
591: tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
592: tau_min = PetscMin(tau_1, tau_2);
593: tau_max = PetscMax(tau_1, tau_2);
595: if (kappa >= 1.0 - tl->mu1) {
596: /* Great agreement; accept step and update radius */
597: if (tau_max < 1.0) {
598: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
599: }
600: else if (tau_max > tl->gamma4) {
601: tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
602: }
603: else {
604: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
605: }
606: }
607: else if (kappa >= 1.0 - tl->mu2) {
608: /* Good agreement */
610: if (tau_max < tl->gamma2) {
611: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
612: }
613: else if (tau_max > tl->gamma3) {
614: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
615: } else if (tau_max < 1.0) {
616: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
617: }
618: else {
619: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
620: }
621: }
622: else {
623: /* Not good agreement */
624: if (tau_min > 1.0) {
625: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
626: }
627: else if (tau_max < tl->gamma1) {
628: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
629: }
630: else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
631: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
632: }
633: else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) &&
634: ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
635: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
636: }
637: else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) &&
638: ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
639: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
640: }
641: else {
642: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
643: }
644: tr_reject = 1;
645: }
646: }
647: }
648: }
650: if (tr_reject) {
651: /* The trust-region constraints rejected the step. Apply a linesearch.
652: Check for descent direction. */
653: VecDot(tao->stepdirection, tao->gradient, &gdx);
654: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
655: /* Newton step is not descent or direction produced Inf or NaN */
656:
657: if (NTL_PC_BFGS != tl->pc_type) {
658: /* We don't have the bfgs matrix around and updated
659: Must use gradient direction in this case */
660: VecCopy(tao->gradient, tao->stepdirection);
661: VecScale(tao->stepdirection, -1.0);
662: ++tl->grad;
663: stepType = NTL_GRADIENT;
664: }
665: else {
666: /* Attempt to use the BFGS direction */
667: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
668: VecScale(tao->stepdirection, -1.0);
669:
670: /* Check for success (descent direction) */
671: VecDot(tao->stepdirection, tao->gradient, &gdx);
672: if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
673: /* BFGS direction is not descent or direction produced not a number
674: We can assert bfgsUpdates > 1 in this case because
675: the first solve produces the scaled gradient direction,
676: which is guaranteed to be descent */
677:
678: /* Use steepest descent direction (scaled) */
679: if (f != 0.0) {
680: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
681: }
682: else {
683: delta = 2.0 / (gnorm*gnorm);
684: }
685: MatLMVMSetDelta(tl->M, delta);
686: MatLMVMReset(tl->M);
687: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
688: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
689: VecScale(tao->stepdirection, -1.0);
690:
691: bfgsUpdates = 1;
692: ++tl->sgrad;
693: stepType = NTL_SCALED_GRADIENT;
694: }
695: else {
696: if (1 == bfgsUpdates) {
697: /* The first BFGS direction is always the scaled gradient */
698: ++tl->sgrad;
699: stepType = NTL_SCALED_GRADIENT;
700: }
701: else {
702: ++tl->bfgs;
703: stepType = NTL_BFGS;
704: }
705: }
706: }
707: }
708: else {
709: /* Computed Newton step is descent */
710: ++tl->newt;
711: stepType = NTL_NEWTON;
712: }
713:
714: /* Perform the linesearch */
715: fold = f;
716: VecCopy(tao->solution, tl->Xold);
717: VecCopy(tao->gradient, tl->Gold);
719: step = 1.0;
720: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
721: TaoAddLineSearchCounts(tao);
723:
724: while (ls_reason != TAOLINESEARCH_SUCCESS &&
725: ls_reason != TAOLINESEARCH_SUCCESS_USER &&
726: stepType != NTL_GRADIENT) { /* Linesearch failed */
727: /* Linesearch failed */
728: f = fold;
729: VecCopy(tl->Xold, tao->solution);
730: VecCopy(tl->Gold, tao->gradient);
731:
732: switch(stepType) {
733: case NTL_NEWTON:
734: /* Failed to obtain acceptable iterate with Newton step */
736: if (NTL_PC_BFGS != tl->pc_type) {
737: /* We don't have the bfgs matrix around and being updated
738: Must use gradient direction in this case */
739: VecCopy(tao->gradient, tao->stepdirection);
740: ++tl->grad;
741: stepType = NTL_GRADIENT;
742: }
743: else {
744: /* Attempt to use the BFGS direction */
745: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
747:
748: /* Check for success (descent direction) */
749: VecDot(tao->stepdirection, tao->gradient, &gdx);
750: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
751: /* BFGS direction is not descent or direction produced
752: not a number. We can assert bfgsUpdates > 1 in this case
753: Use steepest descent direction (scaled) */
754:
755: if (f != 0.0) {
756: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
757: }
758: else {
759: delta = 2.0 / (gnorm*gnorm);
760: }
761: MatLMVMSetDelta(tl->M, delta);
762: MatLMVMReset(tl->M);
763: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
764: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
765:
766: bfgsUpdates = 1;
767: ++tl->sgrad;
768: stepType = NTL_SCALED_GRADIENT;
769: }
770: else {
771: if (1 == bfgsUpdates) {
772: /* The first BFGS direction is always the scaled gradient */
773: ++tl->sgrad;
774: stepType = NTL_SCALED_GRADIENT;
775: }
776: else {
777: ++tl->bfgs;
778: stepType = NTL_BFGS;
779: }
780: }
781: }
782: break;
784: case NTL_BFGS:
785: /* Can only enter if pc_type == NTL_PC_BFGS
786: Failed to obtain acceptable iterate with BFGS step
787: Attempt to use the scaled gradient direction */
788:
789: if (f != 0.0) {
790: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
791: }
792: else {
793: delta = 2.0 / (gnorm*gnorm);
794: }
795: MatLMVMSetDelta(tl->M, delta);
796: MatLMVMReset(tl->M);
797: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
798: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
800: bfgsUpdates = 1;
801: ++tl->sgrad;
802: stepType = NTL_SCALED_GRADIENT;
803: break;
804:
805: case NTL_SCALED_GRADIENT:
806: /* Can only enter if pc_type == NTL_PC_BFGS
807: The scaled gradient step did not produce a new iterate;
808: attemp to use the gradient direction.
809: Need to make sure we are not using a different diagonal scaling */
810: MatLMVMSetScale(tl->M, tl->Diag);
811: MatLMVMSetDelta(tl->M, 1.0);
812: MatLMVMReset(tl->M);
813: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
814: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
815:
816: bfgsUpdates = 1;
817: ++tl->grad;
818: stepType = NTL_GRADIENT;
819: break;
820: }
821: VecScale(tao->stepdirection, -1.0);
823: /* This may be incorrect; linesearch has values for stepmax and stepmin
824: that should be reset. */
825: step = 1.0;
826: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
827: TaoAddLineSearchCounts(tao);
828: }
830: if (ls_reason != TAOLINESEARCH_SUCCESS &&
831: ls_reason != TAOLINESEARCH_SUCCESS_USER) {
832: /* Failed to find an improving point */
833: f = fold;
834: VecCopy(tl->Xold, tao->solution);
835: VecCopy(tl->Gold, tao->gradient);
836: tao->trust = 0.0;
837: step = 0.0;
838: reason = TAO_DIVERGED_LS_FAILURE;
839: tao->reason = TAO_DIVERGED_LS_FAILURE;
840: break;
841: }
842: else if (stepType == NTL_NEWTON) {
843: if (step < tl->nu1) {
844: /* Very bad step taken; reduce radius */
845: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
846: }
847: else if (step < tl->nu2) {
848: /* Reasonably bad step taken; reduce radius */
849: tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
850: }
851: else if (step < tl->nu3) {
852: /* Reasonable step was taken; leave radius alone */
853: if (tl->omega3 < 1.0) {
854: tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
855: }
856: else if (tl->omega3 > 1.0) {
857: tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
858: }
859: }
860: else if (step < tl->nu4) {
861: /* Full step taken; increase the radius */
862: tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
863: }
864: else {
865: /* More than full step taken; increase the radius */
866: tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
867: }
868: }
869: else {
870: /* Newton step was not good; reduce the radius */
871: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
872: }
873: }
874: else {
875: /* Trust-region step is accepted */
876: VecCopy(tl->W, tao->solution);
877: f = ftrial;
878: TaoComputeGradient(tao, tao->solution, tao->gradient);
879: ++tl->ntrust;
880: }
882: /* The radius may have been increased; modify if it is too large */
883: tao->trust = PetscMin(tao->trust, tl->max_radius);
885: /* Check for termination */
886: VecNorm(tao->gradient, NORM_2, &gnorm);
887: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
888: SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
889: }
890: needH = 1;
891:
892: TaoMonitor(tao, iter, f, gnorm, 0.0, tao->trust, &reason);
893: }
894: return(0);
895: }
897: /* ---------------------------------------------------------- */
900: static PetscErrorCode TaoSetUp_NTL(TaoSolver tao)
901: {
902: TAO_NTL *tl = (TAO_NTL *)tao->data;
906: if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient); }
907: if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection); }
908: if (!tl->W) { VecDuplicate(tao->solution, &tl->W); }
909: if (!tl->Xold) { VecDuplicate(tao->solution, &tl->Xold); }
910: if (!tl->Gold) { VecDuplicate(tao->solution, &tl->Gold); }
912: tl->Diag = 0;
913: tl->M = 0;
916: return(0);
917: }
919: /*------------------------------------------------------------*/
922: static PetscErrorCode TaoDestroy_NTL(TaoSolver tao)
923: {
924: TAO_NTL *tl = (TAO_NTL *)tao->data;
928: if (tao->setupcalled) {
929: VecDestroy(&tl->W);
930: VecDestroy(&tl->Xold);
931: VecDestroy(&tl->Gold);
932: }
933: if (tl->Diag) {
934: VecDestroy(&tl->Diag);
935: tl->Diag = PETSC_NULL;
936: }
937: if (tl->M) {
938: MatDestroy(&tl->M);
939: tl->M = PETSC_NULL;
940: }
942: PetscFree(tao->data);
943: tao->data = PETSC_NULL;
944:
945: return(0);
946: }
948: /*------------------------------------------------------------*/
951: static PetscErrorCode TaoSetFromOptions_NTL(TaoSolver tao)
952: {
953: TAO_NTL *tl = (TAO_NTL *)tao->data;
957: PetscOptionsHead("Newton line search method for unconstrained optimization");
958: PetscOptionsEList("-tao_ntl_ksp_type", "ksp type", "", NTL_KSP, NTL_KSP_TYPES, NTL_KSP[tl->ksp_type], &tl->ksp_type, 0);
959: PetscOptionsEList("-tao_ntl_pc_type", "pc type", "", NTL_PC, NTL_PC_TYPES, NTL_PC[tl->pc_type], &tl->pc_type, 0);
960: PetscOptionsEList("-tao_ntl_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tl->bfgs_scale_type], &tl->bfgs_scale_type, 0);
961: PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, 0);
962: PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, 0);
963: PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, 0);
964: PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, 0);
965: PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, 0);
966: PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, 0);
967: PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, 0);
968: PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, 0);
969: PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, 0);
970: PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, 0);
971: PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, 0);
972: PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, 0);
973: PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, 0);
974: PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, 0);
975: PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, 0);
976: PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, 0);
977: PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, 0);
978: PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, 0);
979: PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, 0);
980: PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, 0);
981: PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, 0);
982: PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, 0);
983: PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, 0);
984: PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, 0);
985: PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, 0);
986: PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, 0);
987: PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, 0);
988: PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, 0);
989: PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, 0);
990: PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, 0);
991: PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, 0);
992: PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, 0);
993: PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, 0);
994: PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, 0);
995: PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, 0);
996: PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, 0);
997: PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, 0);
998: PetscOptionsTail();
999: TaoLineSearchSetFromOptions(tao->linesearch);
1000: KSPSetFromOptions(tao->ksp);
1001: return(0);
1002: }
1005: /*------------------------------------------------------------*/
1008: static PetscErrorCode TaoView_NTL(TaoSolver tao, PetscViewer viewer)
1009: {
1010: TAO_NTL *tl = (TAO_NTL *)tao->data;
1011: PetscInt nrejects;
1012: PetscBool isascii;
1016: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1017: if (isascii) {
1018: PetscViewerASCIIPushTab(viewer);
1019: if (NTL_PC_BFGS == tl->pc_type && tl->M) {
1020: MatLMVMGetRejects(tl->M, &nrejects);
1021: PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", &nrejects);
1023: }
1025: PetscViewerASCIIPrintf(viewer, "Trust-region steps: %D\n", tl->ntrust);
1026: PetscViewerASCIIPrintf(viewer, "Newton search steps: %D\n", tl->newt);
1027: PetscViewerASCIIPrintf(viewer, "BFGS search steps: %D\n", tl->bfgs);
1028: PetscViewerASCIIPrintf(viewer, "Scaled gradient search steps: %D\n", tl->sgrad);
1029: PetscViewerASCIIPrintf(viewer, "Gradient search steps: %D\n", tl->grad);
1030: PetscViewerASCIIPopTab(viewer);
1031: } else {
1032: SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NTL",((PetscObject)viewer)->type_name);
1033: }
1034: return(0);
1035: }
1037: /* ---------------------------------------------------------- */
1041: PetscErrorCode TaoCreate_NTL(TaoSolver tao)
1042: {
1043: TAO_NTL *tl;
1045: const char *morethuente_type = TAOLINESEARCH_MT;
1047: PetscNewLog(tao, TAO_NTL, &tl);
1048:
1049: tao->ops->setup = TaoSetUp_NTL;
1050: tao->ops->solve = TaoSolve_NTL;
1051: tao->ops->view = TaoView_NTL;
1052: tao->ops->setfromoptions = TaoSetFromOptions_NTL;
1053: tao->ops->destroy = TaoDestroy_NTL;
1054:
1055: tao->max_it = 50;
1056: tao->fatol = 1e-10;
1057: tao->frtol = 1e-10;
1058: tao->data = (void*)tl;
1059:
1060: tao->trust0 = 100.0;
1063: /* Default values for trust-region radius update based on steplength */
1064: tl->nu1 = 0.25;
1065: tl->nu2 = 0.50;
1066: tl->nu3 = 1.00;
1067: tl->nu4 = 1.25;
1069: tl->omega1 = 0.25;
1070: tl->omega2 = 0.50;
1071: tl->omega3 = 1.00;
1072: tl->omega4 = 2.00;
1073: tl->omega5 = 4.00;
1075: /* Default values for trust-region radius update based on reduction */
1076: tl->eta1 = 1.0e-4;
1077: tl->eta2 = 0.25;
1078: tl->eta3 = 0.50;
1079: tl->eta4 = 0.90;
1081: tl->alpha1 = 0.25;
1082: tl->alpha2 = 0.50;
1083: tl->alpha3 = 1.00;
1084: tl->alpha4 = 2.00;
1085: tl->alpha5 = 4.00;
1087: /* Default values for trust-region radius update based on interpolation */
1088: tl->mu1 = 0.10;
1089: tl->mu2 = 0.50;
1091: tl->gamma1 = 0.25;
1092: tl->gamma2 = 0.50;
1093: tl->gamma3 = 2.00;
1094: tl->gamma4 = 4.00;
1096: tl->theta = 0.05;
1098: /* Default values for trust region initialization based on interpolation */
1099: tl->mu1_i = 0.35;
1100: tl->mu2_i = 0.50;
1102: tl->gamma1_i = 0.0625;
1103: tl->gamma2_i = 0.5;
1104: tl->gamma3_i = 2.0;
1105: tl->gamma4_i = 5.0;
1106:
1107: tl->theta_i = 0.25;
1109: /* Remaining parameters */
1110: tl->min_radius = 1.0e-10;
1111: tl->max_radius = 1.0e10;
1112: tl->epsilon = 1.0e-6;
1114: tl->ksp_type = NTL_KSP_STCG;
1115: tl->pc_type = NTL_PC_BFGS;
1116: tl->bfgs_scale_type = BFGS_SCALE_AHESS;
1117: tl->init_type = NTL_INIT_INTERPOLATION;
1118: tl->update_type = NTL_UPDATE_REDUCTION;
1120: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
1121: TaoLineSearchSetType(tao->linesearch, morethuente_type);
1122: TaoLineSearchUseTaoSolverRoutines(tao->linesearch, tao);
1124: KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
1126: return(0);
1127: }