Actual source code: ssils.c
1: #include "ssls.h"
5: PetscErrorCode TaoSetUp_SSILS(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->ff);
15: VecDuplicate(tao->solution,&ssls->dpsi);
16: VecDuplicate(tao->solution,&ssls->da);
17: VecDuplicate(tao->solution,&ssls->db);
18: VecDuplicate(tao->solution,&ssls->t1);
19: VecDuplicate(tao->solution,&ssls->t2);
20: return(0);
21: }
26: PetscErrorCode TaoDestroy_SSILS(TaoSolver tao)
27: {
28: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
33: VecDestroy(&ssls->ff);
34: VecDestroy(&ssls->dpsi);
35: VecDestroy(&ssls->da);
36: VecDestroy(&ssls->db);
37: VecDestroy(&ssls->t1);
38: VecDestroy(&ssls->t2);
39: PetscFree(tao->data);
40: tao->data = PETSC_NULL;
41: return(0);
42: }
46: static PetscErrorCode TaoSolve_SSILS(TaoSolver tao)
47: {
48: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
49: PetscReal psi, ndpsi, normd, innerd, t=0;
50: PetscReal delta, rho;
51: PetscInt iter=0,kspits;
52: TaoSolverTerminationReason reason;
53: TaoLineSearchTerminationReason ls_reason;
58: /* Assume that Setup has been called!
59: Set the structure for the Jacobian and create a linear solver. */
60:
61: delta = ssls->delta;
62: rho = ssls->rho;
64: TaoComputeVariableBounds(tao);
65: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
66: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
67:
68: TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);
70: /* Calculate the function value and fischer function value at the
71: current iterate */
72: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
73: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
76: while (1) {
77: ierr=PetscInfo3(tao, "iter: %D, merit: %G, ndpsi: %G\n",
78: iter, ssls->merit, ndpsi);
79: /* Check the termination criteria */
80: TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason);
81:
82: if (reason!=TAO_CONTINUE_ITERATING) break;
84: /* Calculate direction. (Really negative of newton direction. Therefore,
85: rest of the code uses -d.) */
86: KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre,ssls->matflag);
87: KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
88: KSPGetIterationNumber(tao->ksp,&kspits);
89: tao->ksp_its+=kspits;
90: VecNorm(tao->stepdirection,NORM_2,&normd);
91: VecDot(tao->stepdirection,ssls->dpsi,&innerd);
93: /* Make sure that we have a descent direction */
94: if (innerd <= delta*pow(normd, rho)) {
95: PetscInfo(tao, "newton direction not descent\n");
96: VecCopy(ssls->dpsi,tao->stepdirection);
97: VecDot(tao->stepdirection,ssls->dpsi,&innerd);
98: }
100: VecScale(tao->stepdirection, -1.0);
101: innerd = -innerd;
103: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
104: TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
105: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
106: }
108: return(0);
109: }
111: /* ---------------------------------------------------------- */
115: PetscErrorCode TaoCreate_SSILS(TaoSolver tao)
116: {
117: TAO_SSLS *ssls;
119: const char *armijo_type = TAOLINESEARCH_ARMIJO;
122: PetscNewLog(tao,TAO_SSLS,&ssls);
123: tao->data = (void*)ssls;
124: tao->ops->solve=TaoSolve_SSILS;
125: tao->ops->setup=TaoSetUp_SSILS;
126: tao->ops->view=TaoView_SSLS;
127: tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
128: tao->ops->destroy = TaoDestroy_SSILS;
130: ssls->delta = 1e-10;
131: ssls->rho = 2.1;
133: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
134: TaoLineSearchSetType(tao->linesearch,armijo_type);
135: TaoLineSearchSetFromOptions(tao->linesearch);
136: /* Note: linesearch objective and objectivegradient routines are set in solve routine */
137: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
138:
139: tao->max_it = 2000;
140: tao->max_funcs = 4000;
141: tao->fatol = 0; tao->frtol = 0; tao->gttol=0; tao->grtol=0;
142: tao->gatol = 1.0e-16;
143: tao->fmin = 1.0e-8;
145: return(0);
146: }