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