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, 1996. 45 . * - Ferris, Kanzow, Munson, "Feasible Descent Algorithms for Mixed 46 Complementarity Problems," Mathematical Programming, 86, 47 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_ASILS(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 asls->fixed = NULL; 69 asls->free = NULL; 70 asls->J_sub = NULL; 71 asls->Jpre_sub = NULL; 72 asls->w = 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 91 CHKERRQ(TaoComputeJacobian(tao,tao->solution,tao->jacobian,tao->jacobian_pre)); 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_ASILS(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 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode TaoSolve_ASILS(Tao tao) 125 { 126 TAO_SSLS *asls = (TAO_SSLS *)tao->data; 127 PetscReal psi,ndpsi, normd, innerd, t=0; 128 PetscInt nf; 129 TaoLineSearchConvergedReason ls_reason; 130 131 PetscFunctionBegin; 132 /* Assume that Setup has been called! 133 Set the structure for the Jacobian and create a linear solver. */ 134 135 CHKERRQ(TaoComputeVariableBounds(tao)); 136 CHKERRQ(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_ASLS_FunctionGradient,tao)); 137 CHKERRQ(TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao)); 138 139 /* Calculate the function value and fischer function value at the 140 current iterate */ 141 CHKERRQ(TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,asls->dpsi)); 142 CHKERRQ(VecNorm(asls->dpsi,NORM_2,&ndpsi)); 143 144 tao->reason = TAO_CONTINUE_ITERATING; 145 while (1) { 146 /* Check the termination criteria */ 147 CHKERRQ(PetscInfo(tao,"iter %D, merit: %g, ||dpsi||: %g\n",tao->niter, (double)asls->merit, (double)ndpsi)); 148 CHKERRQ(TaoLogConvergenceHistory(tao,asls->merit,ndpsi,0.0,tao->ksp_its)); 149 CHKERRQ(TaoMonitor(tao,tao->niter,asls->merit,ndpsi,0.0,t)); 150 CHKERRQ((*tao->ops->convergencetest)(tao,tao->cnvP)); 151 if (TAO_CONTINUE_ITERATING != tao->reason) break; 152 153 /* Call general purpose update function */ 154 if (tao->ops->update) { 155 CHKERRQ((*tao->ops->update)(tao, tao->niter, tao->user_update)); 156 } 157 tao->niter++; 158 159 /* We are going to solve a linear system of equations. We need to 160 set the tolerances for the solve so that we maintain an asymptotic 161 rate of convergence that is superlinear. 162 Note: these tolerances are for the reduced system. We really need 163 to make sure that the full system satisfies the full-space conditions. 164 165 This rule gives superlinear asymptotic convergence 166 asls->atol = min(0.5, asls->merit*sqrt(asls->merit)); 167 asls->rtol = 0.0; 168 169 This rule gives quadratic asymptotic convergence 170 asls->atol = min(0.5, asls->merit*asls->merit); 171 asls->rtol = 0.0; 172 173 Calculate a free and fixed set of variables. The fixed set of 174 variables are those for the d_b is approximately equal to zero. 175 The definition of approximately changes as we approach the solution 176 to the problem. 177 178 No one rule is guaranteed to work in all cases. The following 179 definition is based on the norm of the Jacobian matrix. If the 180 norm is large, the tolerance becomes smaller. */ 181 CHKERRQ(MatNorm(tao->jacobian,NORM_1,&asls->identifier)); 182 asls->identifier = PetscMin(asls->merit, 1e-2) / (1 + asls->identifier); 183 184 CHKERRQ(VecSet(asls->t1,-asls->identifier)); 185 CHKERRQ(VecSet(asls->t2, asls->identifier)); 186 187 CHKERRQ(ISDestroy(&asls->fixed)); 188 CHKERRQ(ISDestroy(&asls->free)); 189 CHKERRQ(VecWhichBetweenOrEqual(asls->t1, asls->db, asls->t2, &asls->fixed)); 190 CHKERRQ(ISComplementVec(asls->fixed,asls->t1, &asls->free)); 191 192 CHKERRQ(ISGetSize(asls->fixed,&nf)); 193 CHKERRQ(PetscInfo(tao,"Number of fixed variables: %D\n", nf)); 194 195 /* We now have our partition. Now calculate the direction in the 196 fixed variable space. */ 197 CHKERRQ(TaoVecGetSubVec(asls->ff, asls->fixed, tao->subset_type, 0.0, &asls->r1)); 198 CHKERRQ(TaoVecGetSubVec(asls->da, asls->fixed, tao->subset_type, 1.0, &asls->r2)); 199 CHKERRQ(VecPointwiseDivide(asls->r1,asls->r1,asls->r2)); 200 CHKERRQ(VecSet(tao->stepdirection,0.0)); 201 CHKERRQ(VecISAXPY(tao->stepdirection, asls->fixed,1.0,asls->r1)); 202 203 /* Our direction in the Fixed Variable Set is fixed. Calculate the 204 information needed for the step in the Free Variable Set. To 205 do this, we need to know the diagonal perturbation and the 206 right hand side. */ 207 208 CHKERRQ(TaoVecGetSubVec(asls->da, asls->free, tao->subset_type, 0.0, &asls->r1)); 209 CHKERRQ(TaoVecGetSubVec(asls->ff, asls->free, tao->subset_type, 0.0, &asls->r2)); 210 CHKERRQ(TaoVecGetSubVec(asls->db, asls->free, tao->subset_type, 1.0, &asls->r3)); 211 CHKERRQ(VecPointwiseDivide(asls->r1,asls->r1, asls->r3)); 212 CHKERRQ(VecPointwiseDivide(asls->r2,asls->r2, asls->r3)); 213 214 /* r1 is the diagonal perturbation 215 r2 is the right hand side 216 r3 is no longer needed 217 218 Now need to modify r2 for our direction choice in the fixed 219 variable set: calculate t1 = J*d, take the reduced vector 220 of t1 and modify r2. */ 221 222 CHKERRQ(MatMult(tao->jacobian, tao->stepdirection, asls->t1)); 223 CHKERRQ(TaoVecGetSubVec(asls->t1,asls->free,tao->subset_type,0.0,&asls->r3)); 224 CHKERRQ(VecAXPY(asls->r2, -1.0, asls->r3)); 225 226 /* Calculate the reduced problem matrix and the direction */ 227 if (!asls->w && (tao->subset_type == TAO_SUBSET_MASK || tao->subset_type == TAO_SUBSET_MATRIXFREE)) { 228 CHKERRQ(VecDuplicate(tao->solution, &asls->w)); 229 } 230 CHKERRQ(TaoMatGetSubMat(tao->jacobian, asls->free, asls->w, tao->subset_type,&asls->J_sub)); 231 if (tao->jacobian != tao->jacobian_pre) { 232 CHKERRQ(TaoMatGetSubMat(tao->jacobian_pre, asls->free, asls->w, tao->subset_type, &asls->Jpre_sub)); 233 } else { 234 CHKERRQ(MatDestroy(&asls->Jpre_sub)); 235 asls->Jpre_sub = asls->J_sub; 236 CHKERRQ(PetscObjectReference((PetscObject)(asls->Jpre_sub))); 237 } 238 CHKERRQ(MatDiagonalSet(asls->J_sub, asls->r1,ADD_VALUES)); 239 CHKERRQ(TaoVecGetSubVec(tao->stepdirection, asls->free, tao->subset_type, 0.0, &asls->dxfree)); 240 CHKERRQ(VecSet(asls->dxfree, 0.0)); 241 242 /* Calculate the reduced direction. (Really negative of Newton 243 direction. Therefore, rest of the code uses -d.) */ 244 CHKERRQ(KSPReset(tao->ksp)); 245 CHKERRQ(KSPSetOperators(tao->ksp, asls->J_sub, asls->Jpre_sub)); 246 CHKERRQ(KSPSolve(tao->ksp, asls->r2, asls->dxfree)); 247 CHKERRQ(KSPGetIterationNumber(tao->ksp,&tao->ksp_its)); 248 tao->ksp_tot_its+=tao->ksp_its; 249 250 /* Add the direction in the free variables back into the real direction. */ 251 CHKERRQ(VecISAXPY(tao->stepdirection, asls->free, 1.0,asls->dxfree)); 252 253 /* Check the real direction for descent and if not, use the negative 254 gradient direction. */ 255 CHKERRQ(VecNorm(tao->stepdirection, NORM_2, &normd)); 256 CHKERRQ(VecDot(tao->stepdirection, asls->dpsi, &innerd)); 257 258 if (innerd <= asls->delta*PetscPowReal(normd, asls->rho)) { 259 CHKERRQ(PetscInfo(tao,"Gradient direction: %5.4e.\n", (double)innerd)); 260 CHKERRQ(PetscInfo(tao, "Iteration %D: newton direction not descent\n", tao->niter)); 261 CHKERRQ(VecCopy(asls->dpsi, tao->stepdirection)); 262 CHKERRQ(VecDot(asls->dpsi, tao->stepdirection, &innerd)); 263 } 264 265 CHKERRQ(VecScale(tao->stepdirection, -1.0)); 266 innerd = -innerd; 267 268 /* We now have a correct descent direction. Apply a linesearch to 269 find the new iterate. */ 270 CHKERRQ(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0)); 271 CHKERRQ(TaoLineSearchApply(tao->linesearch, tao->solution, &psi,asls->dpsi, tao->stepdirection, &t, &ls_reason)); 272 CHKERRQ(VecNorm(asls->dpsi, NORM_2, &ndpsi)); 273 } 274 PetscFunctionReturn(0); 275 } 276 277 /* ---------------------------------------------------------- */ 278 /*MC 279 TAOASILS - Active-set infeasible linesearch algorithm for solving 280 complementarity constraints 281 282 Options Database Keys: 283 + -tao_ssls_delta - descent test fraction 284 - -tao_ssls_rho - descent test power 285 286 Level: beginner 287 M*/ 288 PETSC_EXTERN PetscErrorCode TaoCreate_ASILS(Tao tao) 289 { 290 TAO_SSLS *asls; 291 const char *armijo_type = TAOLINESEARCHARMIJO; 292 293 PetscFunctionBegin; 294 CHKERRQ(PetscNewLog(tao,&asls)); 295 tao->data = (void*)asls; 296 tao->ops->solve = TaoSolve_ASILS; 297 tao->ops->setup = TaoSetUp_ASILS; 298 tao->ops->view = TaoView_SSLS; 299 tao->ops->setfromoptions = TaoSetFromOptions_SSLS; 300 tao->ops->destroy = TaoDestroy_ASILS; 301 tao->subset_type = TAO_SUBSET_SUBVEC; 302 asls->delta = 1e-10; 303 asls->rho = 2.1; 304 asls->fixed = NULL; 305 asls->free = NULL; 306 asls->J_sub = NULL; 307 asls->Jpre_sub = NULL; 308 asls->w = NULL; 309 asls->r1 = NULL; 310 asls->r2 = NULL; 311 asls->r3 = NULL; 312 asls->t1 = NULL; 313 asls->t2 = NULL; 314 asls->dxfree = NULL; 315 316 asls->identifier = 1e-5; 317 318 CHKERRQ(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 319 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 320 CHKERRQ(TaoLineSearchSetType(tao->linesearch, armijo_type)); 321 CHKERRQ(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 322 CHKERRQ(TaoLineSearchSetFromOptions(tao->linesearch)); 323 324 CHKERRQ(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 325 CHKERRQ(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 326 CHKERRQ(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 327 CHKERRQ(KSPSetFromOptions(tao->ksp)); 328 329 /* Override default settings (unless already changed) */ 330 if (!tao->max_it_changed) tao->max_it = 2000; 331 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 332 if (!tao->gttol_changed) tao->gttol = 0; 333 if (!tao->grtol_changed) tao->grtol = 0; 334 #if defined(PETSC_USE_REAL_SINGLE) 335 if (!tao->gatol_changed) tao->gatol = 1.0e-6; 336 if (!tao->fmin_changed) tao->fmin = 1.0e-4; 337 #else 338 if (!tao->gatol_changed) tao->gatol = 1.0e-16; 339 if (!tao->fmin_changed) tao->fmin = 1.0e-8; 340 #endif 341 PetscFunctionReturn(0); 342 } 343