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