Actual source code: ssfls.c
1: #include <../src/tao/complementarity/impls/ssls/ssls.h>
3: PetscErrorCode TaoSetUp_SSFLS(Tao tao)
4: {
5: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
9: VecDuplicate(tao->solution,&tao->gradient);
10: VecDuplicate(tao->solution,&tao->stepdirection);
11: VecDuplicate(tao->solution,&ssls->w);
12: VecDuplicate(tao->solution,&ssls->ff);
13: VecDuplicate(tao->solution,&ssls->dpsi);
14: VecDuplicate(tao->solution,&ssls->da);
15: VecDuplicate(tao->solution,&ssls->db);
16: VecDuplicate(tao->solution,&ssls->t1);
17: VecDuplicate(tao->solution,&ssls->t2);
18: if (!tao->XL) {
19: VecDuplicate(tao->solution,&tao->XL);
20: }
21: if (!tao->XU) {
22: VecDuplicate(tao->solution,&tao->XU);
23: }
24: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
25: return(0);
26: }
28: static PetscErrorCode TaoSolve_SSFLS(Tao tao)
29: {
30: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
31: PetscReal psi, ndpsi, normd, innerd, t=0;
32: PetscReal delta, rho;
33: TaoLineSearchConvergedReason ls_reason;
34: PetscErrorCode ierr;
37: /* Assume that Setup has been called!
38: Set the structure for the Jacobian and create a linear solver. */
39: delta = ssls->delta;
40: rho = ssls->rho;
42: TaoComputeVariableBounds(tao);
43: /* Project solution inside bounds */
44: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
45: TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
46: TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);
48: /* Calculate the function value and fischer function value at the
49: current iterate */
50: TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);
51: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
53: tao->reason = TAO_CONTINUE_ITERATING;
54: while (PETSC_TRUE) {
55: ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",tao->niter, (double)ssls->merit, (double)ndpsi);
56: /* Check the termination criteria */
57: TaoLogConvergenceHistory(tao,ssls->merit,ndpsi,0.0,tao->ksp_its);
58: TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t);
59: (*tao->ops->convergencetest)(tao,tao->cnvP);
60: if (tao->reason!=TAO_CONTINUE_ITERATING) break;
62: /* Call general purpose update function */
63: if (tao->ops->update) {
64: (*tao->ops->update)(tao, tao->niter, tao->user_update);
65: }
66: tao->niter++;
68: /* Calculate direction. (Really negative of newton direction. Therefore,
69: rest of the code uses -d.) */
70: KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);
71: KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);
72: KSPGetIterationNumber(tao->ksp,&tao->ksp_its);
73: tao->ksp_tot_its+=tao->ksp_its;
75: VecCopy(tao->stepdirection,ssls->w);
76: VecScale(ssls->w,-1.0);
77: VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w);
79: VecNorm(ssls->w,NORM_2,&normd);
80: VecDot(ssls->w,ssls->dpsi,&innerd);
82: /* Make sure that we have a descent direction */
83: if (innerd >= -delta*PetscPowReal(normd, rho)) {
84: PetscInfo(tao, "newton direction not descent\n");
85: VecCopy(ssls->dpsi,tao->stepdirection);
86: VecDot(ssls->w,ssls->dpsi,&innerd);
87: }
89: VecScale(tao->stepdirection, -1.0);
90: innerd = -innerd;
92: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
93: TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);
94: VecNorm(ssls->dpsi,NORM_2,&ndpsi);
95: }
96: return(0);
97: }
99: PetscErrorCode TaoDestroy_SSFLS(Tao tao)
100: {
101: TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
105: VecDestroy(&ssls->ff);
106: VecDestroy(&ssls->w);
107: VecDestroy(&ssls->dpsi);
108: VecDestroy(&ssls->da);
109: VecDestroy(&ssls->db);
110: VecDestroy(&ssls->t1);
111: VecDestroy(&ssls->t2);
112: PetscFree(tao->data);
113: return(0);
114: }
116: /* ---------------------------------------------------------- */
117: /*MC
118: TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
119: complementarity constraints
121: Options Database Keys:
122: + -tao_ssls_delta - descent test fraction
123: - -tao_ssls_rho - descent test power
125: Level: beginner
126: M*/
128: PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
129: {
130: TAO_SSLS *ssls;
132: const char *armijo_type = TAOLINESEARCHARMIJO;
135: PetscNewLog(tao,&ssls);
136: tao->data = (void*)ssls;
137: tao->ops->solve=TaoSolve_SSFLS;
138: tao->ops->setup=TaoSetUp_SSFLS;
139: tao->ops->view=TaoView_SSLS;
140: tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
141: tao->ops->destroy = TaoDestroy_SSFLS;
143: ssls->delta = 1e-10;
144: ssls->rho = 2.1;
146: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
147: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
148: TaoLineSearchSetType(tao->linesearch,armijo_type);
149: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
150: TaoLineSearchSetFromOptions(tao->linesearch);
151: /* Linesearch objective and objectivegradient routines are set in solve routine */
152: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
153: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
154: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
156: /* Override default settings (unless already changed) */
157: if (!tao->max_it_changed) tao->max_it = 2000;
158: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
159: if (!tao->gttol_changed) tao->gttol = 0;
160: if (!tao->grtol_changed) tao->grtol = 0;
161: #if defined(PETSC_USE_REAL_SINGLE)
162: if (!tao->gatol_changed) tao->gatol = 1.0e-6;
163: if (!tao->fmin_changed) tao->fmin = 1.0e-4;
164: #else
165: if (!tao->gatol_changed) tao->gatol = 1.0e-16;
166: if (!tao->fmin_changed) tao->fmin = 1.0e-8;
167: #endif
168: return(0);
169: }