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