1 #include <../src/tao/complementarity/impls/ssls/ssls.h> 2 /* 3 Context for ASXLS 4 -- active-set - reduced matrices formed 5 - inherit properties of original system 6 -- semismooth (S) - function not differentiable 7 - merit function continuously differentiable 8 - Fischer-Burmeister reformulation of complementarity 9 - Billups composition for two finite bounds 10 -- infeasible (I) - iterates not guaranteed to remain within bounds 11 -- feasible (F) - iterates guaranteed to remain within bounds 12 -- linesearch (LS) - Armijo rule on direction 13 14 Many other reformulations are possible and combinations of 15 feasible/infeasible and linesearch/trust region are possible. 16 17 Basic theory 18 Fischer-Burmeister reformulation is semismooth with a continuously 19 differentiable merit function and strongly semismooth if the F has 20 lipschitz continuous derivatives. 21 22 Every accumulation point generated by the algorithm is a stationary 23 point for the merit function. Stationary points of the merit function 24 are solutions of the complementarity problem if 25 a. the stationary point has a BD-regular subdifferential, or 26 b. the Schur complement F'/F'_ff is a P_0-matrix where ff is the 27 index set corresponding to the free variables. 28 29 If one of the accumulation points has a BD-regular subdifferential then 30 a. the entire sequence converges to this accumulation point at 31 a local q-superlinear rate 32 b. if in addition the reformulation is strongly semismooth near 33 this accumulation point, then the algorithm converges at a 34 local q-quadratic rate. 35 36 The theory for the feasible version follows from the feasible descent 37 algorithm framework. 38 39 References: 40 + * - Billups, "Algorithms for Complementarity Problems and Generalized 41 Equations," Ph.D thesis, University of Wisconsin Madison, 1995. 42 . * - De Luca, Facchinei, Kanzow, "A Semismooth Equation Approach to the 43 Solution of Nonlinear Complementarity Problems," Mathematical 44 Programming, 75, pages 407439, 1996. 45 . * - Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed 46 Complementarity Problems," Mathematical Programming, 86, 47 pages 475497, 1999. 48 . * - Fischer, "A Special Newton type Optimization Method," Optimization, 49 24, 1992 50 - * - Munson, Facchinei, Ferris, Fischer, Kanzow, "The Semismooth Algorithm 51 for Large Scale Complementarity Problems," Technical Report, 52 University of Wisconsin Madison, 1999. 53 */ 54 55 static PetscErrorCode TaoSetUp_ASFLS(Tao tao) 56 { 57 TAO_SSLS *asls = (TAO_SSLS *)tao->data; 58 59 PetscFunctionBegin; 60 CHKERRQ(VecDuplicate(tao->solution,&tao->gradient)); 61 CHKERRQ(VecDuplicate(tao->solution,&tao->stepdirection)); 62 CHKERRQ(VecDuplicate(tao->solution,&asls->ff)); 63 CHKERRQ(VecDuplicate(tao->solution,&asls->dpsi)); 64 CHKERRQ(VecDuplicate(tao->solution,&asls->da)); 65 CHKERRQ(VecDuplicate(tao->solution,&asls->db)); 66 CHKERRQ(VecDuplicate(tao->solution,&asls->t1)); 67 CHKERRQ(VecDuplicate(tao->solution,&asls->t2)); 68 CHKERRQ(VecDuplicate(tao->solution, &asls->w)); 69 asls->fixed = NULL; 70 asls->free = NULL; 71 asls->J_sub = NULL; 72 asls->Jpre_sub = NULL; 73 asls->r1 = NULL; 74 asls->r2 = NULL; 75 asls->r3 = NULL; 76 asls->dxfree = NULL; 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode Tao_ASLS_FunctionGradient(TaoLineSearch ls, Vec X, PetscReal *fcn, Vec G, void *ptr) 81 { 82 Tao tao = (Tao)ptr; 83 TAO_SSLS *asls = (TAO_SSLS *)tao->data; 84 85 PetscFunctionBegin; 86 CHKERRQ(TaoComputeConstraints(tao, X, tao->constraints)); 87 CHKERRQ(VecFischer(X,tao->constraints,tao->XL,tao->XU,asls->ff)); 88 CHKERRQ(VecNorm(asls->ff,NORM_2,&asls->merit)); 89 *fcn = 0.5*asls->merit*asls->merit; 90 CHKERRQ(TaoComputeJacobian(tao,tao->solution,tao->jacobian,tao->jacobian_pre)); 91 92 CHKERRQ(MatDFischer(tao->jacobian, tao->solution, tao->constraints,tao->XL, tao->XU, asls->t1, asls->t2,asls->da, asls->db)); 93 CHKERRQ(VecPointwiseMult(asls->t1, asls->ff, asls->db)); 94 CHKERRQ(MatMultTranspose(tao->jacobian,asls->t1,G)); 95 CHKERRQ(VecPointwiseMult(asls->t1, asls->ff, asls->da)); 96 CHKERRQ(VecAXPY(G,1.0,asls->t1)); 97 PetscFunctionReturn(0); 98 } 99 100 static PetscErrorCode TaoDestroy_ASFLS(Tao tao) 101 { 102 TAO_SSLS *ssls = (TAO_SSLS *)tao->data; 103 104 PetscFunctionBegin; 105 CHKERRQ(VecDestroy(&ssls->ff)); 106 CHKERRQ(VecDestroy(&ssls->dpsi)); 107 CHKERRQ(VecDestroy(&ssls->da)); 108 CHKERRQ(VecDestroy(&ssls->db)); 109 CHKERRQ(VecDestroy(&ssls->w)); 110 CHKERRQ(VecDestroy(&ssls->t1)); 111 CHKERRQ(VecDestroy(&ssls->t2)); 112 CHKERRQ(VecDestroy(&ssls->r1)); 113 CHKERRQ(VecDestroy(&ssls->r2)); 114 CHKERRQ(VecDestroy(&ssls->r3)); 115 CHKERRQ(VecDestroy(&ssls->dxfree)); 116 CHKERRQ(MatDestroy(&ssls->J_sub)); 117 CHKERRQ(MatDestroy(&ssls->Jpre_sub)); 118 CHKERRQ(ISDestroy(&ssls->fixed)); 119 CHKERRQ(ISDestroy(&ssls->free)); 120 CHKERRQ(PetscFree(tao->data)); 121 tao->data = NULL; 122 PetscFunctionReturn(0); 123 } 124 125 static PetscErrorCode TaoSolve_ASFLS(Tao tao) 126 { 127 TAO_SSLS *asls = (TAO_SSLS *)tao->data; 128 PetscReal psi,ndpsi, normd, innerd, t=0; 129 PetscInt nf; 130 TaoLineSearchConvergedReason ls_reason; 131 132 PetscFunctionBegin; 133 /* Assume that Setup has been called! 134 Set the structure for the Jacobian and create a linear solver. */ 135 136 CHKERRQ(TaoComputeVariableBounds(tao)); 137 CHKERRQ(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_ASLS_FunctionGradient,tao)); 138 CHKERRQ(TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao)); 139 CHKERRQ(TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU)); 140 141 CHKERRQ(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution)); 142 143 /* Calculate the function value and fischer function value at the 144 current iterate */ 145 CHKERRQ(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi)); 146 CHKERRQ(VecNorm(asls->dpsi,NORM_2,&ndpsi)); 147 148 tao->reason = TAO_CONTINUE_ITERATING; 149 while (1) { 150 /* Check the converged criteria */ 151 CHKERRQ(PetscInfo(tao,"iter %D, merit: %g, ||dpsi||: %g\n",tao->niter,(double)asls->merit,(double)ndpsi)); 152 CHKERRQ(TaoLogConvergenceHistory(tao,asls->merit,ndpsi,0.0,tao->ksp_its)); 153 CHKERRQ(TaoMonitor(tao,tao->niter,asls->merit,ndpsi,0.0,t)); 154 CHKERRQ((*tao->ops->convergencetest)(tao,tao->cnvP)); 155 if (TAO_CONTINUE_ITERATING != tao->reason) break; 156 157 /* Call general purpose update function */ 158 if (tao->ops->update) { 159 CHKERRQ((*tao->ops->update)(tao, tao->niter, tao->user_update)); 160 } 161 tao->niter++; 162 163 /* We are going to solve a linear system of equations. We need to 164 set the tolerances for the solve so that we maintain an asymptotic 165 rate of convergence that is superlinear. 166 Note: these tolerances are for the reduced system. We really need 167 to make sure that the full system satisfies the full-space conditions. 168 169 This rule gives superlinear asymptotic convergence 170 asls->atol = min(0.5, asls->merit*sqrt(asls->merit)); 171 asls->rtol = 0.0; 172 173 This rule gives quadratic asymptotic convergence 174 asls->atol = min(0.5, asls->merit*asls->merit); 175 asls->rtol = 0.0; 176 177 Calculate a free and fixed set of variables. The fixed set of 178 variables are those for the d_b is approximately equal to zero. 179 The definition of approximately changes as we approach the solution 180 to the problem. 181 182 No one rule is guaranteed to work in all cases. The following 183 definition is based on the norm of the Jacobian matrix. If the 184 norm is large, the tolerance becomes smaller. */ 185 CHKERRQ(MatNorm(tao->jacobian,NORM_1,&asls->identifier)); 186 asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier); 187 188 CHKERRQ(VecSet(asls->t1,-asls->identifier)); 189 CHKERRQ(VecSet(asls->t2, asls->identifier)); 190 191 CHKERRQ(ISDestroy(&asls->fixed)); 192 CHKERRQ(ISDestroy(&asls->free)); 193 CHKERRQ(VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed)); 194 CHKERRQ(ISComplementVec(asls->fixed,asls->t1, &asls->free)); 195 196 CHKERRQ(ISGetSize(asls->fixed,&nf)); 197 CHKERRQ(PetscInfo(tao,"Number of fixed variables: %D\n", nf)); 198 199 /* We now have our partition. Now calculate the direction in the 200 fixed variable space. */ 201 CHKERRQ(TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1)); 202 CHKERRQ(TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2)); 203 CHKERRQ(VecPointwiseDivide(asls->r1,asls->r1,asls->r2)); 204 CHKERRQ(VecSet(tao->stepdirection,0.0)); 205 CHKERRQ(VecISAXPY(tao->stepdirection, asls->fixed, 1.0,asls->r1)); 206 207 /* Our direction in the Fixed Variable Set is fixed. Calculate the 208 information needed for the step in the Free Variable Set. To 209 do this, we need to know the diagonal perturbation and the 210 right hand side. */ 211 212 CHKERRQ(TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1)); 213 CHKERRQ(TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2)); 214 CHKERRQ(TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3)); 215 CHKERRQ(VecPointwiseDivide(asls->r1,asls->r1, asls->r3)); 216 CHKERRQ(VecPointwiseDivide(asls->r2,asls->r2, asls->r3)); 217 218 /* r1 is the diagonal perturbation 219 r2 is the right hand side 220 r3 is no longer needed 221 222 Now need to modify r2 for our direction choice in the fixed 223 variable set: calculate t1 = J*d, take the reduced vector 224 of t1 and modify r2. */ 225 226 CHKERRQ(MatMult(tao->jacobian, tao->stepdirection, asls->t1)); 227 CHKERRQ(TaoVecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3)); 228 CHKERRQ(VecAXPY(asls->r2, -1.0, asls->r3)); 229 230 /* Calculate the reduced problem matrix and the direction */ 231 CHKERRQ(TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub)); 232 if (tao->jacobian != tao->jacobian_pre) { 233 CHKERRQ(TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub)); 234 } else { 235 CHKERRQ(MatDestroy(&asls->Jpre_sub)); 236 asls->Jpre_sub = asls->J_sub; 237 CHKERRQ(PetscObjectReference((PetscObject)(asls->Jpre_sub))); 238 } 239 CHKERRQ(MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES)); 240 CHKERRQ(TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree)); 241 CHKERRQ(VecSet(asls->dxfree, 0.0)); 242 243 /* Calculate the reduced direction. (Really negative of Newton 244 direction. Therefore, rest of the code uses -d.) */ 245 CHKERRQ(KSPReset(tao->ksp)); 246 CHKERRQ(KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub)); 247 CHKERRQ(KSPSolve(tao->ksp, asls->r2, asls->dxfree)); 248 CHKERRQ(KSPGetIterationNumber(tao->ksp,&tao->ksp_its)); 249 tao->ksp_tot_its+=tao->ksp_its; 250 251 /* Add the direction in the free variables back into the real direction. */ 252 CHKERRQ(VecISAXPY(tao->stepdirection, asls->free, 1.0,asls->dxfree)); 253 254 /* Check the projected real direction for descent and if not, use the negative 255 gradient direction. */ 256 CHKERRQ(VecCopy(tao->stepdirection, asls->w)); 257 CHKERRQ(VecScale(asls->w, -1.0)); 258 CHKERRQ(VecBoundGradientProjection(asls->w, tao->solution, tao->XL, tao->XU, asls->w)); 259 CHKERRQ(VecNorm(asls->w, NORM_2, &normd)); 260 CHKERRQ(VecDot(asls->w, asls->dpsi, &innerd)); 261 262 if (innerd >= -asls->delta*PetscPowReal(normd, asls->rho)) { 263 CHKERRQ(PetscInfo(tao,"Gradient direction: %5.4e.\n", (double)innerd)); 264 CHKERRQ(PetscInfo(tao, "Iteration %D: newton direction not descent\n", tao->niter)); 265 CHKERRQ(VecCopy(asls->dpsi, tao->stepdirection)); 266 CHKERRQ(VecDot(asls->dpsi, tao->stepdirection, &innerd)); 267 } 268 269 CHKERRQ(VecScale(tao->stepdirection, -1.0)); 270 innerd = -innerd; 271 272 /* We now have a correct descent direction. Apply a linesearch to 273 find the new iterate. */ 274 CHKERRQ(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0)); 275 CHKERRQ(TaoLineSearchApply(tao->linesearch, tao->solution, &psi,asls->dpsi, tao->stepdirection, &t, &ls_reason)); 276 CHKERRQ(VecNorm(asls->dpsi, NORM_2, &ndpsi)); 277 } 278 PetscFunctionReturn(0); 279 } 280 281 /* ---------------------------------------------------------- */ 282 /*MC 283 TAOASFLS - Active-set feasible linesearch algorithm for solving 284 complementarity constraints 285 286 Options Database Keys: 287 + -tao_ssls_delta - descent test fraction 288 - -tao_ssls_rho - descent test power 289 290 Level: beginner 291 M*/ 292 PETSC_EXTERN PetscErrorCode TaoCreate_ASFLS(Tao tao) 293 { 294 TAO_SSLS *asls; 295 const char *armijo_type = TAOLINESEARCHARMIJO; 296 297 PetscFunctionBegin; 298 CHKERRQ(PetscNewLog(tao,&asls)); 299 tao->data = (void*)asls; 300 tao->ops->solve = TaoSolve_ASFLS; 301 tao->ops->setup = TaoSetUp_ASFLS; 302 tao->ops->view = TaoView_SSLS; 303 tao->ops->setfromoptions = TaoSetFromOptions_SSLS; 304 tao->ops->destroy = TaoDestroy_ASFLS; 305 tao->subset_type = TAO_SUBSET_SUBVEC; 306 asls->delta = 1e-10; 307 asls->rho = 2.1; 308 asls->fixed = NULL; 309 asls->free = NULL; 310 asls->J_sub = NULL; 311 asls->Jpre_sub = NULL; 312 asls->w = NULL; 313 asls->r1 = NULL; 314 asls->r2 = NULL; 315 asls->r3 = NULL; 316 asls->t1 = NULL; 317 asls->t2 = NULL; 318 asls->dxfree = NULL; 319 asls->identifier = 1e-5; 320 321 CHKERRQ(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 322 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 323 CHKERRQ(TaoLineSearchSetType(tao->linesearch, armijo_type)); 324 CHKERRQ(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 325 CHKERRQ(TaoLineSearchSetFromOptions(tao->linesearch)); 326 327 CHKERRQ(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 328 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 329 CHKERRQ(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 330 CHKERRQ(KSPSetFromOptions(tao->ksp)); 331 332 /* Override default settings (unless already changed) */ 333 if (!tao->max_it_changed) tao->max_it = 2000; 334 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 335 if (!tao->gttol_changed) tao->gttol = 0; 336 if (!tao->grtol_changed) tao->grtol = 0; 337 #if defined(PETSC_USE_REAL_SINGLE) 338 if (!tao->gatol_changed) tao->gatol = 1.0e-6; 339 if (!tao->fmin_changed) tao->fmin = 1.0e-4; 340 #else 341 if (!tao->gatol_changed) tao->gatol = 1.0e-16; 342 if (!tao->fmin_changed) tao->fmin = 1.0e-8; 343 #endif 344 PetscFunctionReturn(0); 345 } 346