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