Actual source code: ssfls.c
1: #include "ssls.h"
5: PetscErrorCode TaoSetUp_SSFLS(TaoSolver tao)
6: {
7: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
12: VecDuplicate(tao->solution,&tao->gradient);
13: VecDuplicate(tao->solution,&tao->stepdirection);
14: VecDuplicate(tao->solution,&ssls->w);
15: VecDuplicate(tao->solution,&ssls->ff);
16: VecDuplicate(tao->solution,&ssls->dpsi);
17: VecDuplicate(tao->solution,&ssls->da);
18: VecDuplicate(tao->solution,&ssls->db);
19: VecDuplicate(tao->solution,&ssls->t1);
20: VecDuplicate(tao->solution,&ssls->t2);
21: if (!tao->XL) {
22: VecDuplicate(tao->solution,&tao->XL);
23: }
24: if (!tao->XU) {
25: VecDuplicate(tao->solution,&tao->XU);
26: }
27: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
28: return(0);
29: }
34: static PetscErrorCode TaoSolve_SSFLS(TaoSolver tao)
35: {
36: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
37: PetscReal psi, ndpsi, normd, innerd, t=0;
38: PetscReal delta, rho;
39: PetscInt iter=0,kspits;
40: TaoSolverTerminationReason reason;
41: TaoLineSearchTerminationReason ls_reason;
46: /* Assume that Setup has been called!
47: Set the structure for the Jacobian and create a linear solver. */
48:
49: delta = ssls->delta;
50: rho = ssls->rho;
52: TaoComputeVariableBounds(tao);
53: /* Project solution inside bounds */
54: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
55: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
56:
57: TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);
59: /* Calculate the function value and fischer function value at the
60: current iterate */
61: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
62: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
65: while (1) {
66: ierr=PetscInfo3(tao, "iter: %D, merit: %G, ndpsi: %G\n",
67: iter, ssls->merit, ndpsi);
68: /* Check the termination criteria */
69: TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);
70:
71: if (reason!=TAO_CONTINUE_ITERATING) break;
73: /* Calculate direction. (Really negative of newton direction. Therefore,
74: rest of the code uses -d.) */
75: KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre,ssls->matflag);
76: KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
77: KSPGetIterationNumber(tao->ksp,&kspits);
78: tao->ksp_its+=kspits;
79:
80: VecCopy(tao->stepdirection,ssls->w);
81: VecScale(ssls->w,-1.0);
82: VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w);
83:
84: VecNorm(ssls->w,NORM_2,&normd);
85: VecDot(ssls->w,ssls->dpsi,&innerd);
87: /* Make sure that we have a descent direction */
88: if (innerd >= -delta*pow(normd, rho)) {
89: PetscInfo(tao, "newton direction not descent\n");
90: VecCopy(ssls->dpsi,tao->stepdirection);
91: VecDot(ssls->w,ssls->dpsi,&innerd);
92: }
94: VecScale(tao->stepdirection, -1.0);
95: innerd = -innerd;
97: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
98: TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
99: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
100: }
102: return(0);
103: }
107: PetscErrorCode TaoDestroy_SSFLS(TaoSolver tao)
108: {
109: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
114: VecDestroy(&ssls->ff);
115: VecDestroy(&ssls->w);
116: VecDestroy(&ssls->dpsi);
117: VecDestroy(&ssls->da);
118: VecDestroy(&ssls->db);
119: VecDestroy(&ssls->t1);
120: VecDestroy(&ssls->t2);
121: PetscFree(tao->data);
122: tao->data = PETSC_NULL;
123: return(0);
124: }
126: /* ---------------------------------------------------------- */
130: PetscErrorCode TaoCreate_SSFLS(TaoSolver tao)
131: {
132: TAO_SSLS *ssls;
134: const char *armijo_type = TAOLINESEARCH_ARMIJO;
137: PetscNewLog(tao,TAO_SSLS,&ssls);
138: tao->data = (void*)ssls;
139: tao->ops->solve=TaoSolve_SSFLS;
140: tao->ops->setup=TaoSetUp_SSFLS;
141: tao->ops->view=TaoView_SSLS;
142: tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
143: tao->ops->destroy = TaoDestroy_SSFLS;
145: ssls->delta = 1e-10;
146: ssls->rho = 2.1;
148: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
149: TaoLineSearchSetType(tao->linesearch,armijo_type);
150: TaoLineSearchSetFromOptions(tao->linesearch);
151: /* Linesearch objective and objectivegradient routines are
152: set in solve routine */
153: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
154:
155: tao->max_it = 2000;
156: tao->max_funcs = 4000;
157: tao->fatol = 0; tao->frtol = 0; tao->grtol=0; tao->grtol=0;
158: tao->gatol = 1.0e-16;
159: tao->fmin = 1.0e-8;
161: return(0);
162: }