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