xref: /petsc/src/tao/complementarity/impls/ssls/ssfls.c (revision 1ebf93c6b7d760d592de6ebe6cdc0debaa3caf75)
1 #include <../src/tao/complementarity/impls/ssls/ssls.h>
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "TaoSetUp_SSFLS"
5 PetscErrorCode TaoSetUp_SSFLS(Tao tao)
6 {
7   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;
8   PetscErrorCode ierr;
9 
10   PetscFunctionBegin;
11   ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);
12   ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);
13   ierr = VecDuplicate(tao->solution,&ssls->w);CHKERRQ(ierr);
14   ierr = VecDuplicate(tao->solution,&ssls->ff);CHKERRQ(ierr);
15   ierr = VecDuplicate(tao->solution,&ssls->dpsi);CHKERRQ(ierr);
16   ierr = VecDuplicate(tao->solution,&ssls->da);CHKERRQ(ierr);
17   ierr = VecDuplicate(tao->solution,&ssls->db);CHKERRQ(ierr);
18   ierr = VecDuplicate(tao->solution,&ssls->t1);CHKERRQ(ierr);
19   ierr = VecDuplicate(tao->solution,&ssls->t2);CHKERRQ(ierr);
20   if (!tao->XL) {
21     ierr = VecDuplicate(tao->solution,&tao->XL);CHKERRQ(ierr);
22   }
23   if (!tao->XU) {
24     ierr = VecDuplicate(tao->solution,&tao->XU);CHKERRQ(ierr);
25   }
26   ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr);
27   PetscFunctionReturn(0);
28 }
29 
30 #undef __FUNCT__
31 #define __FUNCT__ "TaoSolve_SSFLS"
32 static PetscErrorCode TaoSolve_SSFLS(Tao tao)
33 {
34   TAO_SSLS                     *ssls = (TAO_SSLS *)tao->data;
35   PetscReal                    psi, ndpsi, normd, innerd, t=0;
36   PetscReal                    delta, rho;
37   TaoConvergedReason           reason;
38   TaoLineSearchConvergedReason ls_reason;
39   PetscErrorCode               ierr;
40 
41   PetscFunctionBegin;
42   /* Assume that Setup has been called!
43      Set the structure for the Jacobian and create a linear solver. */
44   delta = ssls->delta;
45   rho = ssls->rho;
46 
47   ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
48   /* Project solution inside bounds */
49   ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr);
50   ierr = TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);CHKERRQ(ierr);
51   ierr = TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao);CHKERRQ(ierr);
52 
53   /* Calculate the function value and fischer function value at the
54      current iterate */
55   ierr = TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi);CHKERRQ(ierr);
56   ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
57 
58   while (PETSC_TRUE) {
59     ierr=PetscInfo3(tao, "iter: %D, merit: %g, ndpsi: %g\n",tao->niter, (double)ssls->merit, (double)ndpsi);CHKERRQ(ierr);
60     /* Check the termination criteria */
61     ierr = TaoMonitor(tao,tao->niter,ssls->merit,ndpsi,0.0,t,&reason);CHKERRQ(ierr);
62     if (reason!=TAO_CONTINUE_ITERATING) break;
63     tao->niter++;
64 
65     /* Calculate direction.  (Really negative of newton direction.  Therefore,
66        rest of the code uses -d.) */
67     ierr = KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre);CHKERRQ(ierr);
68     ierr = KSPSolve(tao->ksp,ssls->ff,tao->stepdirection);CHKERRQ(ierr);
69     ierr = KSPGetIterationNumber(tao->ksp,&tao->ksp_its);CHKERRQ(ierr);
70     tao->ksp_tot_its+=tao->ksp_its;
71 
72     ierr = VecCopy(tao->stepdirection,ssls->w);CHKERRQ(ierr);
73     ierr = VecScale(ssls->w,-1.0);CHKERRQ(ierr);
74     ierr = VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w);CHKERRQ(ierr);
75 
76     ierr = VecNorm(ssls->w,NORM_2,&normd);CHKERRQ(ierr);
77     ierr = VecDot(ssls->w,ssls->dpsi,&innerd);CHKERRQ(ierr);
78 
79     /* Make sure that we have a descent direction */
80     if (innerd >= -delta*PetscPowReal(normd, rho)) {
81       ierr = PetscInfo(tao, "newton direction not descent\n");CHKERRQ(ierr);
82       ierr = VecCopy(ssls->dpsi,tao->stepdirection);CHKERRQ(ierr);
83       ierr = VecDot(ssls->w,ssls->dpsi,&innerd);CHKERRQ(ierr);
84     }
85 
86     ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
87     innerd = -innerd;
88 
89     ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);CHKERRQ(ierr);
90     ierr = TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason);CHKERRQ(ierr);
91     ierr = VecNorm(ssls->dpsi,NORM_2,&ndpsi);CHKERRQ(ierr);
92   }
93   PetscFunctionReturn(0);
94 }
95 
96 #undef __FUNCT__
97 #define __FUNCT__ "TaoDestroy_SSFLS"
98 PetscErrorCode TaoDestroy_SSFLS(Tao tao)
99 {
100   TAO_SSLS       *ssls = (TAO_SSLS *)tao->data;
101   PetscErrorCode ierr;
102 
103   PetscFunctionBegin;
104   ierr = VecDestroy(&ssls->ff);CHKERRQ(ierr);
105   ierr = VecDestroy(&ssls->w);CHKERRQ(ierr);
106   ierr = VecDestroy(&ssls->dpsi);CHKERRQ(ierr);
107   ierr = VecDestroy(&ssls->da);CHKERRQ(ierr);
108   ierr = VecDestroy(&ssls->db);CHKERRQ(ierr);
109   ierr = VecDestroy(&ssls->t1);CHKERRQ(ierr);
110   ierr = VecDestroy(&ssls->t2);CHKERRQ(ierr);
111   ierr = PetscFree(tao->data);CHKERRQ(ierr);
112   PetscFunctionReturn(0);
113 }
114 
115 /* ---------------------------------------------------------- */
116 /*MC
117    TAOSSFLS - Semi-smooth feasible linesearch algorithm for solving
118        complementarity constraints
119 
120    Options Database Keys:
121 + -tao_ssls_delta - descent test fraction
122 - -tao_ssls_rho - descent test power
123 
124    Level: beginner
125 M*/
126 
127 #undef __FUNCT__
128 #define __FUNCT__ "TaoCreate_SSFLS"
129 PETSC_EXTERN PetscErrorCode TaoCreate_SSFLS(Tao tao)
130 {
131   TAO_SSLS       *ssls;
132   PetscErrorCode ierr;
133   const char     *armijo_type = TAOLINESEARCHARMIJO;
134 
135   PetscFunctionBegin;
136   ierr = PetscNewLog(tao,&ssls);CHKERRQ(ierr);
137   tao->data = (void*)ssls;
138   tao->ops->solve=TaoSolve_SSFLS;
139   tao->ops->setup=TaoSetUp_SSFLS;
140   tao->ops->view=TaoView_SSLS;
141   tao->ops->setfromoptions = TaoSetFromOptions_SSLS;
142   tao->ops->destroy = TaoDestroy_SSFLS;
143 
144   ssls->delta = 1e-10;
145   ssls->rho = 2.1;
146 
147   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr);
148   ierr = TaoLineSearchSetType(tao->linesearch,armijo_type);CHKERRQ(ierr);
149   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
150   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
151   /* Linesearch objective and objectivegradient routines are  set in solve routine */
152   ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr);
153   ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
154 
155   /* Override default settings (unless already changed) */
156   if (!tao->max_it_changed) tao->max_it = 2000;
157   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
158   if (!tao->gttol_changed) tao->gttol = 0;
159   if (!tao->grtol_changed) tao->grtol = 0;
160 #if defined(PETSC_USE_REAL_SINGLE)
161   if (!tao->gatol_changed) tao->gatol = 1.0e-6;
162   if (!tao->fmin_changed)  tao->fmin = 1.0e-4;
163 #else
164   if (!tao->gatol_changed) tao->gatol = 1.0e-16;
165   if (!tao->fmin_changed)  tao->fmin = 1.0e-8;
166 #endif
167   PetscFunctionReturn(0);
168 }
169 
170 
171