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