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