1 #define PETSCSNES_DLL 2 3 #include "../src/snes/impls/vi/viimpl.h" 4 #include "../include/private/kspimpl.h" 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "SNESMonitorVI" 8 PetscErrorCode PETSCSNES_DLLEXPORT SNESMonitorVI(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy) 9 { 10 PetscErrorCode ierr; 11 SNES_VI *vi = (SNES_VI*)snes->data; 12 PetscViewerASCIIMonitor viewer = (PetscViewerASCIIMonitor) dummy; 13 const PetscScalar *x,*xl,*xu,*f; 14 PetscInt i,n,act = 0; 15 PetscReal rnorm,fnorm; 16 17 PetscFunctionBegin; 18 if (!dummy) { 19 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",0,&viewer);CHKERRQ(ierr); 20 } 21 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 22 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 23 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 24 ierr = VecGetArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 25 ierr = VecGetArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 26 27 rnorm = 0.0; 28 for (i=0; i<n; i++) { 29 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) rnorm += f[i]*f[i]; 30 else act++; 31 } 32 ierr = VecRestoreArrayRead(snes->vec_func,&f);CHKERRQ(ierr); 33 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 34 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 35 ierr = VecRestoreArrayRead(snes->vec_sol,&x);CHKERRQ(ierr); 36 ierr = MPI_Allreduce(&rnorm,&fnorm,1,MPIU_REAL,MPI_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 37 fnorm = sqrt(fnorm); 38 ierr = PetscViewerASCIIMonitorPrintf(viewer,"%3D SNES VI Function norm %14.12e Active constraints %D\n",its,fnorm,act);CHKERRQ(ierr); 39 if (!dummy) { 40 ierr = PetscViewerASCIIMonitorDestroy(viewer);CHKERRQ(ierr); 41 } 42 PetscFunctionReturn(0); 43 } 44 45 /* 46 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function, 47 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that 48 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 49 for this trick. One assumes that the probability that W is in the null space of J is very, very small. 50 */ 51 #undef __FUNCT__ 52 #define __FUNCT__ "SNESVICheckLocalMin_Private" 53 PetscErrorCode SNESVICheckLocalMin_Private(SNES snes,Mat A,Vec F,Vec W,PetscReal fnorm,PetscBool *ismin) 54 { 55 PetscReal a1; 56 PetscErrorCode ierr; 57 PetscBool hastranspose; 58 59 PetscFunctionBegin; 60 *ismin = PETSC_FALSE; 61 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 62 if (hastranspose) { 63 /* Compute || J^T F|| */ 64 ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr); 65 ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr); 66 ierr = PetscInfo1(snes,"|| J^T F|| %G near zero implies found a local minimum\n",a1/fnorm);CHKERRQ(ierr); 67 if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE; 68 } else { 69 Vec work; 70 PetscScalar result; 71 PetscReal wnorm; 72 73 ierr = VecSetRandom(W,PETSC_NULL);CHKERRQ(ierr); 74 ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr); 75 ierr = VecDuplicate(W,&work);CHKERRQ(ierr); 76 ierr = MatMult(A,W,work);CHKERRQ(ierr); 77 ierr = VecDot(F,work,&result);CHKERRQ(ierr); 78 ierr = VecDestroy(work);CHKERRQ(ierr); 79 a1 = PetscAbsScalar(result)/(fnorm*wnorm); 80 ierr = PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %G near zero implies found a local minimum\n",a1);CHKERRQ(ierr); 81 if (a1 < 1.e-4) *ismin = PETSC_TRUE; 82 } 83 PetscFunctionReturn(0); 84 } 85 86 /* 87 Checks if J^T(F - J*X) = 0 88 */ 89 #undef __FUNCT__ 90 #define __FUNCT__ "SNESVICheckResidual_Private" 91 PetscErrorCode SNESVICheckResidual_Private(SNES snes,Mat A,Vec F,Vec X,Vec W1,Vec W2) 92 { 93 PetscReal a1,a2; 94 PetscErrorCode ierr; 95 PetscBool hastranspose; 96 97 PetscFunctionBegin; 98 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 99 if (hastranspose) { 100 ierr = MatMult(A,X,W1);CHKERRQ(ierr); 101 ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr); 102 103 /* Compute || J^T W|| */ 104 ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr); 105 ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr); 106 ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr); 107 if (a1 != 0.0) { 108 ierr = PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %G near zero implies inconsistent rhs\n",a2/a1);CHKERRQ(ierr); 109 } 110 } 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 SNESDefaultConverged_VI - Checks the convergence of the semismooth newton algorithm. 116 117 Notes: 118 The convergence criterion currently implemented is 119 merit < abstol 120 merit < rtol*merit_initial 121 */ 122 #undef __FUNCT__ 123 #define __FUNCT__ "SNESDefaultConverged_VI" 124 PetscErrorCode SNESDefaultConverged_VI(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gradnorm,PetscReal merit,SNESConvergedReason *reason,void *dummy) 125 { 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 130 PetscValidPointer(reason,6); 131 132 *reason = SNES_CONVERGED_ITERATING; 133 134 if (!it) { 135 /* set parameter for default relative tolerance convergence test */ 136 snes->ttol = merit*snes->rtol; 137 } 138 if (merit != merit) { 139 ierr = PetscInfo(snes,"Failed to converged, function norm is NaN\n");CHKERRQ(ierr); 140 *reason = SNES_DIVERGED_FNORM_NAN; 141 } else if (merit < snes->abstol) { 142 ierr = PetscInfo2(snes,"Converged due to merit function %G < %G\n",merit,snes->abstol);CHKERRQ(ierr); 143 *reason = SNES_CONVERGED_FNORM_ABS; 144 } else if (snes->nfuncs >= snes->max_funcs) { 145 ierr = PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);CHKERRQ(ierr); 146 *reason = SNES_DIVERGED_FUNCTION_COUNT; 147 } 148 149 if (it && !*reason) { 150 if (merit < snes->ttol) { 151 ierr = PetscInfo2(snes,"Converged due to merit function %G < %G (relative tolerance)\n",merit,snes->ttol);CHKERRQ(ierr); 152 *reason = SNES_CONVERGED_FNORM_RELATIVE; 153 } 154 } 155 PetscFunctionReturn(0); 156 } 157 158 /* 159 SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem. 160 161 Input Parameter: 162 . phi - the semismooth function 163 164 Output Parameter: 165 . merit - the merit function 166 . phinorm - ||phi|| 167 168 Notes: 169 The merit function for the mixed complementarity problem is defined as 170 merit = 0.5*phi^T*phi 171 */ 172 #undef __FUNCT__ 173 #define __FUNCT__ "SNESVIComputeMeritFunction" 174 static PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal* merit,PetscReal* phinorm) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = VecNormBegin(phi,NORM_2,phinorm); 180 ierr = VecNormEnd(phi,NORM_2,phinorm); 181 182 *merit = 0.5*(*phinorm)*(*phinorm); 183 PetscFunctionReturn(0); 184 } 185 186 /* 187 ComputeFischerFunction - Computes the semismooth fischer burmeister function for a mixed complementarity equation. 188 189 Notes: 190 The Fischer-Burmeister function is defined as 191 ff(a,b) := sqrt(a*a + b*b) - a - b 192 and is used reformulate a complementarity equation as a semismooth equation. 193 */ 194 195 #undef __FUNCT__ 196 #define __FUNCT__ "ComputeFischerFunction" 197 static PetscErrorCode ComputeFischerFunction(PetscScalar a, PetscScalar b, PetscScalar* ff) 198 { 199 PetscFunctionBegin; 200 *ff = sqrt(a*a + b*b) - a - b; 201 PetscFunctionReturn(0); 202 } 203 204 /* 205 SNESVIComputeFunction - Reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear equations in semismooth form. 206 207 Input Parameters: 208 . snes - the SNES context 209 . x - current iterate 210 . functx - user defined function context 211 212 Output Parameters: 213 . phi - Semismooth function 214 215 Notes: 216 The result of this function is done by cases: 217 + l[i] == -infinity, u[i] == infinity -- phi[i] = -f[i] 218 . l[i] == -infinity, u[i] finite -- phi[i] = ss(u[i]-x[i], -f[i]) 219 . l[i] finite, u[i] == infinity -- phi[i] = ss(x[i]-l[i], f[i]) 220 . l[i] finite < u[i] finite -- phi[i] = phi(x[i]-l[i], ss(u[i]-x[i], -f[u])) 221 - otherwise l[i] == u[i] -- phi[i] = l[i] - x[i] 222 ss is the semismoothing function used to reformulate the nonlinear equation in complementarity 223 form to semismooth form 224 225 */ 226 #undef __FUNCT__ 227 #define __FUNCT__ "SNESVIComputeFunction" 228 static PetscErrorCode SNESVIComputeFunction(SNES snes,Vec X,Vec phi,void* functx) 229 { 230 PetscErrorCode ierr; 231 SNES_VI *vi = (SNES_VI*)snes->data; 232 Vec Xl = vi->xl,Xu = vi->xu,F = snes->vec_func; 233 PetscScalar *phi_arr,*x_arr,*f_arr,*l,*u,t; 234 PetscInt i,nlocal; 235 236 PetscFunctionBegin; 237 ierr = (*vi->computeuserfunction)(snes,X,F,functx);CHKERRQ(ierr); 238 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 239 240 ierr = VecGetArray(X,&x_arr);CHKERRQ(ierr); 241 ierr = VecGetArray(F,&f_arr);CHKERRQ(ierr); 242 ierr = VecGetArray(Xl,&l);CHKERRQ(ierr); 243 ierr = VecGetArray(Xu,&u);CHKERRQ(ierr); 244 ierr = VecGetArray(phi,&phi_arr);CHKERRQ(ierr); 245 246 for (i=0;i < nlocal;i++) { 247 if ((l[i] <= PETSC_VI_NINF) && (u[i] >= PETSC_VI_INF)) { 248 phi_arr[i] = -f_arr[i]; 249 } 250 else if (l[i] <= PETSC_VI_NINF) { 251 t = u[i] - x_arr[i]; 252 ierr = ComputeFischerFunction(t,-f_arr[i],&phi_arr[i]);CHKERRQ(ierr); 253 phi_arr[i] = -phi_arr[i]; 254 } 255 else if (u[i] >= PETSC_VI_INF) { 256 t = x_arr[i] - l[i]; 257 ierr = ComputeFischerFunction(t,f_arr[i],&phi_arr[i]);CHKERRQ(ierr); 258 } 259 else if (l[i] == u[i]) { 260 phi_arr[i] = l[i] - x_arr[i]; 261 } 262 else { 263 t = u[i] - x_arr[i]; 264 ierr = ComputeFischerFunction(t,-f_arr[i],&phi_arr[i]); 265 t = x_arr[i] - l[i]; 266 ierr = ComputeFischerFunction(t,phi_arr[i],&phi_arr[i]); 267 } 268 } 269 270 ierr = VecRestoreArray(X,&x_arr);CHKERRQ(ierr); 271 ierr = VecRestoreArray(F,&f_arr);CHKERRQ(ierr); 272 ierr = VecRestoreArray(Xl,&l);CHKERRQ(ierr); 273 ierr = VecRestoreArray(Xu,&u);CHKERRQ(ierr); 274 ierr = VecRestoreArray(phi,&phi_arr);CHKERRQ(ierr); 275 276 PetscFunctionReturn(0); 277 } 278 279 /* 280 SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the 281 the semismooth jacobian. 282 */ 283 #undef __FUNCT__ 284 #define __FUNCT__ "SNESVIComputeBsubdifferentialVectors" 285 PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes,Vec X,Vec F,Mat jac,Vec Da,Vec Db) 286 { 287 PetscErrorCode ierr; 288 SNES_VI *vi = (SNES_VI*)snes->data; 289 PetscScalar *l,*u,*x,*f,*da,*db,*z,*t,t1,t2,ci,di,ei; 290 PetscInt i,nlocal; 291 292 PetscFunctionBegin; 293 294 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 295 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 296 ierr = VecGetArray(vi->xl,&l);CHKERRQ(ierr); 297 ierr = VecGetArray(vi->xu,&u);CHKERRQ(ierr); 298 ierr = VecGetArray(Da,&da);CHKERRQ(ierr); 299 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 300 ierr = VecGetArray(vi->z,&z);CHKERRQ(ierr); 301 302 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 303 /* Set the elements of the vector z: 304 z[i] = 1 if (x[i] - l[i],f[i]) = (0,0) or (u[i] - x[i],f[i]) = (0,0) 305 else z[i] = 0 306 */ 307 for(i=0;i < nlocal;i++) { 308 da[i] = db[i] = z[i] = 0; 309 if(PetscAbsScalar(f[i]) <= vi->const_tol) { 310 if ((l[i] > PETSC_VI_NINF) && (PetscAbsScalar(x[i]-l[i]) <= vi->const_tol)) { 311 da[i] = 1; 312 z[i] = 1; 313 } 314 if ((u[i] < PETSC_VI_INF) && (PetscAbsScalar(u[i]-x[i]) <= vi->const_tol)) { 315 db[i] = 1; 316 z[i] = 1; 317 } 318 } 319 } 320 ierr = VecRestoreArray(vi->z,&z);CHKERRQ(ierr); 321 ierr = MatMult(jac,vi->z,vi->t);CHKERRQ(ierr); 322 ierr = VecGetArray(vi->t,&t);CHKERRQ(ierr); 323 /* Compute the elements of the diagonal perturbation vector Da and row scaling vector Db */ 324 for(i=0;i< nlocal;i++) { 325 /* Free variables */ 326 if ((l[i] <= PETSC_VI_NINF) && (u[i] >= PETSC_VI_INF)) { 327 da[i] = 0; db[i] = -1; 328 } 329 /* lower bounded variables */ 330 else if (u[i] >= PETSC_VI_INF) { 331 if (da[i] >= 1) { 332 t2 = PetscScalarNorm(1,t[i]); 333 da[i] = 1/t2 - 1; 334 db[i] = t[i]/t2 - 1; 335 } else { 336 t1 = x[i] - l[i]; 337 t2 = PetscScalarNorm(t1,f[i]); 338 da[i] = t1/t2 - 1; 339 db[i] = f[i]/t2 - 1; 340 } 341 } 342 /* upper bounded variables */ 343 else if (l[i] <= PETSC_VI_NINF) { 344 if (db[i] >= 1) { 345 t2 = PetscScalarNorm(1,t[i]); 346 da[i] = -1/t2 -1; 347 db[i] = -t[i]/t2 - 1; 348 } 349 else { 350 t1 = u[i] - x[i]; 351 t2 = PetscScalarNorm(t1,f[i]); 352 da[i] = t1/t2 - 1; 353 db[i] = -f[i]/t2 - 1; 354 } 355 } 356 /* Fixed variables */ 357 else if (l[i] == u[i]) { 358 da[i] = -1; 359 db[i] = 0; 360 } 361 /* Box constrained variables */ 362 else { 363 if (db[i] >= 1) { 364 t2 = PetscScalarNorm(1,t[i]); 365 ci = 1/t2 + 1; 366 di = t[i]/t2 + 1; 367 } 368 else { 369 t1 = x[i] - u[i]; 370 t2 = PetscScalarNorm(t1,f[i]); 371 ci = t1/t2 + 1; 372 di = f[i]/t2 + 1; 373 } 374 375 if (da[i] >= 1) { 376 t1 = ci + di*t[i]; 377 t2 = PetscScalarNorm(1,t1); 378 t1 = t1/t2 - 1; 379 t2 = 1/t2 - 1; 380 } 381 else { 382 ierr = ComputeFischerFunction(u[i]-x[i],-f[i],&ei);CHKERRQ(ierr); 383 t2 = PetscScalarNorm(x[i]-l[i],ei); 384 t1 = ei/t2 - 1; 385 t2 = (x[i] - l[i])/t2 - 1; 386 } 387 388 da[i] = t2 + t1*ci; 389 db[i] = t1*di; 390 } 391 } 392 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 393 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 394 ierr = VecRestoreArray(vi->xl,&l);CHKERRQ(ierr); 395 ierr = VecRestoreArray(vi->xu,&u);CHKERRQ(ierr); 396 ierr = VecRestoreArray(Da,&da);CHKERRQ(ierr); 397 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 398 ierr = VecRestoreArray(vi->t,&t);CHKERRQ(ierr); 399 400 PetscFunctionReturn(0); 401 } 402 403 /* 404 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. 405 406 Input Parameters: 407 . Da - Diagonal shift vector for the semismooth jacobian. 408 . Db - Row scaling vector for the semismooth jacobian. 409 410 Output Parameters: 411 . jac - semismooth jacobian 412 . jac_pre - optional preconditioning matrix 413 414 Notes: 415 The semismooth jacobian matrix is given by 416 jac = Da + Db*jacfun 417 where Db is the row scaling matrix stored as a vector, 418 Da is the diagonal perturbation matrix stored as a vector 419 and jacfun is the jacobian of the original nonlinear function. 420 */ 421 #undef __FUNCT__ 422 #define __FUNCT__ "SNESVIComputeJacobian" 423 PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre,Vec Da, Vec Db) 424 { 425 PetscErrorCode ierr; 426 427 /* Do row scaling and add diagonal perturbation */ 428 ierr = MatDiagonalScale(jac,Db,PETSC_NULL);CHKERRQ(ierr); 429 ierr = MatDiagonalSet(jac,Da,ADD_VALUES);CHKERRQ(ierr); 430 if (jac != jac_pre) { /* If jac and jac_pre are different */ 431 ierr = MatDiagonalScale(jac_pre,Db,PETSC_NULL); 432 ierr = MatDiagonalSet(jac_pre,Da,ADD_VALUES);CHKERRQ(ierr); 433 } 434 PetscFunctionReturn(0); 435 } 436 437 /* 438 SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi. 439 440 Input Parameters: 441 phi - semismooth function. 442 H - semismooth jacobian 443 444 Output Parameters: 445 dpsi - merit function gradient 446 447 Notes: 448 The merit function gradient is computed as follows 449 dpsi = H^T*phi 450 */ 451 #undef __FUNCT__ 452 #define __FUNCT__ "SNESVIComputeMeritFunctionGradient" 453 PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi) 454 { 455 PetscErrorCode ierr; 456 457 PetscFunctionBegin; 458 ierr = MatMultTranspose(H,phi,dpsi); 459 PetscFunctionReturn(0); 460 } 461 462 /* 463 SNESVIAdjustInitialGuess - Readjusts the initial guess to the SNES solver supplied by the user so that the initial guess lies inside the feasible region . 464 465 Input Parameters: 466 . lb - lower bound. 467 . ub - upper bound. 468 469 Output Parameters: 470 . X - the readjusted initial guess. 471 472 Notes: 473 The readjusted initial guess X[i] = max(max(min(l[i],X[i]),min(X[i],u[i])),min(l[i],u[i])) 474 */ 475 #undef __FUNCT__ 476 #define __FUNCT__ "SNESVIAdjustInitialGuess" 477 PetscErrorCode SNESVIAdjustInitialGuess(Vec X, Vec lb, Vec ub) 478 { 479 PetscErrorCode ierr; 480 PetscInt i,nlocal; 481 PetscScalar *x,*l,*u; 482 483 PetscFunctionBegin; 484 485 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 486 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 487 ierr = VecGetArray(lb,&l);CHKERRQ(ierr); 488 ierr = VecGetArray(ub,&u);CHKERRQ(ierr); 489 490 for(i = 0; i < nlocal; i++) { 491 x[i] = PetscMax(PetscMax(PetscMin(x[i],l[i]),PetscMin(x[i],u[i])),PetscMin(l[i],u[i])); 492 } 493 494 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 495 ierr = VecRestoreArray(lb,&l);CHKERRQ(ierr); 496 ierr = VecRestoreArray(ub,&u);CHKERRQ(ierr); 497 498 PetscFunctionReturn(0); 499 } 500 501 /* -------------------------------------------------------------------- 502 503 This file implements a semismooth truncated Newton method with a line search, 504 for solving a system of nonlinear equations in complementarity form, using the KSP, Vec, 505 and Mat interfaces for linear solvers, vectors, and matrices, 506 respectively. 507 508 The following basic routines are required for each nonlinear solver: 509 SNESCreate_XXX() - Creates a nonlinear solver context 510 SNESSetFromOptions_XXX() - Sets runtime options 511 SNESSolve_XXX() - Solves the nonlinear system 512 SNESDestroy_XXX() - Destroys the nonlinear solver context 513 The suffix "_XXX" denotes a particular implementation, in this case 514 we use _VI (e.g., SNESCreate_VI, SNESSolve_VI) for solving 515 systems of nonlinear equations with a line search (LS) method. 516 These routines are actually called via the common user interface 517 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 518 SNESDestroy(), so the application code interface remains identical 519 for all nonlinear solvers. 520 521 Another key routine is: 522 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 523 by setting data structures and options. The interface routine SNESSetUp() 524 is not usually called directly by the user, but instead is called by 525 SNESSolve() if necessary. 526 527 Additional basic routines are: 528 SNESView_XXX() - Prints details of runtime options that 529 have actually been used. 530 These are called by application codes via the interface routines 531 SNESView(). 532 533 The various types of solvers (preconditioners, Krylov subspace methods, 534 nonlinear solvers, timesteppers) are all organized similarly, so the 535 above description applies to these categories also. 536 537 -------------------------------------------------------------------- */ 538 /* 539 SNESSolveVI_SS - Solves the complementarity problem with a semismooth Newton 540 method using a line search. 541 542 Input Parameters: 543 . snes - the SNES context 544 545 Output Parameter: 546 . outits - number of iterations until termination 547 548 Application Interface Routine: SNESSolve() 549 550 Notes: 551 This implements essentially a semismooth Newton method with a 552 line search. By default a cubic backtracking line search 553 is employed, as described in the text "Numerical Methods for 554 Unconstrained Optimization and Nonlinear Equations" by Dennis 555 and Schnabel. 556 */ 557 #undef __FUNCT__ 558 #define __FUNCT__ "SNESSolveVI_SS" 559 PetscErrorCode SNESSolveVI_SS(SNES snes) 560 { 561 SNES_VI *vi = (SNES_VI*)snes->data; 562 PetscErrorCode ierr; 563 PetscInt maxits,i,lits; 564 PetscBool lssucceed; 565 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 566 PetscReal gnorm,xnorm=0,ynorm; 567 Vec Y,X,F,G,W; 568 KSPConvergedReason kspreason; 569 570 PetscFunctionBegin; 571 snes->numFailures = 0; 572 snes->numLinearSolveFailures = 0; 573 snes->reason = SNES_CONVERGED_ITERATING; 574 575 maxits = snes->max_its; /* maximum number of iterations */ 576 X = snes->vec_sol; /* solution vector */ 577 F = snes->vec_func; /* residual vector */ 578 Y = snes->work[0]; /* work vectors */ 579 G = snes->work[1]; 580 W = snes->work[2]; 581 582 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 583 snes->iter = 0; 584 snes->norm = 0.0; 585 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 586 587 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 588 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 589 if (snes->domainerror) { 590 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 591 PetscFunctionReturn(0); 592 } 593 /* Compute Merit function */ 594 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 595 596 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 597 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 598 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 599 600 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 601 snes->norm = vi->phinorm; 602 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 603 SNESLogConvHistory(snes,vi->phinorm,0); 604 SNESMonitor(snes,0,vi->phinorm); 605 606 /* set parameter for default relative tolerance convergence test */ 607 snes->ttol = vi->phinorm*snes->rtol; 608 /* test convergence */ 609 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 610 if (snes->reason) PetscFunctionReturn(0); 611 612 for (i=0; i<maxits; i++) { 613 614 /* Call general purpose update function */ 615 if (snes->ops->update) { 616 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 617 } 618 619 /* Solve J Y = Phi, where J is the semismooth jacobian */ 620 /* Get the nonlinear function jacobian */ 621 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 622 /* Get the diagonal shift and row scaling vectors */ 623 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 624 /* Compute the semismooth jacobian */ 625 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 626 /* Compute the merit function gradient */ 627 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 628 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 629 ierr = SNES_KSPSolve(snes,snes->ksp,vi->phi,Y);CHKERRQ(ierr); 630 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 631 632 if (kspreason < 0) { 633 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 634 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 635 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 636 break; 637 } 638 } 639 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 640 snes->linear_its += lits; 641 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 642 /* 643 if (vi->precheckstep) { 644 PetscBool changed_y = PETSC_FALSE; 645 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 646 } 647 648 if (PetscLogPrintInfo){ 649 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 650 } 651 */ 652 /* Compute a (scaled) negative update in the line search routine: 653 Y <- X - lambda*Y 654 and evaluate G = function(Y) (depends on the line search). 655 */ 656 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 657 ynorm = 1; gnorm = vi->phinorm; 658 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 659 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 660 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 661 if (snes->domainerror) { 662 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 663 PetscFunctionReturn(0); 664 } 665 if (!lssucceed) { 666 if (++snes->numFailures >= snes->maxFailures) { 667 PetscBool ismin; 668 snes->reason = SNES_DIVERGED_LINE_SEARCH; 669 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 670 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 671 break; 672 } 673 } 674 /* Update function and solution vectors */ 675 vi->phinorm = gnorm; 676 vi->merit = 0.5*vi->phinorm*vi->phinorm; 677 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 678 ierr = VecCopy(W,X);CHKERRQ(ierr); 679 /* Monitor convergence */ 680 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 681 snes->iter = i+1; 682 snes->norm = vi->phinorm; 683 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 684 SNESLogConvHistory(snes,snes->norm,lits); 685 SNESMonitor(snes,snes->iter,snes->norm); 686 /* Test for convergence, xnorm = || X || */ 687 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 688 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 689 if (snes->reason) break; 690 } 691 if (i == maxits) { 692 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 693 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 694 } 695 PetscFunctionReturn(0); 696 } 697 698 #undef __FUNCT__ 699 #define __FUNCT__ "SNESVICreateIndexSets_AS" 700 PetscErrorCode SNESVICreateIndexSets_AS(SNES snes,Vec Db,PetscReal thresh,IS* ISact,IS* ISinact) 701 { 702 PetscErrorCode ierr; 703 PetscInt i,nlocal,ilow,ihigh,nloc_isact=0,nloc_isinact=0; 704 PetscInt *idx_act,*idx_inact,i1=0,i2=0; 705 PetscScalar *db; 706 707 PetscFunctionBegin; 708 709 ierr = VecGetLocalSize(Db,&nlocal);CHKERRQ(ierr); 710 ierr = VecGetOwnershipRange(Db,&ilow,&ihigh);CHKERRQ(ierr); 711 ierr = VecGetArray(Db,&db);CHKERRQ(ierr); 712 /* Compute the sizes of the active and inactive sets */ 713 for (i=0; i < nlocal;i++) { 714 if (PetscAbsScalar(db[i]) <= thresh) nloc_isact++; 715 else nloc_isinact++; 716 } 717 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 718 ierr = PetscMalloc(nloc_isinact*sizeof(PetscInt),&idx_inact);CHKERRQ(ierr); 719 720 /* Creating the indexing arrays */ 721 for(i=0; i < nlocal; i++) { 722 if (PetscAbsScalar(db[i]) <= thresh) idx_act[i1++] = ilow+i; 723 else idx_inact[i2++] = ilow+i; 724 } 725 726 /* Create the index sets */ 727 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_COPY_VALUES,ISact);CHKERRQ(ierr); 728 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isinact,idx_inact,PETSC_COPY_VALUES,ISinact);CHKERRQ(ierr); 729 730 ierr = VecRestoreArray(Db,&db);CHKERRQ(ierr); 731 ierr = PetscFree(idx_act);CHKERRQ(ierr); 732 ierr = PetscFree(idx_inact);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 #undef __FUNCT__ 737 #define __FUNCT__ "SNESVICreateIndexSets_RS" 738 PetscErrorCode SNESVICreateIndexSets_RS(SNES snes,Vec X,Vec Xl, Vec Xu,IS* ISact,IS* ISinact) 739 { 740 PetscErrorCode ierr; 741 SNES_VI *vi = (SNES_VI*)snes->data; 742 PetscInt i,nlocal,ilow,ihigh,nloc_isact=0,nloc_isinact=0; 743 PetscInt *idx_act,*idx_inact,i1=0,i2=0; 744 PetscScalar *x,*l,*u,*f; 745 Vec F = snes->vec_func; 746 747 PetscFunctionBegin; 748 749 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 750 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 751 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 752 ierr = VecGetArray(Xl,&l);CHKERRQ(ierr); 753 ierr = VecGetArray(Xu,&u);CHKERRQ(ierr); 754 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 755 /* Compute the sizes of the active and inactive sets */ 756 for (i=0; i < nlocal;i++) { 757 if (PetscAbsScalar(x[i] - l[i]) <= vi->const_tol && f[i] > vi->const_tol) nloc_isact++; 758 else if (PetscAbsScalar(u[i] - x[i]) <= vi->const_tol && f[i] < vi->const_tol) nloc_isact++; 759 else nloc_isinact++; 760 } 761 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 762 ierr = PetscMalloc(nloc_isinact*sizeof(PetscInt),&idx_inact);CHKERRQ(ierr); 763 764 /* Creating the indexing arrays */ 765 for(i=0; i < nlocal; i++) { 766 if (PetscAbsScalar(x[i] - l[i]) <= vi->const_tol && f[i] > vi->const_tol) idx_act[i1++] = ilow+i; 767 else if (PetscAbsScalar(u[i] - x[i]) <= vi->const_tol && f[i] < vi->const_tol) idx_act[i1++] = ilow+i; 768 else idx_inact[i2++] = ilow+i; 769 } 770 771 /* Create the index sets */ 772 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_COPY_VALUES,ISact);CHKERRQ(ierr); 773 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isinact,idx_inact,PETSC_COPY_VALUES,ISinact);CHKERRQ(ierr); 774 775 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 776 ierr = VecRestoreArray(Xl,&l);CHKERRQ(ierr); 777 ierr = VecRestoreArray(Xu,&u);CHKERRQ(ierr); 778 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 779 ierr = PetscFree(idx_act);CHKERRQ(ierr); 780 ierr = PetscFree(idx_inact);CHKERRQ(ierr); 781 PetscFunctionReturn(0); 782 } 783 784 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 785 #undef __FUNCT__ 786 #define __FUNCT__ "SNESVICreateVectors_AS" 787 PetscErrorCode SNESVICreateVectors_AS(SNES snes,PetscInt n,Vec* newv) 788 { 789 PetscErrorCode ierr; 790 Vec v; 791 792 PetscFunctionBegin; 793 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 794 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 795 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 796 *newv = v; 797 798 PetscFunctionReturn(0); 799 } 800 801 /* Resets the snes PC and KSP when the active set sizes change */ 802 #undef __FUNCT__ 803 #define __FUNCT__ "SNESVIResetPCandKSP" 804 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 805 { 806 PetscErrorCode ierr; 807 KSP kspnew,snesksp; 808 PC pcnew; 809 const MatSolverPackage stype; 810 811 PetscFunctionBegin; 812 /* The active and inactive set sizes have changed so need to create a new snes->ksp object */ 813 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 814 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 815 /* Copy over snes->ksp info */ 816 kspnew->pc_side = snesksp->pc_side; 817 kspnew->rtol = snesksp->rtol; 818 kspnew->abstol = snesksp->abstol; 819 kspnew->max_it = snesksp->max_it; 820 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 821 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 822 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 823 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 824 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 825 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 826 ierr = KSPDestroy(snesksp);CHKERRQ(ierr); 827 snes->ksp = kspnew; 828 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 829 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 830 PetscFunctionReturn(0); 831 } 832 /* Variational Inequality solver using active set method */ 833 #undef __FUNCT__ 834 #define __FUNCT__ "SNESSolveVI_AS" 835 PetscErrorCode SNESSolveVI_AS(SNES snes) 836 { 837 SNES_VI *vi = (SNES_VI*)snes->data; 838 PetscErrorCode ierr; 839 PetscInt maxits,i,lits,Nis_act=0; 840 PetscBool lssucceed; 841 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 842 PetscReal gnorm,xnorm=0,ynorm; 843 Vec Y,X,F,G,W; 844 KSPConvergedReason kspreason; 845 846 PetscFunctionBegin; 847 snes->numFailures = 0; 848 snes->numLinearSolveFailures = 0; 849 snes->reason = SNES_CONVERGED_ITERATING; 850 851 maxits = snes->max_its; /* maximum number of iterations */ 852 X = snes->vec_sol; /* solution vector */ 853 F = snes->vec_func; /* residual vector */ 854 Y = snes->work[0]; /* work vectors */ 855 G = snes->work[1]; 856 W = snes->work[2]; 857 858 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 859 snes->iter = 0; 860 snes->norm = 0.0; 861 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 862 863 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 864 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 865 if (snes->domainerror) { 866 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 867 PetscFunctionReturn(0); 868 } 869 /* Compute Merit function */ 870 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 871 872 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 873 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 874 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 875 876 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 877 snes->norm = vi->phinorm; 878 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 879 SNESLogConvHistory(snes,vi->phinorm,0); 880 SNESMonitor(snes,0,vi->phinorm); 881 882 /* set parameter for default relative tolerance convergence test */ 883 snes->ttol = vi->phinorm*snes->rtol; 884 /* test convergence */ 885 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 886 if (snes->reason) PetscFunctionReturn(0); 887 888 for (i=0; i<maxits; i++) { 889 890 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 891 PetscReal thresh,J_norm1; 892 VecScatter scat_act,scat_inact; 893 PetscInt nis_act,nis_inact,Nis_act_prev; 894 Vec Da_act,Da_inact,Db_inact; 895 Vec Y_act,Y_inact,phi_act,phi_inact; 896 Mat jac_inact_inact,jac_inact_act,prejac_inact_inact; 897 898 /* Call general purpose update function */ 899 if (snes->ops->update) { 900 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 901 } 902 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 903 /* Compute the threshold value for creating active and inactive sets */ 904 ierr = MatNorm(snes->jacobian,NORM_1,&J_norm1);CHKERRQ(ierr); 905 thresh = PetscMin(vi->merit,1e-2)/(1+J_norm1); 906 907 /* Compute B-subdifferential vectors Da and Db */ 908 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 909 910 /* Create active and inactive index sets */ 911 ierr = SNESVICreateIndexSets_AS(snes,vi->Db,thresh,&IS_act,&IS_inact);CHKERRQ(ierr); 912 913 Nis_act_prev = Nis_act; 914 /* Get sizes of active and inactive sets */ 915 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 916 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 917 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 918 919 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 920 921 /* Create active and inactive set vectors */ 922 ierr = SNESVICreateVectors_AS(snes,nis_act,&Da_act);CHKERRQ(ierr); 923 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Da_inact);CHKERRQ(ierr); 924 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Db_inact);CHKERRQ(ierr); 925 ierr = SNESVICreateVectors_AS(snes,nis_act,&phi_act);CHKERRQ(ierr); 926 ierr = SNESVICreateVectors_AS(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 927 ierr = SNESVICreateVectors_AS(snes,nis_act,&Y_act);CHKERRQ(ierr); 928 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 929 930 /* Create inactive set submatrices */ 931 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_act,MAT_INITIAL_MATRIX,&jac_inact_act);CHKERRQ(ierr); 932 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 933 934 /* Create scatter contexts */ 935 ierr = VecScatterCreate(vi->Da,IS_act,Da_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 936 ierr = VecScatterCreate(vi->Da,IS_inact,Da_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 937 938 /* Do a vec scatter to active and inactive set vectors */ 939 ierr = VecScatterBegin(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 940 ierr = VecScatterEnd(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 941 942 ierr = VecScatterBegin(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 943 ierr = VecScatterEnd(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 944 945 ierr = VecScatterBegin(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 946 ierr = VecScatterEnd(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 947 948 ierr = VecScatterBegin(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 949 ierr = VecScatterEnd(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 950 951 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 952 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 953 954 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 955 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 956 957 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 958 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 959 960 /* Active set direction */ 961 ierr = VecPointwiseDivide(Y_act,phi_act,Da_act);CHKERRQ(ierr); 962 /* inactive set jacobian and preconditioner */ 963 ierr = VecPointwiseDivide(Da_inact,Da_inact,Db_inact);CHKERRQ(ierr); 964 ierr = MatDiagonalSet(jac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 965 if (snes->jacobian != snes->jacobian_pre) { 966 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 967 ierr = MatDiagonalSet(prejac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 968 } else prejac_inact_inact = jac_inact_inact; 969 970 /* right hand side */ 971 ierr = VecPointwiseDivide(phi_inact,phi_inact,Db_inact);CHKERRQ(ierr); 972 ierr = MatMult(jac_inact_act,Y_act,Db_inact);CHKERRQ(ierr); 973 ierr = VecAXPY(phi_inact,-1.0,Db_inact);CHKERRQ(ierr); 974 975 if ((i != 0) && (Nis_act != Nis_act_prev)) { 976 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 977 } 978 979 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 980 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 981 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 982 if (kspreason < 0) { 983 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 984 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 985 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 986 break; 987 } 988 } 989 990 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 991 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 992 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 993 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 994 995 ierr = VecDestroy(Da_act);CHKERRQ(ierr); 996 ierr = VecDestroy(Da_inact);CHKERRQ(ierr); 997 ierr = VecDestroy(Db_inact);CHKERRQ(ierr); 998 ierr = VecDestroy(phi_act);CHKERRQ(ierr); 999 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 1000 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 1001 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 1002 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1003 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1004 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1005 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1006 ierr = MatDestroy(jac_inact_act);CHKERRQ(ierr); 1007 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1008 if (snes->jacobian != snes->jacobian_pre) { 1009 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1010 } 1011 1012 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1013 snes->linear_its += lits; 1014 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1015 /* 1016 if (vi->precheckstep) { 1017 PetscBool changed_y = PETSC_FALSE; 1018 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1019 } 1020 1021 if (PetscLogPrintInfo){ 1022 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1023 } 1024 */ 1025 /* Compute a (scaled) negative update in the line search routine: 1026 Y <- X - lambda*Y 1027 and evaluate G = function(Y) (depends on the line search). 1028 */ 1029 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1030 ynorm = 1; gnorm = vi->phinorm; 1031 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 1032 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 1033 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1034 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1035 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1036 if (snes->domainerror) { 1037 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1038 PetscFunctionReturn(0); 1039 } 1040 if (!lssucceed) { 1041 if (++snes->numFailures >= snes->maxFailures) { 1042 PetscBool ismin; 1043 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1044 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1045 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1046 break; 1047 } 1048 } 1049 /* Update function and solution vectors */ 1050 vi->phinorm = gnorm; 1051 vi->merit = 0.5*vi->phinorm*vi->phinorm; 1052 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 1053 ierr = VecCopy(W,X);CHKERRQ(ierr); 1054 /* Monitor convergence */ 1055 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1056 snes->iter = i+1; 1057 snes->norm = vi->phinorm; 1058 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1059 SNESLogConvHistory(snes,snes->norm,lits); 1060 SNESMonitor(snes,snes->iter,snes->norm); 1061 /* Test for convergence, xnorm = || X || */ 1062 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1063 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1064 if (snes->reason) break; 1065 } 1066 if (i == maxits) { 1067 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1068 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 1073 /* Variational Inequality solver using active set method */ 1074 #undef __FUNCT__ 1075 #define __FUNCT__ "SNESSolveVI_RS" 1076 PetscErrorCode SNESSolveVI_RS(SNES snes) 1077 { 1078 SNES_VI *vi = (SNES_VI*)snes->data; 1079 PetscErrorCode ierr; 1080 PetscInt maxits,i,lits,Nis_act=0; 1081 PetscBool lssucceed; 1082 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1083 PetscReal gnorm,xnorm=0,ynorm; 1084 Vec Y,X,F,G,W; 1085 KSPConvergedReason kspreason; 1086 1087 PetscFunctionBegin; 1088 snes->numFailures = 0; 1089 snes->numLinearSolveFailures = 0; 1090 snes->reason = SNES_CONVERGED_ITERATING; 1091 1092 maxits = snes->max_its; /* maximum number of iterations */ 1093 X = snes->vec_sol; /* solution vector */ 1094 F = snes->vec_func; /* residual vector */ 1095 Y = snes->work[0]; /* work vectors */ 1096 G = snes->work[1]; 1097 W = snes->work[2]; 1098 1099 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1100 snes->iter = 0; 1101 snes->norm = 0.0; 1102 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1103 1104 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 1105 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 1106 if (snes->domainerror) { 1107 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1108 PetscFunctionReturn(0); 1109 } 1110 /* Compute Merit function */ 1111 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 1112 1113 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1114 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1115 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1116 1117 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1118 snes->norm = vi->phinorm; 1119 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1120 SNESLogConvHistory(snes,vi->phinorm,0); 1121 SNESMonitor(snes,0,vi->phinorm); 1122 1123 /* set parameter for default relative tolerance convergence test */ 1124 snes->ttol = vi->phinorm*snes->rtol; 1125 /* test convergence */ 1126 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1127 if (snes->reason) PetscFunctionReturn(0); 1128 1129 for (i=0; i<maxits; i++) { 1130 1131 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1132 VecScatter scat_act,scat_inact; 1133 PetscInt nis_act,nis_inact,Nis_act_prev; 1134 Vec Y_act,Y_inact,phi_inact; 1135 Mat jac_inact_inact,prejac_inact_inact; 1136 1137 /* Call general purpose update function */ 1138 if (snes->ops->update) { 1139 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1140 } 1141 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1142 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 1143 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 1144 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 1145 /* Create active and inactive index sets */ 1146 ierr = SNESVICreateIndexSets_RS(snes,X,vi->xl,vi->xu,&IS_act,&IS_inact);CHKERRQ(ierr); 1147 1148 Nis_act_prev = Nis_act; 1149 /* Get sizes of active and inactive sets */ 1150 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1151 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1152 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 1153 1154 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 1155 1156 /* Create active and inactive set vectors */ 1157 ierr = SNESVICreateVectors_AS(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 1158 ierr = SNESVICreateVectors_AS(snes,nis_act,&Y_act);CHKERRQ(ierr); 1159 ierr = SNESVICreateVectors_AS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1160 1161 /* Create inactive set submatrices */ 1162 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1163 1164 /* Create scatter contexts */ 1165 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1166 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1167 1168 /* Do a vec scatter to active and inactive set vectors */ 1169 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1170 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1171 1172 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1173 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1174 1175 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1176 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1177 1178 /* Active set direction = 0*/ 1179 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1180 if (snes->jacobian != snes->jacobian_pre) { 1181 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1182 } else prejac_inact_inact = jac_inact_inact; 1183 1184 if ((i != 0) && (Nis_act != Nis_act_prev)) { 1185 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1186 } 1187 1188 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1189 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 1190 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1191 if (kspreason < 0) { 1192 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1193 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1194 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1195 break; 1196 } 1197 } 1198 1199 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1200 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1201 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1202 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1203 1204 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 1205 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 1206 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 1207 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1208 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1209 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1210 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1211 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1212 if (snes->jacobian != snes->jacobian_pre) { 1213 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1214 } 1215 1216 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1217 snes->linear_its += lits; 1218 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1219 /* 1220 if (vi->precheckstep) { 1221 PetscBool changed_y = PETSC_FALSE; 1222 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1223 } 1224 1225 if (PetscLogPrintInfo){ 1226 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1227 } 1228 */ 1229 /* Compute a (scaled) negative update in the line search routine: 1230 Y <- X - lambda*Y 1231 and evaluate G = function(Y) (depends on the line search). 1232 */ 1233 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1234 ynorm = 1; gnorm = vi->phinorm; 1235 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1236 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1237 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1238 if (snes->domainerror) { 1239 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1240 PetscFunctionReturn(0); 1241 } 1242 if (!lssucceed) { 1243 if (++snes->numFailures >= snes->maxFailures) { 1244 PetscBool ismin; 1245 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1246 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1247 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1248 break; 1249 } 1250 } 1251 /* Update function and solution vectors */ 1252 vi->phinorm = gnorm; 1253 vi->merit = 0.5*vi->phinorm*vi->phinorm; 1254 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 1255 ierr = VecCopy(W,X);CHKERRQ(ierr); 1256 /* Monitor convergence */ 1257 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1258 snes->iter = i+1; 1259 snes->norm = vi->phinorm; 1260 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1261 SNESLogConvHistory(snes,snes->norm,lits); 1262 SNESMonitor(snes,snes->iter,snes->norm); 1263 /* Test for convergence, xnorm = || X || */ 1264 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1265 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1266 if (snes->reason) break; 1267 } 1268 if (i == maxits) { 1269 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1270 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1271 } 1272 PetscFunctionReturn(0); 1273 } 1274 1275 /* -------------------------------------------------------------------------- */ 1276 /* 1277 SNESSetUp_VI - Sets up the internal data structures for the later use 1278 of the SNESVI nonlinear solver. 1279 1280 Input Parameter: 1281 . snes - the SNES context 1282 . x - the solution vector 1283 1284 Application Interface Routine: SNESSetUp() 1285 1286 Notes: 1287 For basic use of the SNES solvers, the user need not explicitly call 1288 SNESSetUp(), since these actions will automatically occur during 1289 the call to SNESSolve(). 1290 */ 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "SNESSetUp_VI" 1293 PetscErrorCode SNESSetUp_VI(SNES snes) 1294 { 1295 PetscErrorCode ierr; 1296 SNES_VI *vi = (SNES_VI*) snes->data; 1297 PetscInt i_start[3],i_end[3]; 1298 1299 PetscFunctionBegin; 1300 if (!snes->vec_sol_update) { 1301 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1302 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1303 } 1304 if (!snes->work) { 1305 snes->nwork = 3; 1306 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1307 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1308 } 1309 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1310 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1311 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1312 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1313 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1314 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1315 1316 /* If the lower and upper bound on variables are not set, set it to 1317 -Inf and Inf */ 1318 if (!vi->xl && !vi->xu) { 1319 vi->usersetxbounds = PETSC_FALSE; 1320 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1321 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1322 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1323 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1324 } else { 1325 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1326 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1327 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1328 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1329 if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2])) 1330 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Distribution of lower bound, upper bound and the solution vector should be identical across all the processors."); 1331 } 1332 1333 vi->computeuserfunction = snes->ops->computefunction; 1334 snes->ops->computefunction = SNESVIComputeFunction; 1335 1336 PetscFunctionReturn(0); 1337 } 1338 /* -------------------------------------------------------------------------- */ 1339 /* 1340 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1341 with SNESCreate_VI(). 1342 1343 Input Parameter: 1344 . snes - the SNES context 1345 1346 Application Interface Routine: SNESDestroy() 1347 */ 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "SNESDestroy_VI" 1350 PetscErrorCode SNESDestroy_VI(SNES snes) 1351 { 1352 SNES_VI *vi = (SNES_VI*) snes->data; 1353 PetscErrorCode ierr; 1354 1355 PetscFunctionBegin; 1356 if (snes->vec_sol_update) { 1357 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1358 snes->vec_sol_update = PETSC_NULL; 1359 } 1360 if (snes->nwork) { 1361 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1362 snes->nwork = 0; 1363 snes->work = PETSC_NULL; 1364 } 1365 1366 /* clear vectors */ 1367 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1368 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1369 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1370 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1371 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1372 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1373 if (!vi->usersetxbounds) { 1374 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1375 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1376 } 1377 if (vi->lsmonitor) { 1378 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1379 } 1380 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1381 1382 /* clear composed functions */ 1383 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1384 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1385 PetscFunctionReturn(0); 1386 } 1387 /* -------------------------------------------------------------------------- */ 1388 #undef __FUNCT__ 1389 #define __FUNCT__ "SNESLineSearchNo_VI" 1390 1391 /* 1392 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1393 1394 */ 1395 PetscErrorCode SNESLineSearchNo_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1396 { 1397 PetscErrorCode ierr; 1398 SNES_VI *vi = (SNES_VI*)snes->data; 1399 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1400 1401 PetscFunctionBegin; 1402 *flag = PETSC_TRUE; 1403 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1404 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1405 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1406 if (vi->postcheckstep) { 1407 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1408 } 1409 if (changed_y) { 1410 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1411 } 1412 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1413 if (!snes->domainerror) { 1414 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1415 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1416 } 1417 if (vi->lsmonitor) { 1418 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1419 } 1420 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1421 PetscFunctionReturn(0); 1422 } 1423 1424 /* -------------------------------------------------------------------------- */ 1425 #undef __FUNCT__ 1426 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1427 1428 /* 1429 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1430 */ 1431 PetscErrorCode SNESLineSearchNoNorms_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1432 { 1433 PetscErrorCode ierr; 1434 SNES_VI *vi = (SNES_VI*)snes->data; 1435 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1436 1437 PetscFunctionBegin; 1438 *flag = PETSC_TRUE; 1439 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1440 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1441 if (vi->postcheckstep) { 1442 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1443 } 1444 if (changed_y) { 1445 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1446 } 1447 1448 /* don't evaluate function the last time through */ 1449 if (snes->iter < snes->max_its-1) { 1450 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1451 } 1452 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1453 PetscFunctionReturn(0); 1454 } 1455 /* -------------------------------------------------------------------------- */ 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "SNESLineSearchCubic_VI" 1458 /* 1459 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1460 */ 1461 PetscErrorCode SNESLineSearchCubic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1462 { 1463 /* 1464 Note that for line search purposes we work with with the related 1465 minimization problem: 1466 min z(x): R^n -> R, 1467 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1468 This function z(x) is same as the merit function for the complementarity problem. 1469 */ 1470 1471 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1472 PetscReal minlambda,lambda,lambdatemp; 1473 #if defined(PETSC_USE_COMPLEX) 1474 PetscScalar cinitslope; 1475 #endif 1476 PetscErrorCode ierr; 1477 PetscInt count; 1478 SNES_VI *vi = (SNES_VI*)snes->data; 1479 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1480 MPI_Comm comm; 1481 1482 PetscFunctionBegin; 1483 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1484 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1485 *flag = PETSC_TRUE; 1486 1487 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1488 if (*ynorm == 0.0) { 1489 if (vi->lsmonitor) { 1490 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1491 } 1492 *gnorm = fnorm; 1493 ierr = VecCopy(x,w);CHKERRQ(ierr); 1494 ierr = VecCopy(f,g);CHKERRQ(ierr); 1495 *flag = PETSC_FALSE; 1496 goto theend1; 1497 } 1498 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1499 if (vi->lsmonitor) { 1500 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1501 } 1502 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1503 *ynorm = vi->maxstep; 1504 } 1505 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1506 minlambda = vi->minlambda/rellength; 1507 #if defined(PETSC_USE_COMPLEX) 1508 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1509 initslope = PetscRealPart(cinitslope); 1510 #else 1511 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1512 #endif 1513 if (initslope > 0.0) initslope = -initslope; 1514 if (initslope == 0.0) initslope = -1.0; 1515 1516 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1517 if (snes->nfuncs >= snes->max_funcs) { 1518 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1519 *flag = PETSC_FALSE; 1520 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1521 goto theend1; 1522 } 1523 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1524 if (snes->domainerror) { 1525 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1526 PetscFunctionReturn(0); 1527 } 1528 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1529 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1530 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1531 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1532 if (vi->lsmonitor) { 1533 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1534 } 1535 goto theend1; 1536 } 1537 1538 /* Fit points with quadratic */ 1539 lambda = 1.0; 1540 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1541 lambdaprev = lambda; 1542 gnormprev = *gnorm; 1543 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1544 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1545 else lambda = lambdatemp; 1546 1547 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1548 if (snes->nfuncs >= snes->max_funcs) { 1549 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1550 *flag = PETSC_FALSE; 1551 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1552 goto theend1; 1553 } 1554 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1555 if (snes->domainerror) { 1556 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1557 PetscFunctionReturn(0); 1558 } 1559 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1560 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1561 if (vi->lsmonitor) { 1562 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1563 } 1564 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1565 if (vi->lsmonitor) { 1566 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1567 } 1568 goto theend1; 1569 } 1570 1571 /* Fit points with cubic */ 1572 count = 1; 1573 while (PETSC_TRUE) { 1574 if (lambda <= minlambda) { 1575 if (vi->lsmonitor) { 1576 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1577 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,minlambda,lambda,initslope);CHKERRQ(ierr); 1578 } 1579 *flag = PETSC_FALSE; 1580 break; 1581 } 1582 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1583 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1584 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1585 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1586 d = b*b - 3*a*initslope; 1587 if (d < 0.0) d = 0.0; 1588 if (a == 0.0) { 1589 lambdatemp = -initslope/(2.0*b); 1590 } else { 1591 lambdatemp = (-b + sqrt(d))/(3.0*a); 1592 } 1593 lambdaprev = lambda; 1594 gnormprev = *gnorm; 1595 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1596 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1597 else lambda = lambdatemp; 1598 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1599 if (snes->nfuncs >= snes->max_funcs) { 1600 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1601 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1602 *flag = PETSC_FALSE; 1603 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1604 break; 1605 } 1606 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1607 if (snes->domainerror) { 1608 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1609 PetscFunctionReturn(0); 1610 } 1611 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1612 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1613 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1614 if (vi->lsmonitor) { 1615 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1616 } 1617 break; 1618 } else { 1619 if (vi->lsmonitor) { 1620 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1621 } 1622 } 1623 count++; 1624 } 1625 theend1: 1626 /* Optional user-defined check for line search step validity */ 1627 if (vi->postcheckstep && *flag) { 1628 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1629 if (changed_y) { 1630 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1631 } 1632 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1633 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1634 if (snes->domainerror) { 1635 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1636 PetscFunctionReturn(0); 1637 } 1638 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1639 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1640 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1641 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1642 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1643 } 1644 } 1645 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1646 PetscFunctionReturn(0); 1647 } 1648 /* -------------------------------------------------------------------------- */ 1649 #undef __FUNCT__ 1650 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1651 /* 1652 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1653 */ 1654 PetscErrorCode SNESLineSearchQuadratic_VI(SNES snes,void *lsctx,Vec x,Vec f,Vec g,Vec y,Vec w,PetscReal fnorm,PetscReal xnorm,PetscReal *ynorm,PetscReal *gnorm,PetscBool *flag) 1655 { 1656 /* 1657 Note that for line search purposes we work with with the related 1658 minimization problem: 1659 min z(x): R^n -> R, 1660 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1661 z(x) is the same as the merit function for the complementarity problem 1662 */ 1663 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1664 #if defined(PETSC_USE_COMPLEX) 1665 PetscScalar cinitslope; 1666 #endif 1667 PetscErrorCode ierr; 1668 PetscInt count; 1669 SNES_VI *vi = (SNES_VI*)snes->data; 1670 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1671 1672 PetscFunctionBegin; 1673 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1674 *flag = PETSC_TRUE; 1675 1676 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1677 if (*ynorm == 0.0) { 1678 if (vi->lsmonitor) { 1679 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1680 } 1681 *gnorm = fnorm; 1682 ierr = VecCopy(x,w);CHKERRQ(ierr); 1683 ierr = VecCopy(f,g);CHKERRQ(ierr); 1684 *flag = PETSC_FALSE; 1685 goto theend2; 1686 } 1687 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1688 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1689 *ynorm = vi->maxstep; 1690 } 1691 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1692 minlambda = vi->minlambda/rellength; 1693 #if defined(PETSC_USE_COMPLEX) 1694 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1695 initslope = PetscRealPart(cinitslope); 1696 #else 1697 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1698 #endif 1699 if (initslope > 0.0) initslope = -initslope; 1700 if (initslope == 0.0) initslope = -1.0; 1701 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1702 1703 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1704 if (snes->nfuncs >= snes->max_funcs) { 1705 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1706 *flag = PETSC_FALSE; 1707 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1708 goto theend2; 1709 } 1710 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1711 if (snes->domainerror) { 1712 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1713 PetscFunctionReturn(0); 1714 } 1715 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1716 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1717 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1718 if (vi->lsmonitor) { 1719 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1720 } 1721 goto theend2; 1722 } 1723 1724 /* Fit points with quadratic */ 1725 lambda = 1.0; 1726 count = 1; 1727 while (PETSC_TRUE) { 1728 if (lambda <= minlambda) { /* bad luck; use full step */ 1729 if (vi->lsmonitor) { 1730 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1731 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: fnorm=%G, gnorm=%G, ynorm=%G, lambda=%G, initial slope=%G\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1732 } 1733 ierr = VecCopy(x,w);CHKERRQ(ierr); 1734 *flag = PETSC_FALSE; 1735 break; 1736 } 1737 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1738 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1739 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1740 else lambda = lambdatemp; 1741 1742 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1743 if (snes->nfuncs >= snes->max_funcs) { 1744 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1745 ierr = PetscInfo5(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",fnorm,*gnorm,*ynorm,lambda,initslope);CHKERRQ(ierr); 1746 *flag = PETSC_FALSE; 1747 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1748 break; 1749 } 1750 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1751 if (snes->domainerror) { 1752 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1753 PetscFunctionReturn(0); 1754 } 1755 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1756 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1757 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1758 if (vi->lsmonitor) { 1759 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 1760 } 1761 break; 1762 } 1763 count++; 1764 } 1765 theend2: 1766 /* Optional user-defined check for line search step validity */ 1767 if (vi->postcheckstep) { 1768 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1769 if (changed_y) { 1770 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1771 } 1772 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1773 ierr = SNESComputeFunction(snes,w,g); 1774 if (snes->domainerror) { 1775 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1776 PetscFunctionReturn(0); 1777 } 1778 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1779 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1780 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1781 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1782 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1783 } 1784 } 1785 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1786 PetscFunctionReturn(0); 1787 } 1788 1789 typedef PetscErrorCode (*FCN2)(SNES,void*,Vec,Vec,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal*,PetscReal*,PetscBool*); /* force argument to next function to not be extern C*/ 1790 /* -------------------------------------------------------------------------- */ 1791 EXTERN_C_BEGIN 1792 #undef __FUNCT__ 1793 #define __FUNCT__ "SNESLineSearchSet_VI" 1794 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 1795 { 1796 PetscFunctionBegin; 1797 ((SNES_VI *)(snes->data))->LineSearch = func; 1798 ((SNES_VI *)(snes->data))->lsP = lsctx; 1799 PetscFunctionReturn(0); 1800 } 1801 EXTERN_C_END 1802 1803 /* -------------------------------------------------------------------------- */ 1804 EXTERN_C_BEGIN 1805 #undef __FUNCT__ 1806 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 1807 PetscErrorCode PETSCSNES_DLLEXPORT SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 1808 { 1809 SNES_VI *vi = (SNES_VI*)snes->data; 1810 PetscErrorCode ierr; 1811 1812 PetscFunctionBegin; 1813 if (flg && !vi->lsmonitor) { 1814 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 1815 } else if (!flg && vi->lsmonitor) { 1816 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1817 } 1818 PetscFunctionReturn(0); 1819 } 1820 EXTERN_C_END 1821 1822 /* 1823 SNESView_VI - Prints info from the SNESVI data structure. 1824 1825 Input Parameters: 1826 . SNES - the SNES context 1827 . viewer - visualization context 1828 1829 Application Interface Routine: SNESView() 1830 */ 1831 #undef __FUNCT__ 1832 #define __FUNCT__ "SNESView_VI" 1833 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 1834 { 1835 SNES_VI *vi = (SNES_VI *)snes->data; 1836 const char *cstr,*tstr; 1837 PetscErrorCode ierr; 1838 PetscBool iascii; 1839 1840 PetscFunctionBegin; 1841 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1842 if (iascii) { 1843 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 1844 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 1845 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 1846 else cstr = "unknown"; 1847 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 1848 else if (snes->ops->solve == SNESSolveVI_AS) tstr = "Active Set"; 1849 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 1850 else tstr = "unknown"; 1851 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 1852 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 1853 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 1854 } else { 1855 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 1856 } 1857 PetscFunctionReturn(0); 1858 } 1859 1860 /* 1861 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 1862 1863 Input Parameters: 1864 . snes - the SNES context. 1865 . xl - lower bound. 1866 . xu - upper bound. 1867 1868 Notes: 1869 If this routine is not called then the lower and upper bounds are set to 1870 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 1871 1872 */ 1873 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "SNESVISetVariableBounds" 1876 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 1877 { 1878 SNES_VI *vi = (SNES_VI*)snes->data; 1879 1880 PetscFunctionBegin; 1881 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1882 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1883 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 1884 1885 /* Check if SNESSetFunction is called */ 1886 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 1887 1888 /* Check if the vector sizes are compatible for lower and upper bounds */ 1889 if (xl->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths lower bound = %D solution vector = %D",xl->map->N,snes->vec_func->map->N); 1890 if (xu->map->N != snes->vec_func->map->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector lengths: upper bound = %D solution vector = %D",xu->map->N,snes->vec_func->map->N); 1891 vi->usersetxbounds = PETSC_TRUE; 1892 vi->xl = xl; 1893 vi->xu = xu; 1894 1895 PetscFunctionReturn(0); 1896 } 1897 /* -------------------------------------------------------------------------- */ 1898 /* 1899 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 1900 1901 Input Parameter: 1902 . snes - the SNES context 1903 1904 Application Interface Routine: SNESSetFromOptions() 1905 */ 1906 #undef __FUNCT__ 1907 #define __FUNCT__ "SNESSetFromOptions_VI" 1908 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 1909 { 1910 SNES_VI *vi = (SNES_VI *)snes->data; 1911 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 1912 const char *vies[] = {"ss","as","rs"}; 1913 PetscErrorCode ierr; 1914 PetscInt indx; 1915 PetscBool flg,set,flg2; 1916 1917 PetscFunctionBegin; 1918 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 1919 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 1920 if (flg) { 1921 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 1922 } 1923 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 1924 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 1925 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 1926 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 1927 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 1928 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 1929 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 1930 if (flg2) { 1931 switch (indx) { 1932 case 0: 1933 snes->ops->solve = SNESSolveVI_SS; 1934 break; 1935 case 1: 1936 snes->ops->solve = SNESSolveVI_AS; 1937 break; 1938 case 2: 1939 snes->ops->solve = SNESSolveVI_RS; 1940 break; 1941 } 1942 } 1943 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 1944 if (flg) { 1945 switch (indx) { 1946 case 0: 1947 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 1948 break; 1949 case 1: 1950 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 1951 break; 1952 case 2: 1953 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 1954 break; 1955 case 3: 1956 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 1957 break; 1958 } 1959 } 1960 ierr = PetscOptionsTail();CHKERRQ(ierr); 1961 PetscFunctionReturn(0); 1962 } 1963 /* -------------------------------------------------------------------------- */ 1964 /*MC 1965 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 1966 1967 Options Database: 1968 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 1969 . -snes_vi_alpha <alpha> - Sets alpha 1970 . -snes_vi_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y) 1971 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 1972 - -snes_vi_monitor - print information about progress of line searches 1973 1974 1975 Level: beginner 1976 1977 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 1978 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 1979 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 1980 1981 M*/ 1982 EXTERN_C_BEGIN 1983 #undef __FUNCT__ 1984 #define __FUNCT__ "SNESCreate_VI" 1985 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_VI(SNES snes) 1986 { 1987 PetscErrorCode ierr; 1988 SNES_VI *vi; 1989 1990 PetscFunctionBegin; 1991 snes->ops->setup = SNESSetUp_VI; 1992 snes->ops->solve = SNESSolveVI_SS; 1993 snes->ops->destroy = SNESDestroy_VI; 1994 snes->ops->setfromoptions = SNESSetFromOptions_VI; 1995 snes->ops->view = SNESView_VI; 1996 snes->ops->converged = SNESDefaultConverged_VI; 1997 1998 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 1999 snes->data = (void*)vi; 2000 vi->alpha = 1.e-4; 2001 vi->maxstep = 1.e8; 2002 vi->minlambda = 1.e-12; 2003 vi->LineSearch = SNESLineSearchCubic_VI; 2004 vi->lsP = PETSC_NULL; 2005 vi->postcheckstep = PETSC_NULL; 2006 vi->postcheck = PETSC_NULL; 2007 vi->precheckstep = PETSC_NULL; 2008 vi->precheck = PETSC_NULL; 2009 vi->const_tol = 2.2204460492503131e-16; 2010 vi->computessfunction = ComputeFischerFunction; 2011 2012 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2013 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2014 2015 PetscFunctionReturn(0); 2016 } 2017 EXTERN_C_END 2018