1 2 #include <../src/snes/impls/vi/ss/vissimpl.h> /*I "petscsnes.h" I*/ 3 4 /* 5 SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem. 6 7 Input Parameter: 8 . phi - the semismooth function 9 10 Output Parameter: 11 . merit - the merit function 12 . phinorm - ||phi|| 13 14 Notes: 15 The merit function for the mixed complementarity problem is defined as 16 merit = 0.5*phi^T*phi 17 */ 18 static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit, PetscReal *phinorm) { 19 PetscFunctionBegin; 20 PetscCall(VecNormBegin(phi, NORM_2, phinorm)); 21 PetscCall(VecNormEnd(phi, NORM_2, phinorm)); 22 23 *merit = 0.5 * (*phinorm) * (*phinorm); 24 PetscFunctionReturn(0); 25 } 26 27 static inline PetscScalar Phi(PetscScalar a, PetscScalar b) { 28 return a + b - PetscSqrtScalar(a * a + b * b); 29 } 30 31 static inline PetscScalar DPhi(PetscScalar a, PetscScalar b) { 32 if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a / PetscSqrtScalar(a * a + b * b); 33 else return .5; 34 } 35 36 /* 37 SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form. 38 39 Input Parameters: 40 . snes - the SNES context 41 . X - current iterate 42 . functx - user defined function context 43 44 Output Parameters: 45 . phi - Semismooth function 46 47 */ 48 static PetscErrorCode SNESVIComputeFunction(SNES snes, Vec X, Vec phi, void *functx) { 49 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data; 50 Vec Xl = snes->xl, Xu = snes->xu, F = snes->vec_func; 51 PetscScalar *phi_arr, *f_arr, *l, *u; 52 const PetscScalar *x_arr; 53 PetscInt i, nlocal; 54 55 PetscFunctionBegin; 56 PetscCall((*vi->computeuserfunction)(snes, X, F, functx)); 57 PetscCall(VecGetLocalSize(X, &nlocal)); 58 PetscCall(VecGetArrayRead(X, &x_arr)); 59 PetscCall(VecGetArray(F, &f_arr)); 60 PetscCall(VecGetArray(Xl, &l)); 61 PetscCall(VecGetArray(Xu, &u)); 62 PetscCall(VecGetArray(phi, &phi_arr)); 63 64 for (i = 0; i < nlocal; i++) { 65 if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */ 66 phi_arr[i] = f_arr[i]; 67 } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */ 68 phi_arr[i] = -Phi(u[i] - x_arr[i], -f_arr[i]); 69 } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */ 70 phi_arr[i] = Phi(x_arr[i] - l[i], f_arr[i]); 71 } else if (l[i] == u[i]) { 72 phi_arr[i] = l[i] - x_arr[i]; 73 } else { /* both bounds on variable */ 74 phi_arr[i] = Phi(x_arr[i] - l[i], -Phi(u[i] - x_arr[i], -f_arr[i])); 75 } 76 } 77 78 PetscCall(VecRestoreArrayRead(X, &x_arr)); 79 PetscCall(VecRestoreArray(F, &f_arr)); 80 PetscCall(VecRestoreArray(Xl, &l)); 81 PetscCall(VecRestoreArray(Xu, &u)); 82 PetscCall(VecRestoreArray(phi, &phi_arr)); 83 PetscFunctionReturn(0); 84 } 85 86 /* 87 SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the 88 the semismooth jacobian. 89 */ 90 PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes, Vec X, Vec F, Mat jac, Vec Da, Vec Db) { 91 PetscScalar *l, *u, *x, *f, *da, *db, da1, da2, db1, db2; 92 PetscInt i, nlocal; 93 94 PetscFunctionBegin; 95 PetscCall(VecGetArray(X, &x)); 96 PetscCall(VecGetArray(F, &f)); 97 PetscCall(VecGetArray(snes->xl, &l)); 98 PetscCall(VecGetArray(snes->xu, &u)); 99 PetscCall(VecGetArray(Da, &da)); 100 PetscCall(VecGetArray(Db, &db)); 101 PetscCall(VecGetLocalSize(X, &nlocal)); 102 103 for (i = 0; i < nlocal; i++) { 104 if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */ 105 da[i] = 0; 106 db[i] = 1; 107 } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */ 108 da[i] = DPhi(u[i] - x[i], -f[i]); 109 db[i] = DPhi(-f[i], u[i] - x[i]); 110 } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */ 111 da[i] = DPhi(x[i] - l[i], f[i]); 112 db[i] = DPhi(f[i], x[i] - l[i]); 113 } else if (l[i] == u[i]) { /* fixed variable */ 114 da[i] = 1; 115 db[i] = 0; 116 } else { /* upper and lower bounds on variable */ 117 da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i])); 118 db1 = DPhi(-Phi(u[i] - x[i], -f[i]), x[i] - l[i]); 119 da2 = DPhi(u[i] - x[i], -f[i]); 120 db2 = DPhi(-f[i], u[i] - x[i]); 121 da[i] = da1 + db1 * da2; 122 db[i] = db1 * db2; 123 } 124 } 125 126 PetscCall(VecRestoreArray(X, &x)); 127 PetscCall(VecRestoreArray(F, &f)); 128 PetscCall(VecRestoreArray(snes->xl, &l)); 129 PetscCall(VecRestoreArray(snes->xu, &u)); 130 PetscCall(VecRestoreArray(Da, &da)); 131 PetscCall(VecRestoreArray(Db, &db)); 132 PetscFunctionReturn(0); 133 } 134 135 /* 136 SNESVIComputeJacobian - Computes the jacobian of the semismooth function.The Jacobian for the semismooth function is an element of the B-subdifferential of the Fischer-Burmeister function for complementarity problems. 137 138 Input Parameters: 139 . Da - Diagonal shift vector for the semismooth jacobian. 140 . Db - Row scaling vector for the semismooth jacobian. 141 142 Output Parameters: 143 . jac - semismooth jacobian 144 . jac_pre - optional preconditioning matrix 145 146 Notes: 147 The semismooth jacobian matrix is given by 148 jac = Da + Db*jacfun 149 where Db is the row scaling matrix stored as a vector, 150 Da is the diagonal perturbation matrix stored as a vector 151 and jacfun is the jacobian of the original nonlinear function. 152 */ 153 PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre, Vec Da, Vec Db) { 154 /* Do row scaling and add diagonal perturbation */ 155 PetscFunctionBegin; 156 PetscCall(MatDiagonalScale(jac, Db, NULL)); 157 PetscCall(MatDiagonalSet(jac, Da, ADD_VALUES)); 158 if (jac != jac_pre) { /* If jac and jac_pre are different */ 159 PetscCall(MatDiagonalScale(jac_pre, Db, NULL)); 160 PetscCall(MatDiagonalSet(jac_pre, Da, ADD_VALUES)); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 /* 166 SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi. 167 168 Input Parameters: 169 phi - semismooth function. 170 H - semismooth jacobian 171 172 Output Parameters: 173 dpsi - merit function gradient 174 175 Notes: 176 The merit function gradient is computed as follows 177 dpsi = H^T*phi 178 */ 179 PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi) { 180 PetscFunctionBegin; 181 PetscCall(MatMultTranspose(H, phi, dpsi)); 182 PetscFunctionReturn(0); 183 } 184 185 /* 186 SNESSolve_VINEWTONSSLS - Solves the complementarity problem with a semismooth Newton 187 method using a line search. 188 189 Input Parameters: 190 . snes - the SNES context 191 192 Application Interface Routine: SNESSolve() 193 194 Notes: 195 This implements essentially a semismooth Newton method with a 196 line search. The default line search does not do any line search 197 but rather takes a full Newton step. 198 199 Developer Note: the code in this file should be slightly modified so that this routine need not exist and the SNESSolve_NEWTONLS() routine is called directly with the appropriate wrapped function and Jacobian evaluations 200 201 */ 202 PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes) { 203 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data; 204 PetscInt maxits, i, lits; 205 SNESLineSearchReason lssucceed; 206 PetscReal gnorm, xnorm = 0, ynorm; 207 Vec Y, X, F; 208 KSPConvergedReason kspreason; 209 DM dm; 210 DMSNES sdm; 211 212 PetscFunctionBegin; 213 PetscCall(SNESGetDM(snes, &dm)); 214 PetscCall(DMGetDMSNES(dm, &sdm)); 215 216 vi->computeuserfunction = sdm->ops->computefunction; 217 sdm->ops->computefunction = SNESVIComputeFunction; 218 219 snes->numFailures = 0; 220 snes->numLinearSolveFailures = 0; 221 snes->reason = SNES_CONVERGED_ITERATING; 222 223 maxits = snes->max_its; /* maximum number of iterations */ 224 X = snes->vec_sol; /* solution vector */ 225 F = snes->vec_func; /* residual vector */ 226 Y = snes->work[0]; /* work vectors */ 227 228 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 229 snes->iter = 0; 230 snes->norm = 0.0; 231 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 232 233 PetscCall(SNESVIProjectOntoBounds(snes, X)); 234 PetscCall(SNESComputeFunction(snes, X, vi->phi)); 235 if (snes->domainerror) { 236 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 237 sdm->ops->computefunction = vi->computeuserfunction; 238 PetscFunctionReturn(0); 239 } 240 /* Compute Merit function */ 241 PetscCall(SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm)); 242 243 PetscCall(VecNormBegin(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */ 244 PetscCall(VecNormEnd(X, NORM_2, &xnorm)); 245 SNESCheckFunctionNorm(snes, vi->merit); 246 247 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 248 snes->norm = vi->phinorm; 249 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 250 PetscCall(SNESLogConvergenceHistory(snes, vi->phinorm, 0)); 251 PetscCall(SNESMonitor(snes, 0, vi->phinorm)); 252 253 /* test convergence */ 254 PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, vi->phinorm, &snes->reason, snes->cnvP); 255 if (snes->reason) { 256 sdm->ops->computefunction = vi->computeuserfunction; 257 PetscFunctionReturn(0); 258 } 259 260 for (i = 0; i < maxits; i++) { 261 /* Call general purpose update function */ 262 PetscTryTypeMethod(snes, update, snes->iter); 263 264 /* Solve J Y = Phi, where J is the semismooth jacobian */ 265 266 /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/ 267 sdm->ops->computefunction = vi->computeuserfunction; 268 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre)); 269 SNESCheckJacobianDomainerror(snes); 270 sdm->ops->computefunction = SNESVIComputeFunction; 271 272 /* Get the diagonal shift and row scaling vectors */ 273 PetscCall(SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db)); 274 /* Compute the semismooth jacobian */ 275 PetscCall(SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db)); 276 /* Compute the merit function gradient */ 277 PetscCall(SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi)); 278 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre)); 279 PetscCall(KSPSolve(snes->ksp, vi->phi, Y)); 280 PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason)); 281 282 if (kspreason < 0) { 283 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 284 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures)); 285 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 286 break; 287 } 288 } 289 PetscCall(KSPGetIterationNumber(snes->ksp, &lits)); 290 snes->linear_its += lits; 291 PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits)); 292 /* 293 if (snes->ops->precheck) { 294 PetscBool changed_y = PETSC_FALSE; 295 PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y); 296 } 297 298 if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W)); 299 */ 300 /* Compute a (scaled) negative update in the line search routine: 301 Y <- X - lambda*Y 302 and evaluate G = function(Y) (depends on the line search). 303 */ 304 PetscCall(VecCopy(Y, snes->vec_sol_update)); 305 ynorm = 1; 306 gnorm = vi->phinorm; 307 PetscCall(SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y)); 308 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed)); 309 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm)); 310 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)vi->phinorm, (double)gnorm, (double)ynorm, (int)lssucceed)); 311 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 312 if (snes->domainerror) { 313 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 314 sdm->ops->computefunction = vi->computeuserfunction; 315 PetscFunctionReturn(0); 316 } 317 if (lssucceed) { 318 if (++snes->numFailures >= snes->maxFailures) { 319 PetscBool ismin; 320 snes->reason = SNES_DIVERGED_LINE_SEARCH; 321 PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, vi->phi, X, gnorm, &ismin)); 322 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 323 break; 324 } 325 } 326 /* Update function and solution vectors */ 327 vi->phinorm = gnorm; 328 vi->merit = 0.5 * vi->phinorm * vi->phinorm; 329 /* Monitor convergence */ 330 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 331 snes->iter = i + 1; 332 snes->norm = vi->phinorm; 333 snes->xnorm = xnorm; 334 snes->ynorm = ynorm; 335 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 336 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits)); 337 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 338 /* Test for convergence, xnorm = || X || */ 339 if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm)); 340 PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, vi->phinorm, &snes->reason, snes->cnvP); 341 if (snes->reason) break; 342 } 343 if (i == maxits) { 344 PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits)); 345 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 346 } 347 sdm->ops->computefunction = vi->computeuserfunction; 348 PetscFunctionReturn(0); 349 } 350 351 /* -------------------------------------------------------------------------- */ 352 /* 353 SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use 354 of the SNES nonlinear solver. 355 356 Input Parameter: 357 . snes - the SNES context 358 359 Application Interface Routine: SNESSetUp() 360 361 Notes: 362 For basic use of the SNES solvers, the user need not explicitly call 363 SNESSetUp(), since these actions will automatically occur during 364 the call to SNESSolve(). 365 */ 366 PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes) { 367 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data; 368 369 PetscFunctionBegin; 370 PetscCall(SNESSetUp_VI(snes)); 371 PetscCall(VecDuplicate(snes->vec_sol, &vi->dpsi)); 372 PetscCall(VecDuplicate(snes->vec_sol, &vi->phi)); 373 PetscCall(VecDuplicate(snes->vec_sol, &vi->Da)); 374 PetscCall(VecDuplicate(snes->vec_sol, &vi->Db)); 375 PetscCall(VecDuplicate(snes->vec_sol, &vi->z)); 376 PetscCall(VecDuplicate(snes->vec_sol, &vi->t)); 377 PetscFunctionReturn(0); 378 } 379 /* -------------------------------------------------------------------------- */ 380 PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes) { 381 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data; 382 383 PetscFunctionBegin; 384 PetscCall(SNESReset_VI(snes)); 385 PetscCall(VecDestroy(&vi->dpsi)); 386 PetscCall(VecDestroy(&vi->phi)); 387 PetscCall(VecDestroy(&vi->Da)); 388 PetscCall(VecDestroy(&vi->Db)); 389 PetscCall(VecDestroy(&vi->z)); 390 PetscCall(VecDestroy(&vi->t)); 391 PetscFunctionReturn(0); 392 } 393 394 /* -------------------------------------------------------------------------- */ 395 /* 396 SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method. 397 398 Input Parameter: 399 . snes - the SNES context 400 401 Application Interface Routine: SNESSetFromOptions() 402 */ 403 static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject) { 404 PetscFunctionBegin; 405 PetscCall(SNESSetFromOptions_VI(snes, PetscOptionsObject)); 406 PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options"); 407 PetscOptionsHeadEnd(); 408 PetscFunctionReturn(0); 409 } 410 411 /* -------------------------------------------------------------------------- */ 412 /*MC 413 SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method 414 415 Options Database: 416 + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method 417 - -snes_vi_monitor - prints the number of active constraints at each iteration. 418 419 Level: beginner 420 421 References: 422 + * - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth 423 algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001). 424 - * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 425 Applications, Optimization Methods and Software, 21 (2006). 426 427 .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` 428 429 M*/ 430 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes) { 431 SNES_VINEWTONSSLS *vi; 432 SNESLineSearch linesearch; 433 434 PetscFunctionBegin; 435 snes->ops->reset = SNESReset_VINEWTONSSLS; 436 snes->ops->setup = SNESSetUp_VINEWTONSSLS; 437 snes->ops->solve = SNESSolve_VINEWTONSSLS; 438 snes->ops->destroy = SNESDestroy_VI; 439 snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS; 440 snes->ops->view = NULL; 441 442 snes->usesksp = PETSC_TRUE; 443 snes->usesnpc = PETSC_FALSE; 444 445 PetscCall(SNESGetLineSearch(snes, &linesearch)); 446 if (!((PetscObject)linesearch)->type_name) { 447 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 448 PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0)); 449 } 450 451 snes->alwayscomputesfinalresidual = PETSC_FALSE; 452 453 PetscCall(PetscNewLog(snes, &vi)); 454 snes->data = (void *)vi; 455 456 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI)); 457 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI)); 458 PetscFunctionReturn(0); 459 } 460