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