Actual source code: bnk.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/bound/impls/bnk/bnk.h>
3: #include <petscksp.h>
5: static const char *BNK_INIT[64] = {"constant", "direction", "interpolation"};
6: static const char *BNK_UPDATE[64] = {"step", "reduction", "interpolation"};
7: static const char *BNK_AS[64] = {"none", "bertsekas"};
9: /*------------------------------------------------------------*/
11: /* Routine for initializing the KSP solver, the BFGS preconditioner, and the initial trust radius estimation */
13: PetscErrorCode TaoBNKInitialize(Tao tao, PetscInt initType, PetscBool *needH)
14: {
15: PetscErrorCode ierr;
16: TAO_BNK *bnk = (TAO_BNK *)tao->data;
17: PC pc;
19: PetscReal f_min, ftrial, prered, actred, kappa, sigma, resnorm;
20: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
21: PetscBool is_bfgs, is_jacobi, is_symmetric, sym_set;
22: PetscInt n, N, nDiff;
23: PetscInt i_max = 5;
24: PetscInt j_max = 1;
25: PetscInt i, j;
28: /* Project the current point onto the feasible set */
29: TaoComputeVariableBounds(tao);
30: TaoSetVariableBounds(bnk->bncg, tao->XL, tao->XU);
31: if (tao->bounded) {
32: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
33: }
35: /* Project the initial point onto the feasible region */
36: TaoBoundSolution(tao->solution, tao->XL,tao->XU, 0.0, &nDiff, tao->solution);
38: /* Check convergence criteria */
39: TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);
40: TaoBNKEstimateActiveSet(tao, bnk->as_type);
41: VecCopy(bnk->unprojected_gradient, tao->gradient);
42: VecISSet(tao->gradient, bnk->active_idx, 0.0);
43: TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);
45: /* Test the initial point for convergence */
46: VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->W);
47: VecNorm(bnk->W, NORM_2, &resnorm);
48: if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(resnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
49: TaoLogConvergenceHistory(tao,bnk->f,resnorm,0.0,tao->ksp_its);
50: TaoMonitor(tao,tao->niter,bnk->f,resnorm,0.0,1.0);
51: (*tao->ops->convergencetest)(tao,tao->cnvP);
52: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
54: /* Reset KSP stopping reason counters */
55: bnk->ksp_atol = 0;
56: bnk->ksp_rtol = 0;
57: bnk->ksp_dtol = 0;
58: bnk->ksp_ctol = 0;
59: bnk->ksp_negc = 0;
60: bnk->ksp_iter = 0;
61: bnk->ksp_othr = 0;
63: /* Reset accepted step type counters */
64: bnk->tot_cg_its = 0;
65: bnk->newt = 0;
66: bnk->bfgs = 0;
67: bnk->sgrad = 0;
68: bnk->grad = 0;
70: /* Initialize the Hessian perturbation */
71: bnk->pert = bnk->sval;
73: /* Reset initial steplength to zero (this helps BNCG reset its direction internally) */
74: VecSet(tao->stepdirection, 0.0);
76: /* Allocate the vectors needed for the BFGS approximation */
77: KSPGetPC(tao->ksp, &pc);
78: PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);
79: PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);
80: if (is_bfgs) {
81: bnk->bfgs_pre = pc;
82: PCLMVMGetMatLMVM(bnk->bfgs_pre, &bnk->M);
83: VecGetLocalSize(tao->solution, &n);
84: VecGetSize(tao->solution, &N);
85: MatSetSizes(bnk->M, n, n, N, N);
86: MatLMVMAllocate(bnk->M, tao->solution, bnk->unprojected_gradient);
87: MatIsSymmetricKnown(bnk->M, &sym_set, &is_symmetric);
88: if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
89: } else if (is_jacobi) {
90: PCJacobiSetUseAbs(pc,PETSC_TRUE);
91: }
93: /* Prepare the min/max vectors for safeguarding diagonal scales */
94: VecSet(bnk->Diag_min, bnk->dmin);
95: VecSet(bnk->Diag_max, bnk->dmax);
97: /* Initialize trust-region radius. The initialization is only performed
98: when we are using Nash, Steihaug-Toint or the Generalized Lanczos method. */
99: *needH = PETSC_TRUE;
100: if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
101: switch(initType) {
102: case BNK_INIT_CONSTANT:
103: /* Use the initial radius specified */
104: tao->trust = tao->trust0;
105: break;
107: case BNK_INIT_INTERPOLATION:
108: /* Use interpolation based on the initial Hessian */
109: max_radius = 0.0;
110: tao->trust = tao->trust0;
111: for (j = 0; j < j_max; ++j) {
112: f_min = bnk->f;
113: sigma = 0.0;
115: if (*needH) {
116: /* Compute the Hessian at the new step, and extract the inactive subsystem */
117: bnk->computehessian(tao);
118: TaoBNKEstimateActiveSet(tao, BNK_AS_NONE);
119: MatDestroy(&bnk->H_inactive);
120: if (bnk->active_idx) {
121: MatCreateSubMatrix(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->H_inactive);
122: } else {
123: PetscObjectReference((PetscObject)tao->hessian);
124: bnk->H_inactive = tao->hessian;
125: }
126: *needH = PETSC_FALSE;
127: }
129: for (i = 0; i < i_max; ++i) {
130: /* Take a steepest descent step and snap it to bounds */
131: VecCopy(tao->solution, bnk->Xold);
132: VecAXPY(tao->solution, -tao->trust/bnk->gnorm, tao->gradient);
133: TaoBoundSolution(tao->solution, tao->XL,tao->XU, 0.0, &nDiff, tao->solution);
134: /* Compute the step we actually accepted */
135: VecCopy(tao->solution, bnk->W);
136: VecAXPY(bnk->W, -1.0, bnk->Xold);
137: /* Compute the objective at the trial */
138: TaoComputeObjective(tao, tao->solution, &ftrial);
139: if (PetscIsInfOrNanReal(bnk->f)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
140: VecCopy(bnk->Xold, tao->solution);
141: if (PetscIsInfOrNanReal(ftrial)) {
142: tau = bnk->gamma1_i;
143: } else {
144: if (ftrial < f_min) {
145: f_min = ftrial;
146: sigma = -tao->trust / bnk->gnorm;
147: }
149: /* Compute the predicted and actual reduction */
150: if (bnk->active_idx) {
151: VecGetSubVector(bnk->W, bnk->inactive_idx, &bnk->X_inactive);
152: VecGetSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work);
153: } else {
154: bnk->X_inactive = bnk->W;
155: bnk->inactive_work = bnk->Xwork;
156: }
157: MatMult(bnk->H_inactive, bnk->X_inactive, bnk->inactive_work);
158: VecDot(bnk->X_inactive, bnk->inactive_work, &prered);
159: if (bnk->active_idx) {
160: VecRestoreSubVector(bnk->W, bnk->inactive_idx, &bnk->X_inactive);
161: VecRestoreSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work);
162: }
163: prered = tao->trust * (bnk->gnorm - 0.5 * tao->trust * prered / (bnk->gnorm * bnk->gnorm));
164: actred = bnk->f - ftrial;
165: if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) {
166: kappa = 1.0;
167: } else {
168: kappa = actred / prered;
169: }
171: tau_1 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust + (1.0 - bnk->theta_i) * prered - actred);
172: tau_2 = bnk->theta_i * bnk->gnorm * tao->trust / (bnk->theta_i * bnk->gnorm * tao->trust - (1.0 + bnk->theta_i) * prered + actred);
173: tau_min = PetscMin(tau_1, tau_2);
174: tau_max = PetscMax(tau_1, tau_2);
176: if (PetscAbsScalar(kappa - (PetscReal)1.0) <= bnk->mu1_i) {
177: /* Great agreement */
178: max_radius = PetscMax(max_radius, tao->trust);
180: if (tau_max < 1.0) {
181: tau = bnk->gamma3_i;
182: } else if (tau_max > bnk->gamma4_i) {
183: tau = bnk->gamma4_i;
184: } else {
185: tau = tau_max;
186: }
187: } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= bnk->mu2_i) {
188: /* Good agreement */
189: max_radius = PetscMax(max_radius, tao->trust);
191: if (tau_max < bnk->gamma2_i) {
192: tau = bnk->gamma2_i;
193: } else if (tau_max > bnk->gamma3_i) {
194: tau = bnk->gamma3_i;
195: } else {
196: tau = tau_max;
197: }
198: } else {
199: /* Not good agreement */
200: if (tau_min > 1.0) {
201: tau = bnk->gamma2_i;
202: } else if (tau_max < bnk->gamma1_i) {
203: tau = bnk->gamma1_i;
204: } else if ((tau_min < bnk->gamma1_i) && (tau_max >= 1.0)) {
205: tau = bnk->gamma1_i;
206: } else if ((tau_1 >= bnk->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1_i) || (tau_2 >= 1.0))) {
207: tau = tau_1;
208: } else if ((tau_2 >= bnk->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1_i) || (tau_2 >= 1.0))) {
209: tau = tau_2;
210: } else {
211: tau = tau_max;
212: }
213: }
214: }
215: tao->trust = tau * tao->trust;
216: }
218: if (f_min < bnk->f) {
219: /* We accidentally found a solution better than the initial, so accept it */
220: bnk->f = f_min;
221: VecCopy(tao->solution, bnk->Xold);
222: VecAXPY(tao->solution,sigma,tao->gradient);
223: TaoBoundSolution(tao->solution, tao->XL,tao->XU, 0.0, &nDiff, tao->solution);
224: VecCopy(tao->solution, tao->stepdirection);
225: VecAXPY(tao->stepdirection, -1.0, bnk->Xold);
226: TaoComputeGradient(tao,tao->solution,bnk->unprojected_gradient);
227: TaoBNKEstimateActiveSet(tao, bnk->as_type);
228: VecCopy(bnk->unprojected_gradient, tao->gradient);
229: VecISSet(tao->gradient, bnk->active_idx, 0.0);
230: /* Compute gradient at the new iterate and flip switch to compute the Hessian later */
231: TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);
232: *needH = PETSC_TRUE;
233: /* Test the new step for convergence */
234: VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->W);
235: VecNorm(bnk->W, NORM_2, &resnorm);
236: if (PetscIsInfOrNanReal(resnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
237: TaoLogConvergenceHistory(tao,bnk->f,resnorm,0.0,tao->ksp_its);
238: TaoMonitor(tao,tao->niter,bnk->f,resnorm,0.0,1.0);
239: (*tao->ops->convergencetest)(tao,tao->cnvP);
240: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
241: /* active BNCG recycling early because we have a stepdirection computed */
242: TaoSetRecycleHistory(bnk->bncg, PETSC_TRUE);
243: }
244: }
245: tao->trust = PetscMax(tao->trust, max_radius);
247: /* Ensure that the trust radius is within the limits */
248: tao->trust = PetscMax(tao->trust, bnk->min_radius);
249: tao->trust = PetscMin(tao->trust, bnk->max_radius);
250: break;
252: default:
253: /* Norm of the first direction will initialize radius */
254: tao->trust = 0.0;
255: break;
256: }
257: }
258: return(0);
259: }
261: /*------------------------------------------------------------*/
263: /* Routine for computing the exact Hessian and preparing the preconditioner at the new iterate */
265: PetscErrorCode TaoBNKComputeHessian(Tao tao)
266: {
267: PetscErrorCode ierr;
268: TAO_BNK *bnk = (TAO_BNK *)tao->data;
271: /* Compute the Hessian */
272: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
273: /* Add a correction to the BFGS preconditioner */
274: if (bnk->M) {
275: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
276: }
277: /* Prepare the reduced sub-matrices for the inactive set */
278: if (bnk->Hpre_inactive) {
279: MatDestroy(&bnk->Hpre_inactive);
280: }
281: if (bnk->H_inactive) {
282: MatDestroy(&bnk->H_inactive);
283: }
284: if (bnk->active_idx) {
285: MatCreateSubMatrix(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->H_inactive);
286: if (tao->hessian == tao->hessian_pre) {
287: PetscObjectReference((PetscObject)bnk->H_inactive);
288: bnk->Hpre_inactive = bnk->H_inactive;
289: } else {
290: MatCreateSubMatrix(tao->hessian_pre, bnk->inactive_idx, bnk->inactive_idx, MAT_INITIAL_MATRIX, &bnk->Hpre_inactive);
291: }
292: if (bnk->bfgs_pre) {
293: PCLMVMSetIS(bnk->bfgs_pre, bnk->inactive_idx);
294: }
295: } else {
296: PetscObjectReference((PetscObject)tao->hessian);
297: bnk->H_inactive = tao->hessian;
298: if (tao->hessian == tao->hessian_pre) {
299: PetscObjectReference((PetscObject)bnk->H_inactive);
300: bnk->Hpre_inactive = bnk->H_inactive;
301: } else {
302: PetscObjectReference((PetscObject)tao->hessian_pre);
303: bnk->Hpre_inactive = tao->hessian_pre;
304: }
305: if (bnk->bfgs_pre) {
306: PCLMVMClearIS(bnk->bfgs_pre);
307: }
308: }
309: return(0);
310: }
312: /*------------------------------------------------------------*/
314: /* Routine for estimating the active set */
316: PetscErrorCode TaoBNKEstimateActiveSet(Tao tao, PetscInt asType)
317: {
318: PetscErrorCode ierr;
319: TAO_BNK *bnk = (TAO_BNK *)tao->data;
320: PetscBool hessComputed, diagExists;
323: switch (asType) {
324: case BNK_AS_NONE:
325: ISDestroy(&bnk->inactive_idx);
326: VecWhichInactive(tao->XL, tao->solution, bnk->unprojected_gradient, tao->XU, PETSC_TRUE, &bnk->inactive_idx);
327: ISDestroy(&bnk->active_idx);
328: ISComplementVec(bnk->inactive_idx, tao->solution, &bnk->active_idx);
329: break;
331: case BNK_AS_BERTSEKAS:
332: /* Compute the trial step vector with which we will estimate the active set at the next iteration */
333: if (bnk->M) {
334: /* If the BFGS preconditioner matrix is available, we will construct a trial step with it */
335: MatSolve(bnk->M, bnk->unprojected_gradient, bnk->W);
336: } else {
337: hessComputed = diagExists = PETSC_FALSE;
338: if (tao->hessian) {
339: MatAssembled(tao->hessian, &hessComputed);
340: }
341: if (hessComputed) {
342: MatHasOperation(tao->hessian, MATOP_GET_DIAGONAL, &diagExists);
343: }
344: if (diagExists) {
345: /* BFGS preconditioner doesn't exist so let's invert the absolute diagonal of the Hessian instead onto the gradient */
346: MatGetDiagonal(tao->hessian, bnk->Xwork);
347: VecAbs(bnk->Xwork);
348: VecMedian(bnk->Diag_min, bnk->Xwork, bnk->Diag_max, bnk->Xwork);
349: VecReciprocal(bnk->Xwork);
350: VecPointwiseMult(bnk->W, bnk->Xwork, bnk->unprojected_gradient);
351: } else {
352: /* If the Hessian or its diagonal does not exist, we will simply use gradient step */
353: VecCopy(bnk->unprojected_gradient, bnk->W);
354: }
355: }
356: VecScale(bnk->W, -1.0);
357: TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, bnk->unprojected_gradient, bnk->W, bnk->Xwork, bnk->as_step, &bnk->as_tol,
358: &bnk->active_lower, &bnk->active_upper, &bnk->active_fixed, &bnk->active_idx, &bnk->inactive_idx);
359: break;
361: default:
362: break;
363: }
364: return(0);
365: }
367: /*------------------------------------------------------------*/
369: /* Routine for bounding the step direction */
371: PetscErrorCode TaoBNKBoundStep(Tao tao, PetscInt asType, Vec step)
372: {
373: PetscErrorCode ierr;
374: TAO_BNK *bnk = (TAO_BNK *)tao->data;
377: switch (asType) {
378: case BNK_AS_NONE:
379: VecISSet(step, bnk->active_idx, 0.0);
380: break;
382: case BNK_AS_BERTSEKAS:
383: TaoBoundStep(tao->solution, tao->XL, tao->XU, bnk->active_lower, bnk->active_upper, bnk->active_fixed, 1.0, step);
384: break;
386: default:
387: break;
388: }
389: return(0);
390: }
392: /*------------------------------------------------------------*/
394: /* Routine for taking a finite number of BNCG iterations to
395: accelerate Newton convergence.
397: In practice, this approach simply trades off Hessian evaluations
398: for more gradient evaluations.
399: */
401: PetscErrorCode TaoBNKTakeCGSteps(Tao tao, PetscBool *terminate)
402: {
403: TAO_BNK *bnk = (TAO_BNK *)tao->data;
404: PetscErrorCode ierr;
407: *terminate = PETSC_FALSE;
408: if (bnk->max_cg_its > 0) {
409: /* Copy the current function value (important vectors are already shared) */
410: bnk->bncg_ctx->f = bnk->f;
411: /* Take some small finite number of BNCG iterations */
412: TaoSolve(bnk->bncg);
413: /* Add the number of gradient and function evaluations to the total */
414: tao->nfuncs += bnk->bncg->nfuncs;
415: tao->nfuncgrads += bnk->bncg->nfuncgrads;
416: tao->ngrads += bnk->bncg->ngrads;
417: tao->nhess += bnk->bncg->nhess;
418: bnk->tot_cg_its += bnk->bncg->niter;
419: /* Extract the BNCG function value out and save it into BNK */
420: bnk->f = bnk->bncg_ctx->f;
421: if (bnk->bncg->reason == TAO_CONVERGED_GATOL || bnk->bncg->reason == TAO_CONVERGED_GRTOL || bnk->bncg->reason == TAO_CONVERGED_GTTOL || bnk->bncg->reason == TAO_CONVERGED_MINF) {
422: *terminate = PETSC_TRUE;
423: } else {
424: TaoBNKEstimateActiveSet(tao, bnk->as_type);
425: }
426: }
427: return(0);
428: }
430: /*------------------------------------------------------------*/
432: /* Routine for computing the Newton step. */
434: PetscErrorCode TaoBNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type)
435: {
436: PetscErrorCode ierr;
437: TAO_BNK *bnk = (TAO_BNK *)tao->data;
438: PetscInt bfgsUpdates = 0;
439: PetscInt kspits;
440: PetscBool is_lmvm;
443: /* If there are no inactive variables left, save some computation and return an adjusted zero step
444: that has (l-x) and (u-x) for lower and upper bounded variables. */
445: if (!bnk->inactive_idx) {
446: VecSet(tao->stepdirection, 0.0);
447: TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
448: return(0);
449: }
451: /* Shift the reduced Hessian matrix */
452: if ((shift) && (bnk->pert > 0)) {
453: PetscObjectTypeCompare((PetscObject)tao->hessian, MATLMVM, &is_lmvm);
454: if (is_lmvm) {
455: MatShift(tao->hessian, bnk->pert);
456: } else {
457: MatShift(bnk->H_inactive, bnk->pert);
458: if (bnk->H_inactive != bnk->Hpre_inactive) {
459: MatShift(bnk->Hpre_inactive, bnk->pert);
460: }
461: }
462: }
464: /* Solve the Newton system of equations */
465: tao->ksp_its = 0;
466: VecSet(tao->stepdirection, 0.0);
467: KSPReset(tao->ksp);
468: KSPResetFromOptions(tao->ksp);
469: KSPSetOperators(tao->ksp,bnk->H_inactive,bnk->Hpre_inactive);
470: VecCopy(bnk->unprojected_gradient, bnk->Gwork);
471: if (bnk->active_idx) {
472: VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive);
473: VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive);
474: } else {
475: bnk->G_inactive = bnk->unprojected_gradient;
476: bnk->X_inactive = tao->stepdirection;
477: }
478: if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
479: KSPCGSetRadius(tao->ksp,tao->trust);
480: KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive);
481: KSPGetIterationNumber(tao->ksp,&kspits);
482: tao->ksp_its+=kspits;
483: tao->ksp_tot_its+=kspits;
484: KSPCGGetNormD(tao->ksp,&bnk->dnorm);
486: if (0.0 == tao->trust) {
487: /* Radius was uninitialized; use the norm of the direction */
488: if (bnk->dnorm > 0.0) {
489: tao->trust = bnk->dnorm;
491: /* Modify the radius if it is too large or small */
492: tao->trust = PetscMax(tao->trust, bnk->min_radius);
493: tao->trust = PetscMin(tao->trust, bnk->max_radius);
494: } else {
495: /* The direction was bad; set radius to default value and re-solve
496: the trust-region subproblem to get a direction */
497: tao->trust = tao->trust0;
499: /* Modify the radius if it is too large or small */
500: tao->trust = PetscMax(tao->trust, bnk->min_radius);
501: tao->trust = PetscMin(tao->trust, bnk->max_radius);
503: KSPCGSetRadius(tao->ksp,tao->trust);
504: KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive);
505: KSPGetIterationNumber(tao->ksp,&kspits);
506: tao->ksp_its+=kspits;
507: tao->ksp_tot_its+=kspits;
508: KSPCGGetNormD(tao->ksp,&bnk->dnorm);
510: if (bnk->dnorm == 0.0) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_PLIB, "Initial direction zero");
511: }
512: }
513: } else {
514: KSPSolve(tao->ksp, bnk->G_inactive, bnk->X_inactive);
515: KSPGetIterationNumber(tao->ksp, &kspits);
516: tao->ksp_its += kspits;
517: tao->ksp_tot_its+=kspits;
518: }
519: /* Restore sub vectors back */
520: if (bnk->active_idx) {
521: VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive);
522: VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive);
523: }
524: /* Make sure the safeguarded fall-back step is zero for actively bounded variables */
525: VecScale(tao->stepdirection, -1.0);
526: TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
528: /* Record convergence reasons */
529: KSPGetConvergedReason(tao->ksp, ksp_reason);
530: if (KSP_CONVERGED_ATOL == *ksp_reason) {
531: ++bnk->ksp_atol;
532: } else if (KSP_CONVERGED_RTOL == *ksp_reason) {
533: ++bnk->ksp_rtol;
534: } else if (KSP_CONVERGED_CG_CONSTRAINED == *ksp_reason) {
535: ++bnk->ksp_ctol;
536: } else if (KSP_CONVERGED_CG_NEG_CURVE == *ksp_reason) {
537: ++bnk->ksp_negc;
538: } else if (KSP_DIVERGED_DTOL == *ksp_reason) {
539: ++bnk->ksp_dtol;
540: } else if (KSP_DIVERGED_ITS == *ksp_reason) {
541: ++bnk->ksp_iter;
542: } else {
543: ++bnk->ksp_othr;
544: }
546: /* Make sure the BFGS preconditioner is healthy */
547: if (bnk->M) {
548: MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates);
549: if ((KSP_DIVERGED_INDEFINITE_PC == *ksp_reason) && (bfgsUpdates > 0)) {
550: /* Preconditioner is numerically indefinite; reset the approximation. */
551: MatLMVMReset(bnk->M, PETSC_FALSE);
552: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
553: }
554: }
555: *step_type = BNK_NEWTON;
556: return(0);
557: }
559: /*------------------------------------------------------------*/
561: /* Routine for recomputing the predicted reduction for a given step vector */
563: PetscErrorCode TaoBNKRecomputePred(Tao tao, Vec S, PetscReal *prered)
564: {
565: PetscErrorCode ierr;
566: TAO_BNK *bnk = (TAO_BNK *)tao->data;
569: /* Extract subvectors associated with the inactive set */
570: if (bnk->active_idx) {
571: VecGetSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive);
572: VecGetSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work);
573: VecGetSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive);
574: } else {
575: bnk->X_inactive = tao->stepdirection;
576: bnk->inactive_work = bnk->Xwork;
577: bnk->G_inactive = bnk->Gwork;
578: }
579: /* Recompute the predicted decrease based on the quadratic model */
580: MatMult(bnk->H_inactive, bnk->X_inactive, bnk->inactive_work);
581: VecAYPX(bnk->inactive_work, -0.5, bnk->G_inactive);
582: VecDot(bnk->inactive_work, bnk->X_inactive, prered);
583: /* Restore the sub vectors */
584: if (bnk->active_idx) {
585: VecRestoreSubVector(tao->stepdirection, bnk->inactive_idx, &bnk->X_inactive);
586: VecRestoreSubVector(bnk->Xwork, bnk->inactive_idx, &bnk->inactive_work);
587: VecRestoreSubVector(bnk->Gwork, bnk->inactive_idx, &bnk->G_inactive);
588: }
589: return(0);
590: }
592: /*------------------------------------------------------------*/
594: /* Routine for ensuring that the Newton step is a descent direction.
596: The step direction falls back onto BFGS, scaled gradient and gradient steps
597: in the event that the Newton step fails the test.
598: */
600: PetscErrorCode TaoBNKSafeguardStep(Tao tao, KSPConvergedReason ksp_reason, PetscInt *stepType)
601: {
602: PetscErrorCode ierr;
603: TAO_BNK *bnk = (TAO_BNK *)tao->data;
605: PetscReal gdx, e_min;
606: PetscInt bfgsUpdates;
609: switch (*stepType) {
610: case BNK_NEWTON:
611: VecDot(tao->stepdirection, tao->gradient, &gdx);
612: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
613: /* Newton step is not descent or direction produced Inf or NaN
614: Update the perturbation for next time */
615: if (bnk->pert <= 0.0) {
616: /* Initialize the perturbation */
617: bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
618: if (bnk->is_gltr) {
619: KSPGLTRGetMinEig(tao->ksp,&e_min);
620: bnk->pert = PetscMax(bnk->pert, -e_min);
621: }
622: } else {
623: /* Increase the perturbation */
624: bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
625: }
627: if (!bnk->M) {
628: /* We don't have the bfgs matrix around and updated
629: Must use gradient direction in this case */
630: VecCopy(tao->gradient, tao->stepdirection);
631: *stepType = BNK_GRADIENT;
632: } else {
633: /* Attempt to use the BFGS direction */
634: MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);
636: /* Check for success (descent direction)
637: NOTE: Negative gdx here means not a descent direction because
638: the fall-back step is missing a negative sign. */
639: VecDot(tao->gradient, tao->stepdirection, &gdx);
640: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
641: /* BFGS direction is not descent or direction produced not a number
642: We can assert bfgsUpdates > 1 in this case because
643: the first solve produces the scaled gradient direction,
644: which is guaranteed to be descent */
646: /* Use steepest descent direction (scaled) */
647: MatLMVMReset(bnk->M, PETSC_FALSE);
648: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
649: MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);
651: *stepType = BNK_SCALED_GRADIENT;
652: } else {
653: MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates);
654: if (1 == bfgsUpdates) {
655: /* The first BFGS direction is always the scaled gradient */
656: *stepType = BNK_SCALED_GRADIENT;
657: } else {
658: *stepType = BNK_BFGS;
659: }
660: }
661: }
662: /* Make sure the safeguarded fall-back step is zero for actively bounded variables */
663: VecScale(tao->stepdirection, -1.0);
664: TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
665: } else {
666: /* Computed Newton step is descent */
667: switch (ksp_reason) {
668: case KSP_DIVERGED_NANORINF:
669: case KSP_DIVERGED_BREAKDOWN:
670: case KSP_DIVERGED_INDEFINITE_MAT:
671: case KSP_DIVERGED_INDEFINITE_PC:
672: case KSP_CONVERGED_CG_NEG_CURVE:
673: /* Matrix or preconditioner is indefinite; increase perturbation */
674: if (bnk->pert <= 0.0) {
675: /* Initialize the perturbation */
676: bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
677: if (bnk->is_gltr) {
678: KSPGLTRGetMinEig(tao->ksp, &e_min);
679: bnk->pert = PetscMax(bnk->pert, -e_min);
680: }
681: } else {
682: /* Increase the perturbation */
683: bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
684: }
685: break;
687: default:
688: /* Newton step computation is good; decrease perturbation */
689: bnk->pert = PetscMin(bnk->psfac * bnk->pert, bnk->pmsfac * bnk->gnorm);
690: if (bnk->pert < bnk->pmin) {
691: bnk->pert = 0.0;
692: }
693: break;
694: }
695: *stepType = BNK_NEWTON;
696: }
697: break;
699: case BNK_BFGS:
700: /* Check for success (descent direction) */
701: VecDot(tao->stepdirection, tao->gradient, &gdx);
702: if (gdx >= 0 || PetscIsInfOrNanReal(gdx)) {
703: /* Step is not descent or solve was not successful
704: Use steepest descent direction (scaled) */
705: MatLMVMReset(bnk->M, PETSC_FALSE);
706: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
707: MatSolve(bnk->M, tao->gradient, tao->stepdirection);
708: VecScale(tao->stepdirection,-1.0);
709: TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
710: *stepType = BNK_SCALED_GRADIENT;
711: } else {
712: *stepType = BNK_BFGS;
713: }
714: break;
716: case BNK_SCALED_GRADIENT:
717: break;
719: default:
720: break;
721: }
723: return(0);
724: }
726: /*------------------------------------------------------------*/
728: /* Routine for performing a bound-projected More-Thuente line search.
730: Includes fallbacks to BFGS, scaled gradient, and unscaled gradient steps if the
731: Newton step does not produce a valid step length.
732: */
734: PetscErrorCode TaoBNKPerformLineSearch(Tao tao, PetscInt *stepType, PetscReal *steplen, TaoLineSearchConvergedReason *reason)
735: {
736: TAO_BNK *bnk = (TAO_BNK *)tao->data;
738: TaoLineSearchConvergedReason ls_reason;
740: PetscReal e_min, gdx;
741: PetscInt bfgsUpdates;
744: /* Perform the linesearch */
745: TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason);
746: TaoAddLineSearchCounts(tao);
748: while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && *stepType != BNK_SCALED_GRADIENT && *stepType != BNK_GRADIENT) {
749: /* Linesearch failed, revert solution */
750: bnk->f = bnk->fold;
751: VecCopy(bnk->Xold, tao->solution);
752: VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);
754: switch(*stepType) {
755: case BNK_NEWTON:
756: /* Failed to obtain acceptable iterate with Newton step
757: Update the perturbation for next time */
758: if (bnk->pert <= 0.0) {
759: /* Initialize the perturbation */
760: bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
761: if (bnk->is_gltr) {
762: KSPGLTRGetMinEig(tao->ksp,&e_min);
763: bnk->pert = PetscMax(bnk->pert, -e_min);
764: }
765: } else {
766: /* Increase the perturbation */
767: bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
768: }
770: if (!bnk->M) {
771: /* We don't have the bfgs matrix around and being updated
772: Must use gradient direction in this case */
773: VecCopy(bnk->unprojected_gradient, tao->stepdirection);
774: *stepType = BNK_GRADIENT;
775: } else {
776: /* Attempt to use the BFGS direction */
777: MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);
778: /* Check for success (descent direction)
779: NOTE: Negative gdx means not a descent direction because the step here is missing a negative sign. */
780: VecDot(tao->gradient, tao->stepdirection, &gdx);
781: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
782: /* BFGS direction is not descent or direction produced not a number
783: We can assert bfgsUpdates > 1 in this case
784: Use steepest descent direction (scaled) */
785: MatLMVMReset(bnk->M, PETSC_FALSE);
786: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
787: MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);
789: bfgsUpdates = 1;
790: *stepType = BNK_SCALED_GRADIENT;
791: } else {
792: MatLMVMGetUpdateCount(bnk->M, &bfgsUpdates);
793: if (1 == bfgsUpdates) {
794: /* The first BFGS direction is always the scaled gradient */
795: *stepType = BNK_SCALED_GRADIENT;
796: } else {
797: *stepType = BNK_BFGS;
798: }
799: }
800: }
801: break;
803: case BNK_BFGS:
804: /* Can only enter if pc_type == BNK_PC_BFGS
805: Failed to obtain acceptable iterate with BFGS step
806: Attempt to use the scaled gradient direction */
807: MatLMVMReset(bnk->M, PETSC_FALSE);
808: MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);
809: MatSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);
811: bfgsUpdates = 1;
812: *stepType = BNK_SCALED_GRADIENT;
813: break;
814: }
815: /* Make sure the safeguarded fall-back step is zero for actively bounded variables */
816: VecScale(tao->stepdirection, -1.0);
817: TaoBNKBoundStep(tao, bnk->as_type, tao->stepdirection);
819: /* Perform one last line search with the fall-back step */
820: TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, bnk->unprojected_gradient, tao->stepdirection, steplen, &ls_reason);
821: TaoAddLineSearchCounts(tao);
822: }
823: *reason = ls_reason;
824: return(0);
825: }
827: /*------------------------------------------------------------*/
829: /* Routine for updating the trust radius.
831: Function features three different update methods:
832: 1) Line-search step length based
833: 2) Predicted decrease on the CG quadratic model
834: 3) Interpolation
835: */
837: PetscErrorCode TaoBNKUpdateTrustRadius(Tao tao, PetscReal prered, PetscReal actred, PetscInt updateType, PetscInt stepType, PetscBool *accept)
838: {
839: TAO_BNK *bnk = (TAO_BNK *)tao->data;
842: PetscReal step, kappa;
843: PetscReal gdx, tau_1, tau_2, tau_min, tau_max;
846: /* Update trust region radius */
847: *accept = PETSC_FALSE;
848: switch(updateType) {
849: case BNK_UPDATE_STEP:
850: *accept = PETSC_TRUE; /* always accept here because line search succeeded */
851: if (stepType == BNK_NEWTON) {
852: TaoLineSearchGetStepLength(tao->linesearch, &step);
853: if (step < bnk->nu1) {
854: /* Very bad step taken; reduce radius */
855: tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust);
856: } else if (step < bnk->nu2) {
857: /* Reasonably bad step taken; reduce radius */
858: tao->trust = bnk->omega2 * PetscMin(bnk->dnorm, tao->trust);
859: } else if (step < bnk->nu3) {
860: /* Reasonable step was taken; leave radius alone */
861: if (bnk->omega3 < 1.0) {
862: tao->trust = bnk->omega3 * PetscMin(bnk->dnorm, tao->trust);
863: } else if (bnk->omega3 > 1.0) {
864: tao->trust = PetscMax(bnk->omega3 * bnk->dnorm, tao->trust);
865: }
866: } else if (step < bnk->nu4) {
867: /* Full step taken; increase the radius */
868: tao->trust = PetscMax(bnk->omega4 * bnk->dnorm, tao->trust);
869: } else {
870: /* More than full step taken; increase the radius */
871: tao->trust = PetscMax(bnk->omega5 * bnk->dnorm, tao->trust);
872: }
873: } else {
874: /* Newton step was not good; reduce the radius */
875: tao->trust = bnk->omega1 * PetscMin(bnk->dnorm, tao->trust);
876: }
877: break;
879: case BNK_UPDATE_REDUCTION:
880: if (stepType == BNK_NEWTON) {
881: if ((prered < 0.0) || PetscIsInfOrNanReal(prered)) {
882: /* The predicted reduction has the wrong sign. This cannot
883: happen in infinite precision arithmetic. Step should
884: be rejected! */
885: tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
886: } else {
887: if (PetscIsInfOrNanReal(actred)) {
888: tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
889: } else {
890: if ((PetscAbsScalar(actred) <= PetscMax(1.0, PetscAbsScalar(bnk->f))*bnk->epsilon) && (PetscAbsScalar(prered) <= PetscMax(1.0, PetscAbsScalar(bnk->f))*bnk->epsilon)) {
891: kappa = 1.0;
892: } else {
893: kappa = actred / prered;
894: }
895: /* Accept or reject the step and update radius */
896: if (kappa < bnk->eta1) {
897: /* Reject the step */
898: tao->trust = bnk->alpha1 * PetscMin(tao->trust, bnk->dnorm);
899: } else {
900: /* Accept the step */
901: *accept = PETSC_TRUE;
902: /* Update the trust region radius only if the computed step is at the trust radius boundary */
903: if (bnk->dnorm == tao->trust) {
904: if (kappa < bnk->eta2) {
905: /* Marginal bad step */
906: tao->trust = bnk->alpha2 * tao->trust;
907: } else if (kappa < bnk->eta3) {
908: /* Reasonable step */
909: tao->trust = bnk->alpha3 * tao->trust;
910: } else if (kappa < bnk->eta4) {
911: /* Good step */
912: tao->trust = bnk->alpha4 * tao->trust;
913: } else {
914: /* Very good step */
915: tao->trust = bnk->alpha5 * tao->trust;
916: }
917: }
918: }
919: }
920: }
921: } else {
922: /* Newton step was not good; reduce the radius */
923: tao->trust = bnk->alpha1 * PetscMin(bnk->dnorm, tao->trust);
924: }
925: break;
927: default:
928: if (stepType == BNK_NEWTON) {
929: if (prered < 0.0) {
930: /* The predicted reduction has the wrong sign. This cannot */
931: /* happen in infinite precision arithmetic. Step should */
932: /* be rejected! */
933: tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
934: } else {
935: if (PetscIsInfOrNanReal(actred)) {
936: tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
937: } else {
938: if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) {
939: kappa = 1.0;
940: } else {
941: kappa = actred / prered;
942: }
944: VecDot(tao->gradient, tao->stepdirection, &gdx);
945: tau_1 = bnk->theta * gdx / (bnk->theta * gdx - (1.0 - bnk->theta) * prered + actred);
946: tau_2 = bnk->theta * gdx / (bnk->theta * gdx + (1.0 + bnk->theta) * prered - actred);
947: tau_min = PetscMin(tau_1, tau_2);
948: tau_max = PetscMax(tau_1, tau_2);
950: if (kappa >= 1.0 - bnk->mu1) {
951: /* Great agreement */
952: *accept = PETSC_TRUE;
953: if (tau_max < 1.0) {
954: tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm);
955: } else if (tau_max > bnk->gamma4) {
956: tao->trust = PetscMax(tao->trust, bnk->gamma4 * bnk->dnorm);
957: } else {
958: tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm);
959: }
960: } else if (kappa >= 1.0 - bnk->mu2) {
961: /* Good agreement */
962: *accept = PETSC_TRUE;
963: if (tau_max < bnk->gamma2) {
964: tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm);
965: } else if (tau_max > bnk->gamma3) {
966: tao->trust = PetscMax(tao->trust, bnk->gamma3 * bnk->dnorm);
967: } else if (tau_max < 1.0) {
968: tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm);
969: } else {
970: tao->trust = PetscMax(tao->trust, tau_max * bnk->dnorm);
971: }
972: } else {
973: /* Not good agreement */
974: if (tau_min > 1.0) {
975: tao->trust = bnk->gamma2 * PetscMin(tao->trust, bnk->dnorm);
976: } else if (tau_max < bnk->gamma1) {
977: tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
978: } else if ((tau_min < bnk->gamma1) && (tau_max >= 1.0)) {
979: tao->trust = bnk->gamma1 * PetscMin(tao->trust, bnk->dnorm);
980: } else if ((tau_1 >= bnk->gamma1) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1) || (tau_2 >= 1.0))) {
981: tao->trust = tau_1 * PetscMin(tao->trust, bnk->dnorm);
982: } else if ((tau_2 >= bnk->gamma1) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1) || (tau_2 >= 1.0))) {
983: tao->trust = tau_2 * PetscMin(tao->trust, bnk->dnorm);
984: } else {
985: tao->trust = tau_max * PetscMin(tao->trust, bnk->dnorm);
986: }
987: }
988: }
989: }
990: } else {
991: /* Newton step was not good; reduce the radius */
992: tao->trust = bnk->gamma1 * PetscMin(bnk->dnorm, tao->trust);
993: }
994: break;
995: }
996: /* Make sure the radius does not violate min and max settings */
997: tao->trust = PetscMin(tao->trust, bnk->max_radius);
998: tao->trust = PetscMax(tao->trust, bnk->min_radius);
999: return(0);
1000: }
1002: /* ---------------------------------------------------------- */
1004: PetscErrorCode TaoBNKAddStepCounts(Tao tao, PetscInt stepType)
1005: {
1006: TAO_BNK *bnk = (TAO_BNK *)tao->data;
1009: switch (stepType) {
1010: case BNK_NEWTON:
1011: ++bnk->newt;
1012: break;
1013: case BNK_BFGS:
1014: ++bnk->bfgs;
1015: break;
1016: case BNK_SCALED_GRADIENT:
1017: ++bnk->sgrad;
1018: break;
1019: case BNK_GRADIENT:
1020: ++bnk->grad;
1021: break;
1022: default:
1023: break;
1024: }
1025: return(0);
1026: }
1028: /* ---------------------------------------------------------- */
1030: PetscErrorCode TaoSetUp_BNK(Tao tao)
1031: {
1032: TAO_BNK *bnk = (TAO_BNK *)tao->data;
1034: PetscInt i;
1035: KSPType ksp_type;
1038: if (!tao->gradient) {
1039: VecDuplicate(tao->solution,&tao->gradient);
1040: }
1041: if (!tao->stepdirection) {
1042: VecDuplicate(tao->solution,&tao->stepdirection);
1043: }
1044: if (!bnk->W) {
1045: VecDuplicate(tao->solution,&bnk->W);
1046: }
1047: if (!bnk->Xold) {
1048: VecDuplicate(tao->solution,&bnk->Xold);
1049: }
1050: if (!bnk->Gold) {
1051: VecDuplicate(tao->solution,&bnk->Gold);
1052: }
1053: if (!bnk->Xwork) {
1054: VecDuplicate(tao->solution,&bnk->Xwork);
1055: }
1056: if (!bnk->Gwork) {
1057: VecDuplicate(tao->solution,&bnk->Gwork);
1058: }
1059: if (!bnk->unprojected_gradient) {
1060: VecDuplicate(tao->solution,&bnk->unprojected_gradient);
1061: }
1062: if (!bnk->unprojected_gradient_old) {
1063: VecDuplicate(tao->solution,&bnk->unprojected_gradient_old);
1064: }
1065: if (!bnk->Diag_min) {
1066: VecDuplicate(tao->solution,&bnk->Diag_min);
1067: }
1068: if (!bnk->Diag_max) {
1069: VecDuplicate(tao->solution,&bnk->Diag_max);
1070: }
1071: if (bnk->max_cg_its > 0) {
1072: /* Ensure that the important common vectors are shared between BNK and embedded BNCG */
1073: bnk->bncg_ctx = (TAO_BNCG *)bnk->bncg->data;
1074: PetscObjectReference((PetscObject)(bnk->unprojected_gradient_old));
1075: VecDestroy(&bnk->bncg_ctx->unprojected_gradient_old);
1076: bnk->bncg_ctx->unprojected_gradient_old = bnk->unprojected_gradient_old;
1077: PetscObjectReference((PetscObject)(bnk->unprojected_gradient));
1078: VecDestroy(&bnk->bncg_ctx->unprojected_gradient);
1079: bnk->bncg_ctx->unprojected_gradient = bnk->unprojected_gradient;
1080: PetscObjectReference((PetscObject)(bnk->Gold));
1081: VecDestroy(&bnk->bncg_ctx->G_old);
1082: bnk->bncg_ctx->G_old = bnk->Gold;
1083: PetscObjectReference((PetscObject)(tao->gradient));
1084: VecDestroy(&bnk->bncg->gradient);
1085: bnk->bncg->gradient = tao->gradient;
1086: PetscObjectReference((PetscObject)(tao->stepdirection));
1087: VecDestroy(&bnk->bncg->stepdirection);
1088: bnk->bncg->stepdirection = tao->stepdirection;
1089: TaoSetInitialVector(bnk->bncg, tao->solution);
1090: /* Copy over some settings from BNK into BNCG */
1091: TaoSetMaximumIterations(bnk->bncg, bnk->max_cg_its);
1092: TaoSetTolerances(bnk->bncg, tao->gatol, tao->grtol, tao->gttol);
1093: TaoSetFunctionLowerBound(bnk->bncg, tao->fmin);
1094: TaoSetConvergenceTest(bnk->bncg, tao->ops->convergencetest, tao->cnvP);
1095: TaoSetObjectiveRoutine(bnk->bncg, tao->ops->computeobjective, tao->user_objP);
1096: TaoSetGradientRoutine(bnk->bncg, tao->ops->computegradient, tao->user_gradP);
1097: TaoSetObjectiveAndGradientRoutine(bnk->bncg, tao->ops->computeobjectiveandgradient, tao->user_objgradP);
1098: PetscObjectCopyFortranFunctionPointers((PetscObject)tao, (PetscObject)(bnk->bncg));
1099: for (i=0; i<tao->numbermonitors; ++i) {
1100: TaoSetMonitor(bnk->bncg, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);
1101: PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
1102: }
1103: }
1104: KSPGetType(tao->ksp,&ksp_type);
1105: PetscStrcmp(ksp_type,KSPNASH,&bnk->is_nash);
1106: PetscStrcmp(ksp_type,KSPSTCG,&bnk->is_stcg);
1107: PetscStrcmp(ksp_type,KSPGLTR,&bnk->is_gltr);
1108: bnk->X_inactive = NULL;
1109: bnk->G_inactive = NULL;
1110: bnk->inactive_work = NULL;
1111: bnk->active_work = NULL;
1112: bnk->inactive_idx = NULL;
1113: bnk->active_idx = NULL;
1114: bnk->active_lower = NULL;
1115: bnk->active_upper = NULL;
1116: bnk->active_fixed = NULL;
1117: bnk->M = NULL;
1118: bnk->H_inactive = NULL;
1119: bnk->Hpre_inactive = NULL;
1120: return(0);
1121: }
1123: /*------------------------------------------------------------*/
1125: PetscErrorCode TaoDestroy_BNK(Tao tao)
1126: {
1127: TAO_BNK *bnk = (TAO_BNK *)tao->data;
1131: if (tao->setupcalled) {
1132: VecDestroy(&bnk->W);
1133: VecDestroy(&bnk->Xold);
1134: VecDestroy(&bnk->Gold);
1135: VecDestroy(&bnk->Xwork);
1136: VecDestroy(&bnk->Gwork);
1137: VecDestroy(&bnk->unprojected_gradient);
1138: VecDestroy(&bnk->unprojected_gradient_old);
1139: VecDestroy(&bnk->Diag_min);
1140: VecDestroy(&bnk->Diag_max);
1141: }
1142: ISDestroy(&bnk->active_lower);
1143: ISDestroy(&bnk->active_upper);
1144: ISDestroy(&bnk->active_fixed);
1145: ISDestroy(&bnk->active_idx);
1146: ISDestroy(&bnk->inactive_idx);
1147: MatDestroy(&bnk->Hpre_inactive);
1148: MatDestroy(&bnk->H_inactive);
1149: TaoDestroy(&bnk->bncg);
1150: PetscFree(tao->data);
1151: return(0);
1152: }
1154: /*------------------------------------------------------------*/
1156: PetscErrorCode TaoSetFromOptions_BNK(PetscOptionItems *PetscOptionsObject,Tao tao)
1157: {
1158: TAO_BNK *bnk = (TAO_BNK *)tao->data;
1162: PetscOptionsHead(PetscOptionsObject,"Newton-Krylov method for bound constrained optimization");
1163: PetscOptionsEList("-tao_bnk_init_type", "radius initialization type", "", BNK_INIT, BNK_INIT_TYPES, BNK_INIT[bnk->init_type], &bnk->init_type, NULL);
1164: PetscOptionsEList("-tao_bnk_update_type", "radius update type", "", BNK_UPDATE, BNK_UPDATE_TYPES, BNK_UPDATE[bnk->update_type], &bnk->update_type, NULL);
1165: PetscOptionsEList("-tao_bnk_as_type", "active set estimation method", "", BNK_AS, BNK_AS_TYPES, BNK_AS[bnk->as_type], &bnk->as_type, NULL);
1166: PetscOptionsReal("-tao_bnk_sval", "(developer) Hessian perturbation starting value", "", bnk->sval, &bnk->sval,NULL);
1167: PetscOptionsReal("-tao_bnk_imin", "(developer) minimum initial Hessian perturbation", "", bnk->imin, &bnk->imin,NULL);
1168: PetscOptionsReal("-tao_bnk_imax", "(developer) maximum initial Hessian perturbation", "", bnk->imax, &bnk->imax,NULL);
1169: PetscOptionsReal("-tao_bnk_imfac", "(developer) initial merit factor for Hessian perturbation", "", bnk->imfac, &bnk->imfac,NULL);
1170: PetscOptionsReal("-tao_bnk_pmin", "(developer) minimum Hessian perturbation", "", bnk->pmin, &bnk->pmin,NULL);
1171: PetscOptionsReal("-tao_bnk_pmax", "(developer) maximum Hessian perturbation", "", bnk->pmax, &bnk->pmax,NULL);
1172: PetscOptionsReal("-tao_bnk_pgfac", "(developer) Hessian perturbation growth factor", "", bnk->pgfac, &bnk->pgfac,NULL);
1173: PetscOptionsReal("-tao_bnk_psfac", "(developer) Hessian perturbation shrink factor", "", bnk->psfac, &bnk->psfac,NULL);
1174: PetscOptionsReal("-tao_bnk_pmgfac", "(developer) merit growth factor for Hessian perturbation", "", bnk->pmgfac, &bnk->pmgfac,NULL);
1175: PetscOptionsReal("-tao_bnk_pmsfac", "(developer) merit shrink factor for Hessian perturbation", "", bnk->pmsfac, &bnk->pmsfac,NULL);
1176: PetscOptionsReal("-tao_bnk_eta1", "(developer) threshold for rejecting step (-tao_bnk_update_type reduction)", "", bnk->eta1, &bnk->eta1,NULL);
1177: PetscOptionsReal("-tao_bnk_eta2", "(developer) threshold for accepting marginal step (-tao_bnk_update_type reduction)", "", bnk->eta2, &bnk->eta2,NULL);
1178: PetscOptionsReal("-tao_bnk_eta3", "(developer) threshold for accepting reasonable step (-tao_bnk_update_type reduction)", "", bnk->eta3, &bnk->eta3,NULL);
1179: PetscOptionsReal("-tao_bnk_eta4", "(developer) threshold for accepting good step (-tao_bnk_update_type reduction)", "", bnk->eta4, &bnk->eta4,NULL);
1180: PetscOptionsReal("-tao_bnk_alpha1", "(developer) radius reduction factor for rejected step (-tao_bnk_update_type reduction)", "", bnk->alpha1, &bnk->alpha1,NULL);
1181: PetscOptionsReal("-tao_bnk_alpha2", "(developer) radius reduction factor for marginally accepted bad step (-tao_bnk_update_type reduction)", "", bnk->alpha2, &bnk->alpha2,NULL);
1182: PetscOptionsReal("-tao_bnk_alpha3", "(developer) radius increase factor for reasonable accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha3, &bnk->alpha3,NULL);
1183: PetscOptionsReal("-tao_bnk_alpha4", "(developer) radius increase factor for good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha4, &bnk->alpha4,NULL);
1184: PetscOptionsReal("-tao_bnk_alpha5", "(developer) radius increase factor for very good accepted step (-tao_bnk_update_type reduction)", "", bnk->alpha5, &bnk->alpha5,NULL);
1185: PetscOptionsReal("-tao_bnk_nu1", "(developer) threshold for small line-search step length (-tao_bnk_update_type step)", "", bnk->nu1, &bnk->nu1,NULL);
1186: PetscOptionsReal("-tao_bnk_nu2", "(developer) threshold for reasonable line-search step length (-tao_bnk_update_type step)", "", bnk->nu2, &bnk->nu2,NULL);
1187: PetscOptionsReal("-tao_bnk_nu3", "(developer) threshold for large line-search step length (-tao_bnk_update_type step)", "", bnk->nu3, &bnk->nu3,NULL);
1188: PetscOptionsReal("-tao_bnk_nu4", "(developer) threshold for very large line-search step length (-tao_bnk_update_type step)", "", bnk->nu4, &bnk->nu4,NULL);
1189: PetscOptionsReal("-tao_bnk_omega1", "(developer) radius reduction factor for very small line-search step length (-tao_bnk_update_type step)", "", bnk->omega1, &bnk->omega1,NULL);
1190: PetscOptionsReal("-tao_bnk_omega2", "(developer) radius reduction factor for small line-search step length (-tao_bnk_update_type step)", "", bnk->omega2, &bnk->omega2,NULL);
1191: PetscOptionsReal("-tao_bnk_omega3", "(developer) radius factor for decent line-search step length (-tao_bnk_update_type step)", "", bnk->omega3, &bnk->omega3,NULL);
1192: PetscOptionsReal("-tao_bnk_omega4", "(developer) radius increase factor for large line-search step length (-tao_bnk_update_type step)", "", bnk->omega4, &bnk->omega4,NULL);
1193: PetscOptionsReal("-tao_bnk_omega5", "(developer) radius increase factor for very large line-search step length (-tao_bnk_update_type step)", "", bnk->omega5, &bnk->omega5,NULL);
1194: PetscOptionsReal("-tao_bnk_mu1_i", "(developer) threshold for accepting very good step (-tao_bnk_init_type interpolation)", "", bnk->mu1_i, &bnk->mu1_i,NULL);
1195: PetscOptionsReal("-tao_bnk_mu2_i", "(developer) threshold for accepting good step (-tao_bnk_init_type interpolation)", "", bnk->mu2_i, &bnk->mu2_i,NULL);
1196: PetscOptionsReal("-tao_bnk_gamma1_i", "(developer) radius reduction factor for rejected very bad step (-tao_bnk_init_type interpolation)", "", bnk->gamma1_i, &bnk->gamma1_i,NULL);
1197: PetscOptionsReal("-tao_bnk_gamma2_i", "(developer) radius reduction factor for rejected bad step (-tao_bnk_init_type interpolation)", "", bnk->gamma2_i, &bnk->gamma2_i,NULL);
1198: PetscOptionsReal("-tao_bnk_gamma3_i", "(developer) radius increase factor for accepted good step (-tao_bnk_init_type interpolation)", "", bnk->gamma3_i, &bnk->gamma3_i,NULL);
1199: PetscOptionsReal("-tao_bnk_gamma4_i", "(developer) radius increase factor for accepted very good step (-tao_bnk_init_type interpolation)", "", bnk->gamma4_i, &bnk->gamma4_i,NULL);
1200: PetscOptionsReal("-tao_bnk_theta_i", "(developer) trust region interpolation factor (-tao_bnk_init_type interpolation)", "", bnk->theta_i, &bnk->theta_i,NULL);
1201: PetscOptionsReal("-tao_bnk_mu1", "(developer) threshold for accepting very good step (-tao_bnk_update_type interpolation)", "", bnk->mu1, &bnk->mu1,NULL);
1202: PetscOptionsReal("-tao_bnk_mu2", "(developer) threshold for accepting good step (-tao_bnk_update_type interpolation)", "", bnk->mu2, &bnk->mu2,NULL);
1203: PetscOptionsReal("-tao_bnk_gamma1", "(developer) radius reduction factor for rejected very bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma1, &bnk->gamma1,NULL);
1204: PetscOptionsReal("-tao_bnk_gamma2", "(developer) radius reduction factor for rejected bad step (-tao_bnk_update_type interpolation)", "", bnk->gamma2, &bnk->gamma2,NULL);
1205: PetscOptionsReal("-tao_bnk_gamma3", "(developer) radius increase factor for accepted good step (-tao_bnk_update_type interpolation)", "", bnk->gamma3, &bnk->gamma3,NULL);
1206: PetscOptionsReal("-tao_bnk_gamma4", "(developer) radius increase factor for accepted very good step (-tao_bnk_update_type interpolation)", "", bnk->gamma4, &bnk->gamma4,NULL);
1207: PetscOptionsReal("-tao_bnk_theta", "(developer) trust region interpolation factor (-tao_bnk_update_type interpolation)", "", bnk->theta, &bnk->theta,NULL);
1208: PetscOptionsReal("-tao_bnk_min_radius", "(developer) lower bound on initial radius", "", bnk->min_radius, &bnk->min_radius,NULL);
1209: PetscOptionsReal("-tao_bnk_max_radius", "(developer) upper bound on radius", "", bnk->max_radius, &bnk->max_radius,NULL);
1210: PetscOptionsReal("-tao_bnk_epsilon", "(developer) tolerance used when computing actual and predicted reduction", "", bnk->epsilon, &bnk->epsilon,NULL);
1211: PetscOptionsReal("-tao_bnk_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", bnk->as_tol, &bnk->as_tol,NULL);
1212: PetscOptionsReal("-tao_bnk_as_step", "(developer) step length used when estimating actively bounded variables", "", bnk->as_step, &bnk->as_step,NULL);
1213: PetscOptionsInt("-tao_bnk_max_cg_its", "number of BNCG iterations to take for each Newton step", "", bnk->max_cg_its, &bnk->max_cg_its,NULL);
1214: PetscOptionsTail();
1215: TaoSetFromOptions(bnk->bncg);
1216: TaoLineSearchSetFromOptions(tao->linesearch);
1217: KSPSetFromOptions(tao->ksp);
1218: return(0);
1219: }
1221: /*------------------------------------------------------------*/
1223: PetscErrorCode TaoView_BNK(Tao tao, PetscViewer viewer)
1224: {
1225: TAO_BNK *bnk = (TAO_BNK *)tao->data;
1226: PetscInt nrejects;
1227: PetscBool isascii;
1231: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1232: if (isascii) {
1233: PetscViewerASCIIPushTab(viewer);
1234: if (bnk->M) {
1235: MatLMVMGetRejectCount(bnk->M,&nrejects);
1236: PetscViewerASCIIPrintf(viewer, "Rejected BFGS updates: %D\n",nrejects);
1237: }
1238: PetscViewerASCIIPrintf(viewer, "CG steps: %D\n", bnk->tot_cg_its);
1239: PetscViewerASCIIPrintf(viewer, "Newton steps: %D\n", bnk->newt);
1240: if (bnk->M) {
1241: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", bnk->bfgs);
1242: }
1243: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", bnk->sgrad);
1244: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", bnk->grad);
1245: PetscViewerASCIIPrintf(viewer, "KSP termination reasons:\n");
1246: PetscViewerASCIIPrintf(viewer, " atol: %D\n", bnk->ksp_atol);
1247: PetscViewerASCIIPrintf(viewer, " rtol: %D\n", bnk->ksp_rtol);
1248: PetscViewerASCIIPrintf(viewer, " ctol: %D\n", bnk->ksp_ctol);
1249: PetscViewerASCIIPrintf(viewer, " negc: %D\n", bnk->ksp_negc);
1250: PetscViewerASCIIPrintf(viewer, " dtol: %D\n", bnk->ksp_dtol);
1251: PetscViewerASCIIPrintf(viewer, " iter: %D\n", bnk->ksp_iter);
1252: PetscViewerASCIIPrintf(viewer, " othr: %D\n", bnk->ksp_othr);
1253: PetscViewerASCIIPopTab(viewer);
1254: }
1255: return(0);
1256: }
1258: /* ---------------------------------------------------------- */
1260: /*MC
1261: TAOBNK - Shared base-type for Bounded Newton-Krylov type algorithms.
1262: At each iteration, the BNK methods solve the symmetric
1263: system of equations to obtain the step diretion dk:
1264: Hk dk = -gk
1265: for free variables only. The step can be globalized either through
1266: trust-region methods, or a line search, or a heuristic mixture of both.
1268: Options Database Keys:
1269: + -max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop
1270: . -init_type - trust radius initialization method ("constant", "direction", "interpolation")
1271: . -update_type - trust radius update method ("step", "direction", "interpolation")
1272: . -as_type - active-set estimation method ("none", "bertsekas")
1273: . -as_tol - (developer) initial tolerance used in estimating bounded active variables (-as_type bertsekas)
1274: . -as_step - (developer) trial step length used in estimating bounded active variables (-as_type bertsekas)
1275: . -sval - (developer) Hessian perturbation starting value
1276: . -imin - (developer) minimum initial Hessian perturbation
1277: . -imax - (developer) maximum initial Hessian perturbation
1278: . -pmin - (developer) minimum Hessian perturbation
1279: . -pmax - (developer) aximum Hessian perturbation
1280: . -pgfac - (developer) Hessian perturbation growth factor
1281: . -psfac - (developer) Hessian perturbation shrink factor
1282: . -imfac - (developer) initial merit factor for Hessian perturbation
1283: . -pmgfac - (developer) merit growth factor for Hessian perturbation
1284: . -pmsfac - (developer) merit shrink factor for Hessian perturbation
1285: . -eta1 - (developer) threshold for rejecting step (-update_type reduction)
1286: . -eta2 - (developer) threshold for accepting marginal step (-update_type reduction)
1287: . -eta3 - (developer) threshold for accepting reasonable step (-update_type reduction)
1288: . -eta4 - (developer) threshold for accepting good step (-update_type reduction)
1289: . -alpha1 - (developer) radius reduction factor for rejected step (-update_type reduction)
1290: . -alpha2 - (developer) radius reduction factor for marginally accepted bad step (-update_type reduction)
1291: . -alpha3 - (developer) radius increase factor for reasonable accepted step (-update_type reduction)
1292: . -alpha4 - (developer) radius increase factor for good accepted step (-update_type reduction)
1293: . -alpha5 - (developer) radius increase factor for very good accepted step (-update_type reduction)
1294: . -epsilon - (developer) tolerance for small pred/actual ratios that trigger automatic step acceptance (-update_type reduction)
1295: . -mu1 - (developer) threshold for accepting very good step (-update_type interpolation)
1296: . -mu2 - (developer) threshold for accepting good step (-update_type interpolation)
1297: . -gamma1 - (developer) radius reduction factor for rejected very bad step (-update_type interpolation)
1298: . -gamma2 - (developer) radius reduction factor for rejected bad step (-update_type interpolation)
1299: . -gamma3 - (developer) radius increase factor for accepted good step (-update_type interpolation)
1300: . -gamma4 - (developer) radius increase factor for accepted very good step (-update_type interpolation)
1301: . -theta - (developer) trust region interpolation factor (-update_type interpolation)
1302: . -nu1 - (developer) threshold for small line-search step length (-update_type step)
1303: . -nu2 - (developer) threshold for reasonable line-search step length (-update_type step)
1304: . -nu3 - (developer) threshold for large line-search step length (-update_type step)
1305: . -nu4 - (developer) threshold for very large line-search step length (-update_type step)
1306: . -omega1 - (developer) radius reduction factor for very small line-search step length (-update_type step)
1307: . -omega2 - (developer) radius reduction factor for small line-search step length (-update_type step)
1308: . -omega3 - (developer) radius factor for decent line-search step length (-update_type step)
1309: . -omega4 - (developer) radius increase factor for large line-search step length (-update_type step)
1310: . -omega5 - (developer) radius increase factor for very large line-search step length (-update_type step)
1311: . -mu1_i - (developer) threshold for accepting very good step (-init_type interpolation)
1312: . -mu2_i - (developer) threshold for accepting good step (-init_type interpolation)
1313: . -gamma1_i - (developer) radius reduction factor for rejected very bad step (-init_type interpolation)
1314: . -gamma2_i - (developer) radius reduction factor for rejected bad step (-init_type interpolation)
1315: . -gamma3_i - (developer) radius increase factor for accepted good step (-init_type interpolation)
1316: . -gamma4_i - (developer) radius increase factor for accepted very good step (-init_type interpolation)
1317: - -theta_i - (developer) trust region interpolation factor (-init_type interpolation)
1319: Level: beginner
1320: M*/
1322: PetscErrorCode TaoCreate_BNK(Tao tao)
1323: {
1324: TAO_BNK *bnk;
1325: const char *morethuente_type = TAOLINESEARCHMT;
1327: PC pc;
1330: PetscNewLog(tao,&bnk);
1332: tao->ops->setup = TaoSetUp_BNK;
1333: tao->ops->view = TaoView_BNK;
1334: tao->ops->setfromoptions = TaoSetFromOptions_BNK;
1335: tao->ops->destroy = TaoDestroy_BNK;
1337: /* Override default settings (unless already changed) */
1338: if (!tao->max_it_changed) tao->max_it = 50;
1339: if (!tao->trust0_changed) tao->trust0 = 100.0;
1341: tao->data = (void*)bnk;
1343: /* Hessian shifting parameters */
1344: bnk->computehessian = TaoBNKComputeHessian;
1345: bnk->computestep = TaoBNKComputeStep;
1347: bnk->sval = 0.0;
1348: bnk->imin = 1.0e-4;
1349: bnk->imax = 1.0e+2;
1350: bnk->imfac = 1.0e-1;
1352: bnk->pmin = 1.0e-12;
1353: bnk->pmax = 1.0e+2;
1354: bnk->pgfac = 1.0e+1;
1355: bnk->psfac = 4.0e-1;
1356: bnk->pmgfac = 1.0e-1;
1357: bnk->pmsfac = 1.0e-1;
1359: /* Default values for trust-region radius update based on steplength */
1360: bnk->nu1 = 0.25;
1361: bnk->nu2 = 0.50;
1362: bnk->nu3 = 1.00;
1363: bnk->nu4 = 1.25;
1365: bnk->omega1 = 0.25;
1366: bnk->omega2 = 0.50;
1367: bnk->omega3 = 1.00;
1368: bnk->omega4 = 2.00;
1369: bnk->omega5 = 4.00;
1371: /* Default values for trust-region radius update based on reduction */
1372: bnk->eta1 = 1.0e-4;
1373: bnk->eta2 = 0.25;
1374: bnk->eta3 = 0.50;
1375: bnk->eta4 = 0.90;
1377: bnk->alpha1 = 0.25;
1378: bnk->alpha2 = 0.50;
1379: bnk->alpha3 = 1.00;
1380: bnk->alpha4 = 2.00;
1381: bnk->alpha5 = 4.00;
1383: /* Default values for trust-region radius update based on interpolation */
1384: bnk->mu1 = 0.10;
1385: bnk->mu2 = 0.50;
1387: bnk->gamma1 = 0.25;
1388: bnk->gamma2 = 0.50;
1389: bnk->gamma3 = 2.00;
1390: bnk->gamma4 = 4.00;
1392: bnk->theta = 0.05;
1394: /* Default values for trust region initialization based on interpolation */
1395: bnk->mu1_i = 0.35;
1396: bnk->mu2_i = 0.50;
1398: bnk->gamma1_i = 0.0625;
1399: bnk->gamma2_i = 0.5;
1400: bnk->gamma3_i = 2.0;
1401: bnk->gamma4_i = 5.0;
1403: bnk->theta_i = 0.25;
1405: /* Remaining parameters */
1406: bnk->max_cg_its = 0;
1407: bnk->min_radius = 1.0e-10;
1408: bnk->max_radius = 1.0e10;
1409: bnk->epsilon = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
1410: bnk->as_tol = 1.0e-3;
1411: bnk->as_step = 1.0e-3;
1412: bnk->dmin = 1.0e-6;
1413: bnk->dmax = 1.0e6;
1415: bnk->M = NULL;
1416: bnk->bfgs_pre = NULL;
1417: bnk->init_type = BNK_INIT_INTERPOLATION;
1418: bnk->update_type = BNK_UPDATE_REDUCTION;
1419: bnk->as_type = BNK_AS_BERTSEKAS;
1421: bnk->is_stcg = PETSC_FALSE;
1422: bnk->is_gltr = PETSC_FALSE;
1423: bnk->is_nash = PETSC_FALSE;
1425: /* Create the embedded BNCG solver */
1426: TaoCreate(PetscObjectComm((PetscObject)tao), &bnk->bncg);
1427: PetscObjectIncrementTabLevel((PetscObject)bnk->bncg, (PetscObject)tao, 1);
1428: TaoSetOptionsPrefix(bnk->bncg, "tao_bnk_");
1429: TaoSetType(bnk->bncg, TAOBNCG);
1431: /* Create the line search */
1432: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
1433: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
1434: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
1435: TaoLineSearchSetType(tao->linesearch,morethuente_type);
1436: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
1438: /* Set linear solver to default for symmetric matrices */
1439: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
1440: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
1441: KSPSetOptionsPrefix(tao->ksp,"tao_bnk_");
1442: KSPSetType(tao->ksp,KSPSTCG);
1443: KSPGetPC(tao->ksp, &pc);
1444: PCSetType(pc, PCLMVM);
1445: return(0);
1446: }