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