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) { 273 PetscCall((*snes->ops->update)(snes, snes->iter)); 274 } 275 276 /* Solve J Y = Phi, where J is the semismooth jacobian */ 277 278 /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/ 279 sdm->ops->computefunction = vi->computeuserfunction; 280 PetscCall(SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre)); 281 SNESCheckJacobianDomainerror(snes); 282 sdm->ops->computefunction = SNESVIComputeFunction; 283 284 /* Get the diagonal shift and row scaling vectors */ 285 PetscCall(SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db)); 286 /* Compute the semismooth jacobian */ 287 PetscCall(SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db)); 288 /* Compute the merit function gradient */ 289 PetscCall(SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi)); 290 PetscCall(KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre)); 291 PetscCall(KSPSolve(snes->ksp,vi->phi,Y)); 292 PetscCall(KSPGetConvergedReason(snes->ksp,&kspreason)); 293 294 if (kspreason < 0) { 295 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 296 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures)); 297 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 298 break; 299 } 300 } 301 PetscCall(KSPGetIterationNumber(snes->ksp,&lits)); 302 snes->linear_its += lits; 303 PetscCall(PetscInfo(snes,"iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n",snes->iter,lits)); 304 /* 305 if (snes->ops->precheck) { 306 PetscBool changed_y = PETSC_FALSE; 307 PetscCall((*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y)); 308 } 309 310 if (PetscLogPrintInfo) { 311 PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W)); 312 } 313 */ 314 /* Compute a (scaled) negative update in the line search routine: 315 Y <- X - lambda*Y 316 and evaluate G = function(Y) (depends on the line search). 317 */ 318 PetscCall(VecCopy(Y,snes->vec_sol_update)); 319 ynorm = 1; gnorm = vi->phinorm; 320 PetscCall(SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y)); 321 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed)); 322 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm)); 323 PetscCall(PetscInfo(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)vi->phinorm,(double)gnorm,(double)ynorm,(int)lssucceed)); 324 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 325 if (snes->domainerror) { 326 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 327 sdm->ops->computefunction = vi->computeuserfunction; 328 PetscFunctionReturn(0); 329 } 330 if (lssucceed) { 331 if (++snes->numFailures >= snes->maxFailures) { 332 PetscBool ismin; 333 snes->reason = SNES_DIVERGED_LINE_SEARCH; 334 PetscCall(SNESVICheckLocalMin_Private(snes,snes->jacobian,vi->phi,X,gnorm,&ismin)); 335 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 336 break; 337 } 338 } 339 /* Update function and solution vectors */ 340 vi->phinorm = gnorm; 341 vi->merit = 0.5*vi->phinorm*vi->phinorm; 342 /* Monitor convergence */ 343 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 344 snes->iter = i+1; 345 snes->norm = vi->phinorm; 346 snes->xnorm = xnorm; 347 snes->ynorm = ynorm; 348 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 349 PetscCall(SNESLogConvergenceHistory(snes,snes->norm,lits)); 350 PetscCall(SNESMonitor(snes,snes->iter,snes->norm)); 351 /* Test for convergence, xnorm = || X || */ 352 if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X,NORM_2,&xnorm)); 353 PetscCall((*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP)); 354 if (snes->reason) break; 355 } 356 if (i == maxits) { 357 PetscCall(PetscInfo(snes,"Maximum number of iterations has been reached: %" PetscInt_FMT "\n",maxits)); 358 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 359 } 360 sdm->ops->computefunction = vi->computeuserfunction; 361 PetscFunctionReturn(0); 362 } 363 364 /* -------------------------------------------------------------------------- */ 365 /* 366 SNESSetUp_VINEWTONSSLS - Sets up the internal data structures for the later use 367 of the SNES nonlinear solver. 368 369 Input Parameter: 370 . snes - the SNES context 371 372 Application Interface Routine: SNESSetUp() 373 374 Notes: 375 For basic use of the SNES solvers, the user need not explicitly call 376 SNESSetUp(), since these actions will automatically occur during 377 the call to SNESSolve(). 378 */ 379 PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes) 380 { 381 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data; 382 383 PetscFunctionBegin; 384 PetscCall(SNESSetUp_VI(snes)); 385 PetscCall(VecDuplicate(snes->vec_sol, &vi->dpsi)); 386 PetscCall(VecDuplicate(snes->vec_sol, &vi->phi)); 387 PetscCall(VecDuplicate(snes->vec_sol, &vi->Da)); 388 PetscCall(VecDuplicate(snes->vec_sol, &vi->Db)); 389 PetscCall(VecDuplicate(snes->vec_sol, &vi->z)); 390 PetscCall(VecDuplicate(snes->vec_sol, &vi->t)); 391 PetscFunctionReturn(0); 392 } 393 /* -------------------------------------------------------------------------- */ 394 PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes) 395 { 396 SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS*) snes->data; 397 398 PetscFunctionBegin; 399 PetscCall(SNESReset_VI(snes)); 400 PetscCall(VecDestroy(&vi->dpsi)); 401 PetscCall(VecDestroy(&vi->phi)); 402 PetscCall(VecDestroy(&vi->Da)); 403 PetscCall(VecDestroy(&vi->Db)); 404 PetscCall(VecDestroy(&vi->z)); 405 PetscCall(VecDestroy(&vi->t)); 406 PetscFunctionReturn(0); 407 } 408 409 /* -------------------------------------------------------------------------- */ 410 /* 411 SNESSetFromOptions_VINEWTONSSLS - Sets various parameters for the SNESVI method. 412 413 Input Parameter: 414 . snes - the SNES context 415 416 Application Interface Routine: SNESSetFromOptions() 417 */ 418 static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(PetscOptionItems *PetscOptionsObject,SNES snes) 419 { 420 PetscFunctionBegin; 421 PetscCall(SNESSetFromOptions_VI(PetscOptionsObject,snes)); 422 PetscOptionsHeadBegin(PetscOptionsObject,"SNES semismooth method options"); 423 PetscOptionsHeadEnd(); 424 PetscFunctionReturn(0); 425 } 426 427 /* -------------------------------------------------------------------------- */ 428 /*MC 429 SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method 430 431 Options Database: 432 + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method 433 - -snes_vi_monitor - prints the number of active constraints at each iteration. 434 435 Level: beginner 436 437 References: 438 + * - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth 439 algorithm for large scale complementarity problems. INFORMS Journal on Computing, 13 (2001). 440 - * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale 441 Applications, Optimization Methods and Software, 21 (2006). 442 443 .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()` 444 445 M*/ 446 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes) 447 { 448 SNES_VINEWTONSSLS *vi; 449 SNESLineSearch linesearch; 450 451 PetscFunctionBegin; 452 snes->ops->reset = SNESReset_VINEWTONSSLS; 453 snes->ops->setup = SNESSetUp_VINEWTONSSLS; 454 snes->ops->solve = SNESSolve_VINEWTONSSLS; 455 snes->ops->destroy = SNESDestroy_VI; 456 snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS; 457 snes->ops->view = NULL; 458 459 snes->usesksp = PETSC_TRUE; 460 snes->usesnpc = PETSC_FALSE; 461 462 PetscCall(SNESGetLineSearch(snes, &linesearch)); 463 if (!((PetscObject)linesearch)->type_name) { 464 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT)); 465 PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0)); 466 } 467 468 snes->alwayscomputesfinalresidual = PETSC_FALSE; 469 470 PetscCall(PetscNewLog(snes,&vi)); 471 snes->data = (void*)vi; 472 473 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI)); 474 PetscCall(PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI)); 475 PetscFunctionReturn(0); 476 } 477