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 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 PetscInt i,nlocal,ilow,ihigh,nloc_isact=0,nloc_isinact=0; 742 PetscInt *idx_act,*idx_inact,i1=0,i2=0; 743 const PetscScalar *x,*xl,*xu,*f; 744 Vec F = snes->vec_func; 745 746 PetscFunctionBegin; 747 748 ierr = VecGetLocalSize(X,&nlocal);CHKERRQ(ierr); 749 ierr = VecGetOwnershipRange(X,&ilow,&ihigh);CHKERRQ(ierr); 750 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 751 ierr = VecGetArrayRead(Xl,&xl);CHKERRQ(ierr); 752 ierr = VecGetArrayRead(Xu,&xu);CHKERRQ(ierr); 753 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 754 /* Compute the sizes of the active and inactive sets */ 755 for (i=0; i < nlocal;i++) { 756 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) nloc_isinact++; 757 else nloc_isact++; 758 } 759 ierr = PetscMalloc(nloc_isact*sizeof(PetscInt),&idx_act);CHKERRQ(ierr); 760 ierr = PetscMalloc(nloc_isinact*sizeof(PetscInt),&idx_inact);CHKERRQ(ierr); 761 762 /* Creating the indexing arrays */ 763 for(i=0; i < nlocal; i++) { 764 if (((x[i] > xl[i] + 1.e-8 || (f[i] < 0.0)) && ((x[i] < xu[i] - 1.e-8) || f[i] > 0.0))) idx_inact[i2++] = ilow+i; 765 else idx_act[i1++] = ilow+i; 766 } 767 768 /* Create the index sets */ 769 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isact,idx_act,PETSC_COPY_VALUES,ISact);CHKERRQ(ierr); 770 ierr = ISCreateGeneral(((PetscObject)snes)->comm,nloc_isinact,idx_inact,PETSC_COPY_VALUES,ISinact);CHKERRQ(ierr); 771 772 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 773 ierr = VecRestoreArrayRead(Xl,&xl);CHKERRQ(ierr); 774 ierr = VecRestoreArrayRead(Xu,&xu);CHKERRQ(ierr); 775 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 776 ierr = PetscFree(idx_act);CHKERRQ(ierr); 777 ierr = PetscFree(idx_inact);CHKERRQ(ierr); 778 PetscFunctionReturn(0); 779 } 780 781 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 782 #undef __FUNCT__ 783 #define __FUNCT__ "SNESVICreateSubVectors" 784 PetscErrorCode SNESVICreateSubVectors(SNES snes,PetscInt n,Vec* newv) 785 { 786 PetscErrorCode ierr; 787 Vec v; 788 789 PetscFunctionBegin; 790 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 791 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 792 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 793 *newv = v; 794 795 PetscFunctionReturn(0); 796 } 797 798 /* Resets the snes PC and KSP when the active set sizes change */ 799 #undef __FUNCT__ 800 #define __FUNCT__ "SNESVIResetPCandKSP" 801 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 802 { 803 PetscErrorCode ierr; 804 KSP kspnew,snesksp; 805 PC pcnew; 806 const MatSolverPackage stype; 807 808 PetscFunctionBegin; 809 /* The active and inactive set sizes have changed so need to create a new snes->ksp object */ 810 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 811 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 812 /* Copy over snes->ksp info */ 813 kspnew->pc_side = snesksp->pc_side; 814 kspnew->rtol = snesksp->rtol; 815 kspnew->abstol = snesksp->abstol; 816 kspnew->max_it = snesksp->max_it; 817 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 818 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 819 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 820 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 821 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 822 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 823 ierr = KSPDestroy(snesksp);CHKERRQ(ierr); 824 snes->ksp = kspnew; 825 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 826 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr); 827 PetscFunctionReturn(0); 828 } 829 /* Variational Inequality solver using active set method */ 830 #undef __FUNCT__ 831 #define __FUNCT__ "SNESSolveVI_AS" 832 PetscErrorCode SNESSolveVI_AS(SNES snes) 833 { 834 SNES_VI *vi = (SNES_VI*)snes->data; 835 PetscErrorCode ierr; 836 PetscInt maxits,i,lits,Nis_act=0; 837 PetscBool lssucceed; 838 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 839 PetscReal gnorm,xnorm=0,ynorm; 840 Vec Y,X,F,G,W; 841 KSPConvergedReason kspreason; 842 843 PetscFunctionBegin; 844 snes->numFailures = 0; 845 snes->numLinearSolveFailures = 0; 846 snes->reason = SNES_CONVERGED_ITERATING; 847 848 maxits = snes->max_its; /* maximum number of iterations */ 849 X = snes->vec_sol; /* solution vector */ 850 F = snes->vec_func; /* residual vector */ 851 Y = snes->work[0]; /* work vectors */ 852 G = snes->work[1]; 853 W = snes->work[2]; 854 855 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 856 snes->iter = 0; 857 snes->norm = 0.0; 858 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 859 860 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 861 ierr = SNESComputeFunction(snes,X,vi->phi);CHKERRQ(ierr); 862 if (snes->domainerror) { 863 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 864 PetscFunctionReturn(0); 865 } 866 /* Compute Merit function */ 867 ierr = SNESVIComputeMeritFunction(vi->phi,&vi->merit,&vi->phinorm);CHKERRQ(ierr); 868 869 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 870 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 871 if PetscIsInfOrNanReal(vi->merit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 872 873 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 874 snes->norm = vi->phinorm; 875 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 876 SNESLogConvHistory(snes,vi->phinorm,0); 877 SNESMonitor(snes,0,vi->phinorm); 878 879 /* set parameter for default relative tolerance convergence test */ 880 snes->ttol = vi->phinorm*snes->rtol; 881 /* test convergence */ 882 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 883 if (snes->reason) PetscFunctionReturn(0); 884 885 for (i=0; i<maxits; i++) { 886 887 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 888 PetscReal thresh,J_norm1; 889 VecScatter scat_act,scat_inact; 890 PetscInt nis_act,nis_inact,Nis_act_prev; 891 Vec Da_act,Da_inact,Db_inact; 892 Vec Y_act,Y_inact,phi_act,phi_inact; 893 Mat jac_inact_inact,jac_inact_act,prejac_inact_inact; 894 895 /* Call general purpose update function */ 896 if (snes->ops->update) { 897 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 898 } 899 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 900 /* Compute the threshold value for creating active and inactive sets */ 901 ierr = MatNorm(snes->jacobian,NORM_1,&J_norm1);CHKERRQ(ierr); 902 thresh = PetscMin(vi->merit,1e-2)/(1+J_norm1); 903 904 /* Compute B-subdifferential vectors Da and Db */ 905 ierr = SNESVIComputeBsubdifferentialVectors(snes,X,F,snes->jacobian,vi->Da,vi->Db);CHKERRQ(ierr); 906 907 /* Create active and inactive index sets */ 908 ierr = SNESVICreateIndexSets_AS(snes,vi->Db,thresh,&IS_act,&IS_inact);CHKERRQ(ierr); 909 910 Nis_act_prev = Nis_act; 911 /* Get sizes of active and inactive sets */ 912 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 913 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 914 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 915 916 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 917 918 /* Create active and inactive set vectors */ 919 ierr = SNESVICreateSubVectors(snes,nis_act,&Da_act);CHKERRQ(ierr); 920 ierr = SNESVICreateSubVectors(snes,nis_inact,&Da_inact);CHKERRQ(ierr); 921 ierr = SNESVICreateSubVectors(snes,nis_inact,&Db_inact);CHKERRQ(ierr); 922 ierr = SNESVICreateSubVectors(snes,nis_act,&phi_act);CHKERRQ(ierr); 923 ierr = SNESVICreateSubVectors(snes,nis_inact,&phi_inact);CHKERRQ(ierr); 924 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 925 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 926 927 /* Create inactive set submatrices */ 928 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_act,MAT_INITIAL_MATRIX,&jac_inact_act);CHKERRQ(ierr); 929 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 930 931 /* Create scatter contexts */ 932 ierr = VecScatterCreate(vi->Da,IS_act,Da_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 933 ierr = VecScatterCreate(vi->Da,IS_inact,Da_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 934 935 /* Do a vec scatter to active and inactive set vectors */ 936 ierr = VecScatterBegin(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 937 ierr = VecScatterEnd(scat_act,vi->Da,Da_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 938 939 ierr = VecScatterBegin(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 940 ierr = VecScatterEnd(scat_inact,vi->Da,Da_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 941 942 ierr = VecScatterBegin(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 943 ierr = VecScatterEnd(scat_inact,vi->Db,Db_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 944 945 ierr = VecScatterBegin(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 946 ierr = VecScatterEnd(scat_act,vi->phi,phi_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 947 948 ierr = VecScatterBegin(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 949 ierr = VecScatterEnd(scat_inact,vi->phi,phi_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 950 951 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 952 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 953 954 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 955 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 956 957 /* Active set direction */ 958 ierr = VecPointwiseDivide(Y_act,phi_act,Da_act);CHKERRQ(ierr); 959 /* inactive set jacobian and preconditioner */ 960 ierr = VecPointwiseDivide(Da_inact,Da_inact,Db_inact);CHKERRQ(ierr); 961 ierr = MatDiagonalSet(jac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 962 if (snes->jacobian != snes->jacobian_pre) { 963 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 964 ierr = MatDiagonalSet(prejac_inact_inact,Da_inact,ADD_VALUES);CHKERRQ(ierr); 965 } else prejac_inact_inact = jac_inact_inact; 966 967 /* right hand side */ 968 ierr = VecPointwiseDivide(phi_inact,phi_inact,Db_inact);CHKERRQ(ierr); 969 ierr = MatMult(jac_inact_act,Y_act,Db_inact);CHKERRQ(ierr); 970 ierr = VecAXPY(phi_inact,-1.0,Db_inact);CHKERRQ(ierr); 971 972 if ((i != 0) && (Nis_act != Nis_act_prev)) { 973 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 974 } 975 976 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 977 ierr = SNES_KSPSolve(snes,snes->ksp,phi_inact,Y_inact);CHKERRQ(ierr); 978 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 979 if (kspreason < 0) { 980 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 981 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 982 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 983 break; 984 } 985 } 986 987 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 988 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 989 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 990 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 991 992 ierr = VecDestroy(Da_act);CHKERRQ(ierr); 993 ierr = VecDestroy(Da_inact);CHKERRQ(ierr); 994 ierr = VecDestroy(Db_inact);CHKERRQ(ierr); 995 ierr = VecDestroy(phi_act);CHKERRQ(ierr); 996 ierr = VecDestroy(phi_inact);CHKERRQ(ierr); 997 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 998 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 999 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1000 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1001 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1002 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1003 ierr = MatDestroy(jac_inact_act);CHKERRQ(ierr); 1004 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1005 if (snes->jacobian != snes->jacobian_pre) { 1006 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1007 } 1008 1009 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1010 snes->linear_its += lits; 1011 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1012 /* 1013 if (vi->precheckstep) { 1014 PetscBool changed_y = PETSC_FALSE; 1015 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1016 } 1017 1018 if (PetscLogPrintInfo){ 1019 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1020 } 1021 */ 1022 /* Compute a (scaled) negative update in the line search routine: 1023 Y <- X - lambda*Y 1024 and evaluate G = function(Y) (depends on the line search). 1025 */ 1026 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1027 ynorm = 1; gnorm = vi->phinorm; 1028 ierr = SNESVIComputeJacobian(snes->jacobian,snes->jacobian_pre,vi->Da,vi->Db);CHKERRQ(ierr); 1029 ierr = SNESVIComputeMeritFunctionGradient(snes->jacobian,vi->phi,vi->dpsi);CHKERRQ(ierr); 1030 ierr = (*vi->LineSearch)(snes,vi->lsP,X,vi->phi,G,Y,W,vi->phinorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1031 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",vi->phinorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1032 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1033 if (snes->domainerror) { 1034 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1035 PetscFunctionReturn(0); 1036 } 1037 if (!lssucceed) { 1038 if (++snes->numFailures >= snes->maxFailures) { 1039 PetscBool ismin; 1040 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1041 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1042 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1043 break; 1044 } 1045 } 1046 /* Update function and solution vectors */ 1047 vi->phinorm = gnorm; 1048 vi->merit = 0.5*vi->phinorm*vi->phinorm; 1049 ierr = VecCopy(G,vi->phi);CHKERRQ(ierr); 1050 ierr = VecCopy(W,X);CHKERRQ(ierr); 1051 /* Monitor convergence */ 1052 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1053 snes->iter = i+1; 1054 snes->norm = vi->phinorm; 1055 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1056 SNESLogConvHistory(snes,snes->norm,lits); 1057 SNESMonitor(snes,snes->iter,snes->norm); 1058 /* Test for convergence, xnorm = || X || */ 1059 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1060 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,vi->phinorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1061 if (snes->reason) break; 1062 } 1063 if (i == maxits) { 1064 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1065 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1066 } 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "SNESVIComputeInactiveSetFnorm" 1072 PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes,Vec F,Vec X,PetscScalar *fnorm) 1073 { 1074 PetscErrorCode ierr; 1075 SNES_VI *vi = (SNES_VI*)snes->data; 1076 const PetscScalar *x,*xl,*xu,*f; 1077 PetscInt i,n; 1078 PetscReal rnorm; 1079 1080 PetscFunctionBegin; 1081 1082 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 1083 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1084 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1085 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 1086 ierr = VecGetArrayRead(F,&f);CHKERRQ(ierr); 1087 1088 rnorm = 0.0; 1089 for (i=0; i<n; i++) { 1090 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]; 1091 } 1092 ierr = VecRestoreArrayRead(F,&f);CHKERRQ(ierr); 1093 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1094 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1095 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 1096 ierr = MPI_Allreduce(&rnorm,fnorm,1,MPIU_REAL,MPI_SUM,((PetscObject)snes)->comm);CHKERRQ(ierr); 1097 *fnorm = sqrt(*fnorm); 1098 PetscFunctionReturn(0); 1099 } 1100 1101 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 1102 implemented in this algorithm. It basically identifies the active variables and does 1103 a linear solve on the inactive variables. */ 1104 #undef __FUNCT__ 1105 #define __FUNCT__ "SNESSolveVI_RS" 1106 PetscErrorCode SNESSolveVI_RS(SNES snes) 1107 { 1108 SNES_VI *vi = (SNES_VI*)snes->data; 1109 PetscErrorCode ierr; 1110 PetscInt maxits,i,lits,Nis_act=0; 1111 PetscBool lssucceed; 1112 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1113 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1114 Vec Y,X,F,G,W; 1115 KSPConvergedReason kspreason; 1116 1117 PetscFunctionBegin; 1118 snes->numFailures = 0; 1119 snes->numLinearSolveFailures = 0; 1120 snes->reason = SNES_CONVERGED_ITERATING; 1121 1122 maxits = snes->max_its; /* maximum number of iterations */ 1123 X = snes->vec_sol; /* solution vector */ 1124 F = snes->vec_func; /* residual vector */ 1125 Y = snes->work[0]; /* work vectors */ 1126 G = snes->work[1]; 1127 W = snes->work[2]; 1128 1129 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1130 snes->iter = 0; 1131 snes->norm = 0.0; 1132 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1133 1134 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 1135 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1136 if (snes->domainerror) { 1137 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1138 PetscFunctionReturn(0); 1139 } 1140 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1141 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1142 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1143 if PetscIsInfOrNanReal(fnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1144 1145 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1146 snes->norm = fnorm; 1147 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1148 SNESLogConvHistory(snes,fnorm,0); 1149 SNESMonitor(snes,0,fnorm); 1150 1151 /* set parameter for default relative tolerance convergence test */ 1152 snes->ttol = fnorm*snes->rtol; 1153 /* test convergence */ 1154 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1155 if (snes->reason) PetscFunctionReturn(0); 1156 1157 for (i=0; i<maxits; i++) { 1158 1159 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1160 VecScatter scat_act,scat_inact; 1161 PetscInt nis_act,nis_inact,Nis_act_prev; 1162 Vec Y_act,Y_inact,F_inact; 1163 Mat jac_inact_inact,prejac_inact_inact; 1164 1165 /* Call general purpose update function */ 1166 if (snes->ops->update) { 1167 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1168 } 1169 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1170 /* Create active and inactive index sets */ 1171 ierr = SNESVICreateIndexSets_RS(snes,X,vi->xl,vi->xu,&IS_act,&IS_inact);CHKERRQ(ierr); 1172 1173 Nis_act_prev = Nis_act; 1174 /* Get sizes of active and inactive sets */ 1175 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1176 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1177 ierr = ISGetSize(IS_act,&Nis_act);CHKERRQ(ierr); 1178 1179 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 1180 1181 /* Create active and inactive set vectors */ 1182 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 1183 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 1184 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1185 1186 /* Create inactive set submatrices */ 1187 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1188 1189 /* Create scatter contexts */ 1190 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1191 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1192 1193 /* Do a vec scatter to active and inactive set vectors */ 1194 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1195 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1196 1197 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1198 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1199 1200 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1201 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1202 1203 /* Active set direction = 0 */ 1204 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1205 if (snes->jacobian != snes->jacobian_pre) { 1206 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1207 } else prejac_inact_inact = jac_inact_inact; 1208 1209 if ((i != 0) && (Nis_act != Nis_act_prev)) { 1210 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1211 } 1212 1213 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1214 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 1215 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1216 if (kspreason < 0) { 1217 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1218 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1219 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1220 break; 1221 } 1222 } 1223 1224 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1225 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1226 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1227 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1228 1229 ierr = VecDestroy(F_inact);CHKERRQ(ierr); 1230 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 1231 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 1232 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1233 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1234 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1235 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1236 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1237 if (snes->jacobian != snes->jacobian_pre) { 1238 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1239 } 1240 1241 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1242 snes->linear_its += lits; 1243 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1244 /* 1245 if (vi->precheckstep) { 1246 PetscBool changed_y = PETSC_FALSE; 1247 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1248 } 1249 1250 if (PetscLogPrintInfo){ 1251 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1252 } 1253 */ 1254 /* Compute a (scaled) negative update in the line search routine: 1255 Y <- X - lambda*Y 1256 and evaluate G = function(Y) (depends on the line search). 1257 */ 1258 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1259 ynorm = 1; gnorm = fnorm; 1260 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1261 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1262 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1263 if (snes->domainerror) { 1264 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1265 PetscFunctionReturn(0); 1266 } 1267 if (!lssucceed) { 1268 if (++snes->numFailures >= snes->maxFailures) { 1269 PetscBool ismin; 1270 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1271 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1272 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1273 break; 1274 } 1275 } 1276 /* Update function and solution vectors */ 1277 fnorm = gnorm; 1278 ierr = VecCopy(G,F);CHKERRQ(ierr); 1279 ierr = VecCopy(W,X);CHKERRQ(ierr); 1280 /* Monitor convergence */ 1281 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1282 snes->iter = i+1; 1283 snes->norm = fnorm; 1284 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1285 SNESLogConvHistory(snes,snes->norm,lits); 1286 SNESMonitor(snes,snes->iter,snes->norm); 1287 /* Test for convergence, xnorm = || X || */ 1288 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1289 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1290 if (snes->reason) break; 1291 } 1292 if (i == maxits) { 1293 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1294 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1295 } 1296 PetscFunctionReturn(0); 1297 } 1298 1299 /* -------------------------------------------------------------------------- */ 1300 /* 1301 SNESSetUp_VI - Sets up the internal data structures for the later use 1302 of the SNESVI nonlinear solver. 1303 1304 Input Parameter: 1305 . snes - the SNES context 1306 . x - the solution vector 1307 1308 Application Interface Routine: SNESSetUp() 1309 1310 Notes: 1311 For basic use of the SNES solvers, the user need not explicitly call 1312 SNESSetUp(), since these actions will automatically occur during 1313 the call to SNESSolve(). 1314 */ 1315 #undef __FUNCT__ 1316 #define __FUNCT__ "SNESSetUp_VI" 1317 PetscErrorCode SNESSetUp_VI(SNES snes) 1318 { 1319 PetscErrorCode ierr; 1320 SNES_VI *vi = (SNES_VI*) snes->data; 1321 PetscInt i_start[3],i_end[3]; 1322 1323 PetscFunctionBegin; 1324 if (!snes->vec_sol_update) { 1325 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1326 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1327 } 1328 if (!snes->work) { 1329 snes->nwork = 3; 1330 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1331 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1332 } 1333 1334 /* If the lower and upper bound on variables are not set, set it to 1335 -Inf and Inf */ 1336 if (!vi->xl && !vi->xu) { 1337 vi->usersetxbounds = PETSC_FALSE; 1338 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1339 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1340 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1341 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1342 } else { 1343 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1344 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1345 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1346 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1347 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])) 1348 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."); 1349 } 1350 1351 if (snes->ops->solve != SNESSolveVI_RS) { 1352 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1353 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1354 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1355 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1356 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1357 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1358 1359 1360 vi->computeuserfunction = snes->ops->computefunction; 1361 snes->ops->computefunction = SNESVIComputeFunction; 1362 } 1363 PetscFunctionReturn(0); 1364 } 1365 /* -------------------------------------------------------------------------- */ 1366 /* 1367 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1368 with SNESCreate_VI(). 1369 1370 Input Parameter: 1371 . snes - the SNES context 1372 1373 Application Interface Routine: SNESDestroy() 1374 */ 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "SNESDestroy_VI" 1377 PetscErrorCode SNESDestroy_VI(SNES snes) 1378 { 1379 SNES_VI *vi = (SNES_VI*) snes->data; 1380 PetscErrorCode ierr; 1381 1382 PetscFunctionBegin; 1383 if (snes->vec_sol_update) { 1384 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1385 snes->vec_sol_update = PETSC_NULL; 1386 } 1387 if (snes->nwork) { 1388 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1389 snes->nwork = 0; 1390 snes->work = PETSC_NULL; 1391 } 1392 1393 if (snes->ops->solve != SNESSolveVI_RS) { 1394 /* clear vectors */ 1395 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1396 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1397 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1398 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1399 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1400 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1401 if (!vi->usersetxbounds) { 1402 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1403 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1404 } 1405 } 1406 if (vi->lsmonitor) { 1407 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1408 } 1409 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1410 1411 /* clear composed functions */ 1412 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1413 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1414 PetscFunctionReturn(0); 1415 } 1416 /* -------------------------------------------------------------------------- */ 1417 #undef __FUNCT__ 1418 #define __FUNCT__ "SNESLineSearchNo_VI" 1419 1420 /* 1421 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1422 1423 */ 1424 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) 1425 { 1426 PetscErrorCode ierr; 1427 SNES_VI *vi = (SNES_VI*)snes->data; 1428 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1429 1430 PetscFunctionBegin; 1431 *flag = PETSC_TRUE; 1432 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1433 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1434 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1435 if (vi->postcheckstep) { 1436 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1437 } 1438 if (changed_y) { 1439 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1440 } 1441 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1442 if (!snes->domainerror) { 1443 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1444 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1445 } 1446 if (vi->lsmonitor) { 1447 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1448 } 1449 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1450 PetscFunctionReturn(0); 1451 } 1452 1453 /* -------------------------------------------------------------------------- */ 1454 #undef __FUNCT__ 1455 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1456 1457 /* 1458 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1459 */ 1460 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) 1461 { 1462 PetscErrorCode ierr; 1463 SNES_VI *vi = (SNES_VI*)snes->data; 1464 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1465 1466 PetscFunctionBegin; 1467 *flag = PETSC_TRUE; 1468 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1469 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1470 if (vi->postcheckstep) { 1471 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1472 } 1473 if (changed_y) { 1474 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1475 } 1476 1477 /* don't evaluate function the last time through */ 1478 if (snes->iter < snes->max_its-1) { 1479 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1480 } 1481 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1482 PetscFunctionReturn(0); 1483 } 1484 /* -------------------------------------------------------------------------- */ 1485 #undef __FUNCT__ 1486 #define __FUNCT__ "SNESLineSearchCubic_VI" 1487 /* 1488 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1489 */ 1490 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) 1491 { 1492 /* 1493 Note that for line search purposes we work with with the related 1494 minimization problem: 1495 min z(x): R^n -> R, 1496 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1497 This function z(x) is same as the merit function for the complementarity problem. 1498 */ 1499 1500 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1501 PetscReal minlambda,lambda,lambdatemp; 1502 #if defined(PETSC_USE_COMPLEX) 1503 PetscScalar cinitslope; 1504 #endif 1505 PetscErrorCode ierr; 1506 PetscInt count; 1507 SNES_VI *vi = (SNES_VI*)snes->data; 1508 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1509 MPI_Comm comm; 1510 1511 PetscFunctionBegin; 1512 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1513 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1514 *flag = PETSC_TRUE; 1515 1516 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1517 if (*ynorm == 0.0) { 1518 if (vi->lsmonitor) { 1519 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1520 } 1521 *gnorm = fnorm; 1522 ierr = VecCopy(x,w);CHKERRQ(ierr); 1523 ierr = VecCopy(f,g);CHKERRQ(ierr); 1524 *flag = PETSC_FALSE; 1525 goto theend1; 1526 } 1527 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1528 if (vi->lsmonitor) { 1529 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1530 } 1531 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1532 *ynorm = vi->maxstep; 1533 } 1534 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1535 minlambda = vi->minlambda/rellength; 1536 #if defined(PETSC_USE_COMPLEX) 1537 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1538 initslope = PetscRealPart(cinitslope); 1539 #else 1540 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1541 #endif 1542 if (initslope > 0.0) initslope = -initslope; 1543 if (initslope == 0.0) initslope = -1.0; 1544 1545 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1546 if (snes->nfuncs >= snes->max_funcs) { 1547 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1548 *flag = PETSC_FALSE; 1549 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1550 goto theend1; 1551 } 1552 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1553 if (snes->domainerror) { 1554 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1555 PetscFunctionReturn(0); 1556 } 1557 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1558 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1559 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1560 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1561 if (vi->lsmonitor) { 1562 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1563 } 1564 goto theend1; 1565 } 1566 1567 /* Fit points with quadratic */ 1568 lambda = 1.0; 1569 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1570 lambdaprev = lambda; 1571 gnormprev = *gnorm; 1572 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1573 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1574 else lambda = lambdatemp; 1575 1576 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1577 if (snes->nfuncs >= snes->max_funcs) { 1578 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1579 *flag = PETSC_FALSE; 1580 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1581 goto theend1; 1582 } 1583 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1584 if (snes->domainerror) { 1585 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1586 PetscFunctionReturn(0); 1587 } 1588 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1589 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1590 if (vi->lsmonitor) { 1591 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1592 } 1593 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1594 if (vi->lsmonitor) { 1595 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1596 } 1597 goto theend1; 1598 } 1599 1600 /* Fit points with cubic */ 1601 count = 1; 1602 while (PETSC_TRUE) { 1603 if (lambda <= minlambda) { 1604 if (vi->lsmonitor) { 1605 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1606 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); 1607 } 1608 *flag = PETSC_FALSE; 1609 break; 1610 } 1611 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1612 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1613 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1614 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1615 d = b*b - 3*a*initslope; 1616 if (d < 0.0) d = 0.0; 1617 if (a == 0.0) { 1618 lambdatemp = -initslope/(2.0*b); 1619 } else { 1620 lambdatemp = (-b + sqrt(d))/(3.0*a); 1621 } 1622 lambdaprev = lambda; 1623 gnormprev = *gnorm; 1624 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1625 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1626 else lambda = lambdatemp; 1627 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1628 if (snes->nfuncs >= snes->max_funcs) { 1629 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1630 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); 1631 *flag = PETSC_FALSE; 1632 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1633 break; 1634 } 1635 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1636 if (snes->domainerror) { 1637 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1638 PetscFunctionReturn(0); 1639 } 1640 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1641 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1642 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1643 if (vi->lsmonitor) { 1644 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1645 } 1646 break; 1647 } else { 1648 if (vi->lsmonitor) { 1649 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1650 } 1651 } 1652 count++; 1653 } 1654 theend1: 1655 /* Optional user-defined check for line search step validity */ 1656 if (vi->postcheckstep && *flag) { 1657 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1658 if (changed_y) { 1659 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1660 } 1661 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1662 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1663 if (snes->domainerror) { 1664 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1668 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1669 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1670 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1671 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1672 } 1673 } 1674 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1675 PetscFunctionReturn(0); 1676 } 1677 /* -------------------------------------------------------------------------- */ 1678 /* 1679 Projects vector X onto the variable bounds such that 1680 Xl[i] <= X[i] <= Xu[i] for i = 1..n 1681 */ 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "SNESVIProjectOntoBounds" 1685 PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X) 1686 { 1687 PetscErrorCode ierr; 1688 SNES_VI *vi = (SNES_VI*)snes->data; 1689 const PetscScalar *xl,*xu; 1690 PetscScalar *x; 1691 PetscInt i,n; 1692 1693 PetscFunctionBegin; 1694 1695 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 1696 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 1697 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1698 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1699 1700 for(i = 0;i<n;i++) { 1701 if (x[i] < xl[i]) x[i] = xl[i]; 1702 else if (x[i] > xu[i]) x[i] = xu[i]; 1703 } 1704 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 1705 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1706 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1707 1708 PetscFunctionReturn(0); 1709 } 1710 1711 /* -------------------------------------------------------------------------- */ 1712 #undef __FUNCT__ 1713 #define __FUNCT__ "SNESLineSearchBoundProjectedCubic_VI" 1714 /* 1715 This routine implements a bound projected cubic line search 1716 */ 1717 PetscErrorCode SNESLineSearchBoundProjectedCubic_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) 1718 { 1719 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1720 PetscReal minlambda,lambda,lambdatemp; 1721 #if defined(PETSC_USE_COMPLEX) 1722 PetscScalar cinitslope; 1723 #endif 1724 PetscErrorCode ierr; 1725 PetscInt count; 1726 SNES_VI *vi = (SNES_VI*)snes->data; 1727 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1728 MPI_Comm comm; 1729 1730 PetscFunctionBegin; 1731 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1732 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1733 *flag = PETSC_TRUE; 1734 1735 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1736 if (*ynorm == 0.0) { 1737 if (vi->lsmonitor) { 1738 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1739 } 1740 *gnorm = fnorm; 1741 ierr = VecCopy(x,w);CHKERRQ(ierr); 1742 ierr = VecCopy(f,g);CHKERRQ(ierr); 1743 *flag = PETSC_FALSE; 1744 goto theend1; 1745 } 1746 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1747 if (vi->lsmonitor) { 1748 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1749 } 1750 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1751 *ynorm = vi->maxstep; 1752 } 1753 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1754 minlambda = vi->minlambda/rellength; 1755 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1756 #if defined(PETSC_USE_COMPLEX) 1757 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1758 initslope = PetscRealPart(cinitslope); 1759 #else 1760 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1761 #endif 1762 if (initslope > 0.0) initslope = -initslope; 1763 if (initslope == 0.0) initslope = -1.0; 1764 1765 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1766 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1767 if (snes->nfuncs >= snes->max_funcs) { 1768 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1769 *flag = PETSC_FALSE; 1770 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1771 goto theend1; 1772 } 1773 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1774 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1775 if (snes->domainerror) { 1776 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1777 PetscFunctionReturn(0); 1778 } 1779 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1780 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1781 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1782 if (vi->lsmonitor) { 1783 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1784 } 1785 goto theend1; 1786 } 1787 1788 /* Fit points with quadratic */ 1789 lambda = 1.0; 1790 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1791 lambdaprev = lambda; 1792 gnormprev = *gnorm; 1793 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1794 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1795 else lambda = lambdatemp; 1796 1797 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1798 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1799 if (snes->nfuncs >= snes->max_funcs) { 1800 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1801 *flag = PETSC_FALSE; 1802 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1803 goto theend1; 1804 } 1805 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1806 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1807 if (snes->domainerror) { 1808 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1809 PetscFunctionReturn(0); 1810 } 1811 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1812 if (vi->lsmonitor) { 1813 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1814 } 1815 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1816 if (vi->lsmonitor) { 1817 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1818 } 1819 goto theend1; 1820 } 1821 1822 /* Fit points with cubic */ 1823 count = 1; 1824 while (PETSC_TRUE) { 1825 if (lambda <= minlambda) { 1826 if (vi->lsmonitor) { 1827 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1828 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); 1829 } 1830 *flag = PETSC_FALSE; 1831 break; 1832 } 1833 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1834 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1835 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1836 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1837 d = b*b - 3*a*initslope; 1838 if (d < 0.0) d = 0.0; 1839 if (a == 0.0) { 1840 lambdatemp = -initslope/(2.0*b); 1841 } else { 1842 lambdatemp = (-b + sqrt(d))/(3.0*a); 1843 } 1844 lambdaprev = lambda; 1845 gnormprev = *gnorm; 1846 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1847 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1848 else lambda = lambdatemp; 1849 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1850 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1851 if (snes->nfuncs >= snes->max_funcs) { 1852 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1853 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); 1854 *flag = PETSC_FALSE; 1855 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1856 break; 1857 } 1858 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1859 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1860 if (snes->domainerror) { 1861 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1865 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1866 if (vi->lsmonitor) { 1867 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1868 } 1869 break; 1870 } else { 1871 if (vi->lsmonitor) { 1872 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1873 } 1874 } 1875 count++; 1876 } 1877 theend1: 1878 /* Optional user-defined check for line search step validity */ 1879 if (vi->postcheckstep && *flag) { 1880 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1881 if (changed_y) { 1882 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1883 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1884 } 1885 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1886 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1887 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1888 if (snes->domainerror) { 1889 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1893 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1894 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1895 } 1896 } 1897 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1898 PetscFunctionReturn(0); 1899 } 1900 1901 #undef __FUNCT__ 1902 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1903 /* 1904 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1905 */ 1906 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) 1907 { 1908 /* 1909 Note that for line search purposes we work with with the related 1910 minimization problem: 1911 min z(x): R^n -> R, 1912 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1913 z(x) is the same as the merit function for the complementarity problem 1914 */ 1915 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1916 #if defined(PETSC_USE_COMPLEX) 1917 PetscScalar cinitslope; 1918 #endif 1919 PetscErrorCode ierr; 1920 PetscInt count; 1921 SNES_VI *vi = (SNES_VI*)snes->data; 1922 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1923 1924 PetscFunctionBegin; 1925 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1926 *flag = PETSC_TRUE; 1927 1928 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1929 if (*ynorm == 0.0) { 1930 if (vi->lsmonitor) { 1931 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1932 } 1933 *gnorm = fnorm; 1934 ierr = VecCopy(x,w);CHKERRQ(ierr); 1935 ierr = VecCopy(f,g);CHKERRQ(ierr); 1936 *flag = PETSC_FALSE; 1937 goto theend2; 1938 } 1939 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1940 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1941 *ynorm = vi->maxstep; 1942 } 1943 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1944 minlambda = vi->minlambda/rellength; 1945 #if defined(PETSC_USE_COMPLEX) 1946 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1947 initslope = PetscRealPart(cinitslope); 1948 #else 1949 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1950 #endif 1951 if (initslope > 0.0) initslope = -initslope; 1952 if (initslope == 0.0) initslope = -1.0; 1953 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1954 1955 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1956 if (snes->nfuncs >= snes->max_funcs) { 1957 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1958 *flag = PETSC_FALSE; 1959 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1960 goto theend2; 1961 } 1962 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1963 if (snes->domainerror) { 1964 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1965 PetscFunctionReturn(0); 1966 } 1967 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1968 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1969 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1970 if (vi->lsmonitor) { 1971 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1972 } 1973 goto theend2; 1974 } 1975 1976 /* Fit points with quadratic */ 1977 lambda = 1.0; 1978 count = 1; 1979 while (PETSC_TRUE) { 1980 if (lambda <= minlambda) { /* bad luck; use full step */ 1981 if (vi->lsmonitor) { 1982 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1983 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); 1984 } 1985 ierr = VecCopy(x,w);CHKERRQ(ierr); 1986 *flag = PETSC_FALSE; 1987 break; 1988 } 1989 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1990 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1991 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1992 else lambda = lambdatemp; 1993 1994 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1995 if (snes->nfuncs >= snes->max_funcs) { 1996 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1997 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); 1998 *flag = PETSC_FALSE; 1999 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2000 break; 2001 } 2002 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2003 if (snes->domainerror) { 2004 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2008 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2009 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 2010 if (vi->lsmonitor) { 2011 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 2012 } 2013 break; 2014 } 2015 count++; 2016 } 2017 theend2: 2018 /* Optional user-defined check for line search step validity */ 2019 if (vi->postcheckstep) { 2020 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2021 if (changed_y) { 2022 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2023 } 2024 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2025 ierr = SNESComputeFunction(snes,w,g); 2026 if (snes->domainerror) { 2027 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2028 PetscFunctionReturn(0); 2029 } 2030 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 2031 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2032 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 2033 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2034 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2035 } 2036 } 2037 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2038 PetscFunctionReturn(0); 2039 } 2040 2041 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*/ 2042 /* -------------------------------------------------------------------------- */ 2043 EXTERN_C_BEGIN 2044 #undef __FUNCT__ 2045 #define __FUNCT__ "SNESLineSearchSet_VI" 2046 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 2047 { 2048 PetscFunctionBegin; 2049 ((SNES_VI *)(snes->data))->LineSearch = func; 2050 ((SNES_VI *)(snes->data))->lsP = lsctx; 2051 PetscFunctionReturn(0); 2052 } 2053 EXTERN_C_END 2054 2055 /* -------------------------------------------------------------------------- */ 2056 EXTERN_C_BEGIN 2057 #undef __FUNCT__ 2058 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 2059 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 2060 { 2061 SNES_VI *vi = (SNES_VI*)snes->data; 2062 PetscErrorCode ierr; 2063 2064 PetscFunctionBegin; 2065 if (flg && !vi->lsmonitor) { 2066 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 2067 } else if (!flg && vi->lsmonitor) { 2068 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 2069 } 2070 PetscFunctionReturn(0); 2071 } 2072 EXTERN_C_END 2073 2074 /* 2075 SNESView_VI - Prints info from the SNESVI data structure. 2076 2077 Input Parameters: 2078 . SNES - the SNES context 2079 . viewer - visualization context 2080 2081 Application Interface Routine: SNESView() 2082 */ 2083 #undef __FUNCT__ 2084 #define __FUNCT__ "SNESView_VI" 2085 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 2086 { 2087 SNES_VI *vi = (SNES_VI *)snes->data; 2088 const char *cstr,*tstr; 2089 PetscErrorCode ierr; 2090 PetscBool iascii; 2091 2092 PetscFunctionBegin; 2093 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2094 if (iascii) { 2095 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2096 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2097 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2098 else if (vi->LineSearch == SNESLineSearchBoundProjectedCubic_VI) cstr = "SNESLineSearchBoundProjectedCubic"; 2099 else cstr = "unknown"; 2100 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2101 else if (snes->ops->solve == SNESSolveVI_AS) tstr = "Active Set"; 2102 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2103 else tstr = "unknown"; 2104 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2105 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2106 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2107 } else { 2108 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2109 } 2110 PetscFunctionReturn(0); 2111 } 2112 2113 /* 2114 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2115 2116 Input Parameters: 2117 . snes - the SNES context. 2118 . xl - lower bound. 2119 . xu - upper bound. 2120 2121 Notes: 2122 If this routine is not called then the lower and upper bounds are set to 2123 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 2124 2125 */ 2126 2127 #undef __FUNCT__ 2128 #define __FUNCT__ "SNESVISetVariableBounds" 2129 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2130 { 2131 SNES_VI *vi = (SNES_VI*)snes->data; 2132 2133 PetscFunctionBegin; 2134 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2135 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2136 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2137 2138 /* Check if SNESSetFunction is called */ 2139 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2140 2141 /* Check if the vector sizes are compatible for lower and upper bounds */ 2142 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); 2143 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); 2144 vi->usersetxbounds = PETSC_TRUE; 2145 vi->xl = xl; 2146 vi->xu = xu; 2147 2148 PetscFunctionReturn(0); 2149 } 2150 /* -------------------------------------------------------------------------- */ 2151 /* 2152 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2153 2154 Input Parameter: 2155 . snes - the SNES context 2156 2157 Application Interface Routine: SNESSetFromOptions() 2158 */ 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "SNESSetFromOptions_VI" 2161 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2162 { 2163 SNES_VI *vi = (SNES_VI *)snes->data; 2164 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2165 const char *vies[] = {"ss","as","rs"}; 2166 PetscErrorCode ierr; 2167 PetscInt indx; 2168 PetscBool flg,set,flg2; 2169 2170 PetscFunctionBegin; 2171 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2172 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2173 if (flg) { 2174 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2175 } 2176 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2177 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2178 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2179 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2180 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2181 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2182 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2183 if (flg2) { 2184 switch (indx) { 2185 case 0: 2186 snes->ops->solve = SNESSolveVI_SS; 2187 break; 2188 case 1: 2189 snes->ops->solve = SNESSolveVI_AS; 2190 break; 2191 case 2: 2192 snes->ops->solve = SNESSolveVI_RS; 2193 break; 2194 } 2195 } 2196 if (snes->ops->solve == SNESSolveVI_RS) { 2197 ierr = SNESLineSearchSet(snes,SNESLineSearchBoundProjectedCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2198 } 2199 else { 2200 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 2201 if (flg) { 2202 switch (indx) { 2203 case 0: 2204 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2205 break; 2206 case 1: 2207 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2208 break; 2209 case 2: 2210 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2211 break; 2212 case 3: 2213 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2214 break; 2215 } 2216 } 2217 } 2218 ierr = PetscOptionsTail();CHKERRQ(ierr); 2219 PetscFunctionReturn(0); 2220 } 2221 /* -------------------------------------------------------------------------- */ 2222 /*MC 2223 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 2224 2225 Options Database: 2226 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 2227 . -snes_vi_alpha <alpha> - Sets alpha 2228 . -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) 2229 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 2230 - -snes_vi_monitor - print information about progress of line searches 2231 2232 2233 Level: beginner 2234 2235 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2236 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2237 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2238 2239 M*/ 2240 EXTERN_C_BEGIN 2241 #undef __FUNCT__ 2242 #define __FUNCT__ "SNESCreate_VI" 2243 PetscErrorCode SNESCreate_VI(SNES snes) 2244 { 2245 PetscErrorCode ierr; 2246 SNES_VI *vi; 2247 2248 PetscFunctionBegin; 2249 snes->ops->setup = SNESSetUp_VI; 2250 snes->ops->solve = SNESSolveVI_SS; 2251 snes->ops->destroy = SNESDestroy_VI; 2252 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2253 snes->ops->view = SNESView_VI; 2254 snes->ops->converged = SNESDefaultConverged_VI; 2255 2256 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2257 snes->data = (void*)vi; 2258 vi->alpha = 1.e-4; 2259 vi->maxstep = 1.e8; 2260 vi->minlambda = 1.e-12; 2261 vi->LineSearch = SNESLineSearchCubic_VI; 2262 vi->lsP = PETSC_NULL; 2263 vi->postcheckstep = PETSC_NULL; 2264 vi->postcheck = PETSC_NULL; 2265 vi->precheckstep = PETSC_NULL; 2266 vi->precheck = PETSC_NULL; 2267 vi->const_tol = 2.2204460492503131e-16; 2268 vi->computessfunction = ComputeFischerFunction; 2269 2270 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2271 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2272 2273 PetscFunctionReturn(0); 2274 } 2275 EXTERN_C_END 2276