1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2
TaoSetUp_SSILS(Tao tao)3 static PetscErrorCode TaoSetUp_SSILS(Tao tao)
4 {
5 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
6
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 }
18
TaoDestroy_SSILS(Tao tao)19 static PetscErrorCode TaoDestroy_SSILS(Tao tao)
20 {
21 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
22
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 }
34
TaoSolve_SSILS(Tao tao)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;
41
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;
47
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));
52
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));
57
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;
66
67 /* Call general purpose update function */
68 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
69 tao->niter++;
70
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));
79
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 }
86
87 PetscCall(VecScale(tao->stepdirection, -1.0));
88 innerd = -innerd;
89
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 }
96
97 /*MC
98 TAOSSILS - semi-smooth infeasible linesearch algorithm for solving
99 complementarity constraints
100
101 Options Database Keys:
102 + -tao_ssls_delta - descent test fraction
103 - -tao_ssls_rho - descent test power
104
105 Level: beginner
106 M*/
TaoCreate_SSILS(Tao tao)107 PETSC_EXTERN PetscErrorCode TaoCreate_SSILS(Tao tao)
108 {
109 TAO_SSLS *ssls;
110 const char *armijo_type = TAOLINESEARCHARMIJO;
111
112 PetscFunctionBegin;
113 PetscCall(PetscNew(&ssls));
114 tao->data = (void *)ssls;
115 tao->ops->solve = TaoSolve_SSILS;
116 tao->ops->setup = TaoSetUp_SSILS;
117 tao->ops->view = TaoView_SSLS;
118 tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
119 tao->ops->destroy = TaoDestroy_SSILS;
120
121 ssls->delta = 1e-10;
122 ssls->rho = 2.1;
123
124 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
125 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
126 PetscCall(TaoLineSearchSetType(tao->linesearch, armijo_type));
127 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
128 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
129 /* Note: linesearch objective and objectivegradient routines are set in solve routine */
130 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
131 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
132 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
133
134 /* Override default settings (unless already changed) */
135 PetscCall(TaoParametersInitialize(tao));
136 PetscObjectParameterSetDefault(tao, max_it, 2000);
137 PetscObjectParameterSetDefault(tao, max_funcs, 4000);
138 PetscObjectParameterSetDefault(tao, gttol, 0);
139 PetscObjectParameterSetDefault(tao, grtol, 0);
140 PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1.0e-6 : 1.0e-16);
141 PetscObjectParameterSetDefault(tao, fmin, PetscDefined(USE_REAL_SINGLE) ? 1.0e-4 : 1.0e-8);
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144