1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2
TaoSetUp_SSFLS(Tao tao)3 static PetscErrorCode TaoSetUp_SSFLS(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->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 }
20
TaoSolve_SSFLS(Tao tao)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;
27
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;
33
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));
39
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));
44
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;
53
54 /* Call general purpose update function */
55 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
56 tao->niter++;
57
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;
64
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));
68
69 PetscCall(VecNorm(ssls->w, NORM_2, &normd));
70 PetscCall(VecDot(ssls->w, ssls->dpsi, &innerd));
71
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 }
78
79 PetscCall(VecScale(tao->stepdirection, -1.0));
80 innerd = -innerd;
81
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 }
88
TaoDestroy_SSFLS(Tao tao)89 static PetscErrorCode TaoDestroy_SSFLS(Tao tao)
90 {
91 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
92
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 }
105
106 /*MC
107 TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
108 complementarity constraints
109
110 Options Database Keys:
111 + -tao_ssls_delta - descent test fraction
112 - -tao_ssls_rho - descent test power
113
114 Level: beginner
115 M*/
116
TaoCreate_SSFLS(Tao tao)117 PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
118 {
119 TAO_SSLS *ssls;
120 const char *armijo_type = TAOLINESEARCHARMIJO;
121
122 PetscFunctionBegin;
123 PetscCall(PetscNew(&ssls));
124 tao->data = (void *)ssls;
125 tao->ops->solve = TaoSolve_SSFLS;
126 tao->ops->setup = TaoSetUp_SSFLS;
127 tao->ops->view = TaoView_SSLS;
128 tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
129 tao->ops->destroy = TaoDestroy_SSFLS;
130
131 ssls->delta = 1e-10;
132 ssls->rho = 2.1;
133
134 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
135 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
136 PetscCall(TaoLineSearchSetType(tao->linesearch, armijo_type));
137 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
138 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
139 /* Linesearch objective and objectivegradient routines are set in solve routine */
140 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
141 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
142 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
143
144 /* Override default settings (unless already changed) */
145 PetscCall(TaoParametersInitialize(tao));
146 PetscObjectParameterSetDefault(tao, max_it, 2000);
147 PetscObjectParameterSetDefault(tao, max_funcs, 4000);
148 PetscObjectParameterSetDefault(tao, gttol, 0);
149 PetscObjectParameterSetDefault(tao, grtol, 0);
150 PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1.0e-6 : 1.0e-16);
151 PetscObjectParameterSetDefault(tao, fmin, PetscDefined(USE_REAL_SINGLE) ? 1.0e-4 : 1.0e-8);
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154