1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2
TaoSetUp_SSILS(Tao tao)3 PetscErrorCode TaoSetUp_SSILS(Tao tao)
4 {
5 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
6 PetscErrorCode ierr;
7
8 PetscFunctionBegin;
9 ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);
10 ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);
11 ierr = VecDuplicate(tao->solution,&ssls->ff);CHKERRQ(ierr);
12 ierr = VecDuplicate(tao->solution,&ssls->dpsi);CHKERRQ(ierr);
13 ierr = VecDuplicate(tao->solution,&ssls->da);CHKERRQ(ierr);
14 ierr = VecDuplicate(tao->solution,&ssls->db);CHKERRQ(ierr);
15 ierr = VecDuplicate(tao->solution,&ssls->t1);CHKERRQ(ierr);
16 ierr = VecDuplicate(tao->solution,&ssls->t2);CHKERRQ(ierr);
17 PetscFunctionReturn(0);
18 }
19
TaoDestroy_SSILS(Tao tao)20 PetscErrorCode TaoDestroy_SSILS(Tao tao)
21 {
22 TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
23 PetscErrorCode ierr;
24
25 PetscFunctionBegin;
26 ierr = VecDestroy(&ssls->ff);CHKERRQ(ierr);
27 ierr = VecDestroy(&ssls->dpsi);CHKERRQ(ierr);
28 ierr = VecDestroy(&ssls->da);CHKERRQ(ierr);
29 ierr = VecDestroy(&ssls->db);CHKERRQ(ierr);
30 ierr = VecDestroy(&ssls->t1);CHKERRQ(ierr);
31 ierr = VecDestroy(&ssls->t2);CHKERRQ(ierr);
32 ierr = PetscFree(tao->data);CHKERRQ(ierr);
33 PetscFunctionReturn(0);
34 }
35
TaoSolve_SSILS(Tao tao)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;
43
44 PetscFunctionBegin;
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;
49
50 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
51 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr);
52 ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);CHKERRQ(ierr);
53 ierr = TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);CHKERRQ(ierr);
54
55 /* Calculate the function value and fischer function value at the
56 current iterate */
57 ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);CHKERRQ(ierr);
58 ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
59
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);CHKERRQ(ierr);
63 /* Check the termination criteria */
64 ierr = TaoLogConvergenceHistory(tao,ssls->merit,ndpsi,0.0,tao->ksp_its);CHKERRQ(ierr);
65 ierr = TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t);CHKERRQ(ierr);
66 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
67 if (tao->reason!=TAO_CONTINUE_ITERATING) break;
68
69 /* Call general purpose update function */
70 if (tao->ops->update) {
71 ierr = (*tao->ops->update)(tao, tao->niter, tao->user_update);CHKERRQ(ierr);
72 }
73 tao->niter++;
74
75 /* Calculate direction. (Really negative of newton direction. Therefore,
76 rest of the code uses -d.) */
77 ierr = KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);CHKERRQ(ierr);
78 ierr = KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);CHKERRQ(ierr);
79 ierr = KSPGetIterationNumber(tao->ksp,&tao->ksp_its);CHKERRQ(ierr);
80 tao->ksp_tot_its+=tao->ksp_its;
81 ierr = VecNorm(tao->stepdirection,NORM_2,&normd);CHKERRQ(ierr);
82 ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr);
83
84 /* Make sure that we have a descent direction */
85 if (innerd <= delta*PetscPowReal(normd, rho)) {
86 ierr = PetscInfo(tao, "newton direction not descent\n");CHKERRQ(ierr);
87 ierr = VecCopy(ssls->dpsi,tao->stepdirection);CHKERRQ(ierr);
88 ierr = VecDot(tao->stepdirection,ssls->dpsi,&innerd);CHKERRQ(ierr);
89 }
90
91 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
92 innerd = -innerd;
93
94 ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr);
95 ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);CHKERRQ(ierr);
96 ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
97 }
98 PetscFunctionReturn(0);
99 }
100
101 /* ---------------------------------------------------------- */
102 /*MC
103 TAOSSILS - semi-smooth infeasible linesearch algorithm for solving
104 complementarity constraints
105
106 Options Database Keys:
107 + -tao_ssls_delta - descent test fraction
108 - -tao_ssls_rho - descent test power
109
110 Level: beginner
111 M*/
TaoCreate_SSILS(Tao tao)112 PETSC_EXTERN PetscErrorCode TaoCreate_SSILS(Tao tao)
113 {
114 TAO_SSLS *ssls;
115 PetscErrorCode ierr;
116 const char *armijo_type = TAOLINESEARCHARMIJO;
117
118 PetscFunctionBegin;
119 ierr = PetscNewLog(tao,&ssls);CHKERRQ(ierr);
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;
126
127 ssls->delta = 1e-10;
128 ssls->rho = 2.1;
129
130 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr);
131 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr);
132 ierr = TaoLineSearchSetType(tao->linesearch,armijo_type);CHKERRQ(ierr);
133 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
134 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
135 /* Note: linesearch objective and objectivegradient routines are set in solve routine */
136 ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr);
137 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);CHKERRQ(ierr);
138 ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
139
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 PetscFunctionReturn(0);
153 }
154