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