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