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