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