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_inact; 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 Nis_inact = F->map->N; 1130 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1131 snes->iter = 0; 1132 snes->norm = 0.0; 1133 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1134 1135 ierr = SNESVIAdjustInitialGuess(X,vi->xl,vi->xu);CHKERRQ(ierr); 1136 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1137 if (snes->domainerror) { 1138 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1139 PetscFunctionReturn(0); 1140 } 1141 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1142 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1143 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1144 if PetscIsInfOrNanReal(fnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1145 1146 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1147 snes->norm = fnorm; 1148 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1149 SNESLogConvHistory(snes,fnorm,0); 1150 SNESMonitor(snes,0,fnorm); 1151 1152 /* set parameter for default relative tolerance convergence test */ 1153 snes->ttol = fnorm*snes->rtol; 1154 /* test convergence */ 1155 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1156 if (snes->reason) PetscFunctionReturn(0); 1157 1158 for (i=0; i<maxits; i++) { 1159 1160 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1161 VecScatter scat_act,scat_inact; 1162 PetscInt nis_act,nis_inact,Nis_inact_prev; 1163 Vec Y_act,Y_inact,F_inact; 1164 Mat jac_inact_inact,prejac_inact_inact; 1165 1166 /* Call general purpose update function */ 1167 if (snes->ops->update) { 1168 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1169 } 1170 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1171 /* Create active and inactive index sets */ 1172 ierr = SNESVICreateIndexSets_RS(snes,X,vi->xl,vi->xu,&IS_act,&IS_inact);CHKERRQ(ierr); 1173 1174 Nis_inact_prev = Nis_inact; 1175 /* Get sizes of active and inactive sets */ 1176 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1177 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1178 ierr = ISGetSize(IS_inact,&Nis_inact);CHKERRQ(ierr); 1179 1180 ierr = PetscPrintf(PETSC_COMM_WORLD,"Size of active set = %d, size of inactive set = %d\n",nis_act,nis_inact);CHKERRQ(ierr); 1181 1182 /* Create active and inactive set vectors */ 1183 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 1184 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 1185 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1186 1187 /* Create inactive set submatrices */ 1188 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1189 1190 /* Create scatter contexts */ 1191 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1192 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1193 1194 /* Do a vec scatter to active and inactive set vectors */ 1195 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1196 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1197 1198 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1199 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1200 1201 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1202 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1203 1204 /* Active set direction = 0 */ 1205 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1206 if (snes->jacobian != snes->jacobian_pre) { 1207 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1208 } else prejac_inact_inact = jac_inact_inact; 1209 1210 if (Nis_inact != Nis_inact_prev) { 1211 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1212 } 1213 1214 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1215 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 1216 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1217 if (kspreason < 0) { 1218 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1219 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1220 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1221 break; 1222 } 1223 } 1224 1225 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1226 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1227 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1228 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1229 1230 ierr = VecDestroy(F_inact);CHKERRQ(ierr); 1231 ierr = VecDestroy(Y_act);CHKERRQ(ierr); 1232 ierr = VecDestroy(Y_inact);CHKERRQ(ierr); 1233 ierr = VecScatterDestroy(scat_act);CHKERRQ(ierr); 1234 ierr = VecScatterDestroy(scat_inact);CHKERRQ(ierr); 1235 ierr = ISDestroy(IS_act);CHKERRQ(ierr); 1236 ierr = ISDestroy(IS_inact);CHKERRQ(ierr); 1237 ierr = MatDestroy(jac_inact_inact);CHKERRQ(ierr); 1238 if (snes->jacobian != snes->jacobian_pre) { 1239 ierr = MatDestroy(prejac_inact_inact);CHKERRQ(ierr); 1240 } 1241 1242 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1243 snes->linear_its += lits; 1244 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1245 /* 1246 if (vi->precheckstep) { 1247 PetscBool changed_y = PETSC_FALSE; 1248 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1249 } 1250 1251 if (PetscLogPrintInfo){ 1252 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1253 } 1254 */ 1255 /* Compute a (scaled) negative update in the line search routine: 1256 Y <- X - lambda*Y 1257 and evaluate G = function(Y) (depends on the line search). 1258 */ 1259 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1260 ynorm = 1; gnorm = fnorm; 1261 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1262 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1263 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1264 if (snes->domainerror) { 1265 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1266 PetscFunctionReturn(0); 1267 } 1268 if (!lssucceed) { 1269 if (++snes->numFailures >= snes->maxFailures) { 1270 PetscBool ismin; 1271 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1272 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1273 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1274 break; 1275 } 1276 } 1277 /* Update function and solution vectors */ 1278 fnorm = gnorm; 1279 ierr = VecCopy(G,F);CHKERRQ(ierr); 1280 ierr = VecCopy(W,X);CHKERRQ(ierr); 1281 /* Monitor convergence */ 1282 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1283 snes->iter = i+1; 1284 snes->norm = fnorm; 1285 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1286 SNESLogConvHistory(snes,snes->norm,lits); 1287 SNESMonitor(snes,snes->iter,snes->norm); 1288 /* Test for convergence, xnorm = || X || */ 1289 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1290 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1291 if (snes->reason) break; 1292 } 1293 if (i == maxits) { 1294 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1295 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1296 } 1297 PetscFunctionReturn(0); 1298 } 1299 1300 /* -------------------------------------------------------------------------- */ 1301 /* 1302 SNESSetUp_VI - Sets up the internal data structures for the later use 1303 of the SNESVI nonlinear solver. 1304 1305 Input Parameter: 1306 . snes - the SNES context 1307 . x - the solution vector 1308 1309 Application Interface Routine: SNESSetUp() 1310 1311 Notes: 1312 For basic use of the SNES solvers, the user need not explicitly call 1313 SNESSetUp(), since these actions will automatically occur during 1314 the call to SNESSolve(). 1315 */ 1316 #undef __FUNCT__ 1317 #define __FUNCT__ "SNESSetUp_VI" 1318 PetscErrorCode SNESSetUp_VI(SNES snes) 1319 { 1320 PetscErrorCode ierr; 1321 SNES_VI *vi = (SNES_VI*) snes->data; 1322 PetscInt i_start[3],i_end[3]; 1323 1324 PetscFunctionBegin; 1325 if (!snes->vec_sol_update) { 1326 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 1327 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 1328 } 1329 if (!snes->work) { 1330 snes->nwork = 3; 1331 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 1332 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 1333 } 1334 1335 /* If the lower and upper bound on variables are not set, set it to 1336 -Inf and Inf */ 1337 if (!vi->xl && !vi->xu) { 1338 vi->usersetxbounds = PETSC_FALSE; 1339 ierr = VecDuplicate(snes->vec_sol, &vi->xl); CHKERRQ(ierr); 1340 ierr = VecSet(vi->xl,PETSC_VI_NINF);CHKERRQ(ierr); 1341 ierr = VecDuplicate(snes->vec_sol, &vi->xu); CHKERRQ(ierr); 1342 ierr = VecSet(vi->xu,PETSC_VI_INF);CHKERRQ(ierr); 1343 } else { 1344 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1345 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1346 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1347 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1348 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])) 1349 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."); 1350 } 1351 1352 if (snes->ops->solve != SNESSolveVI_RS) { 1353 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1354 ierr = VecDuplicate(snes->vec_sol, &vi->phi); CHKERRQ(ierr); 1355 ierr = VecDuplicate(snes->vec_sol, &vi->Da); CHKERRQ(ierr); 1356 ierr = VecDuplicate(snes->vec_sol, &vi->Db); CHKERRQ(ierr); 1357 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1358 ierr = VecDuplicate(snes->vec_sol, &vi->t); CHKERRQ(ierr); 1359 1360 1361 vi->computeuserfunction = snes->ops->computefunction; 1362 snes->ops->computefunction = SNESVIComputeFunction; 1363 } 1364 PetscFunctionReturn(0); 1365 } 1366 /* -------------------------------------------------------------------------- */ 1367 /* 1368 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1369 with SNESCreate_VI(). 1370 1371 Input Parameter: 1372 . snes - the SNES context 1373 1374 Application Interface Routine: SNESDestroy() 1375 */ 1376 #undef __FUNCT__ 1377 #define __FUNCT__ "SNESDestroy_VI" 1378 PetscErrorCode SNESDestroy_VI(SNES snes) 1379 { 1380 SNES_VI *vi = (SNES_VI*) snes->data; 1381 PetscErrorCode ierr; 1382 1383 PetscFunctionBegin; 1384 if (snes->vec_sol_update) { 1385 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 1386 snes->vec_sol_update = PETSC_NULL; 1387 } 1388 if (snes->nwork) { 1389 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 1390 snes->nwork = 0; 1391 snes->work = PETSC_NULL; 1392 } 1393 1394 if (snes->ops->solve != SNESSolveVI_RS) { 1395 /* clear vectors */ 1396 ierr = VecDestroy(vi->dpsi);CHKERRQ(ierr); 1397 ierr = VecDestroy(vi->phi); CHKERRQ(ierr); 1398 ierr = VecDestroy(vi->Da); CHKERRQ(ierr); 1399 ierr = VecDestroy(vi->Db); CHKERRQ(ierr); 1400 ierr = VecDestroy(vi->z); CHKERRQ(ierr); 1401 ierr = VecDestroy(vi->t); CHKERRQ(ierr); 1402 if (!vi->usersetxbounds) { 1403 ierr = VecDestroy(vi->xl); CHKERRQ(ierr); 1404 ierr = VecDestroy(vi->xu); CHKERRQ(ierr); 1405 } 1406 } 1407 if (vi->lsmonitor) { 1408 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 1409 } 1410 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1411 1412 /* clear composed functions */ 1413 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1414 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1415 PetscFunctionReturn(0); 1416 } 1417 /* -------------------------------------------------------------------------- */ 1418 #undef __FUNCT__ 1419 #define __FUNCT__ "SNESLineSearchNo_VI" 1420 1421 /* 1422 This routine is a copy of SNESLineSearchNo routine in snes/impls/ls/ls.c 1423 1424 */ 1425 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) 1426 { 1427 PetscErrorCode ierr; 1428 SNES_VI *vi = (SNES_VI*)snes->data; 1429 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1430 1431 PetscFunctionBegin; 1432 *flag = PETSC_TRUE; 1433 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1434 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1435 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1436 if (vi->postcheckstep) { 1437 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1438 } 1439 if (changed_y) { 1440 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1441 } 1442 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1443 if (!snes->domainerror) { 1444 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1445 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1446 } 1447 if (vi->lsmonitor) { 1448 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1449 } 1450 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 /* -------------------------------------------------------------------------- */ 1455 #undef __FUNCT__ 1456 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1457 1458 /* 1459 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1460 */ 1461 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) 1462 { 1463 PetscErrorCode ierr; 1464 SNES_VI *vi = (SNES_VI*)snes->data; 1465 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1466 1467 PetscFunctionBegin; 1468 *flag = PETSC_TRUE; 1469 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1470 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1471 if (vi->postcheckstep) { 1472 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1473 } 1474 if (changed_y) { 1475 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1476 } 1477 1478 /* don't evaluate function the last time through */ 1479 if (snes->iter < snes->max_its-1) { 1480 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1481 } 1482 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1483 PetscFunctionReturn(0); 1484 } 1485 /* -------------------------------------------------------------------------- */ 1486 #undef __FUNCT__ 1487 #define __FUNCT__ "SNESLineSearchCubic_VI" 1488 /* 1489 This routine is a copy of SNESLineSearchCubic in snes/impls/ls/ls.c 1490 */ 1491 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) 1492 { 1493 /* 1494 Note that for line search purposes we work with with the related 1495 minimization problem: 1496 min z(x): R^n -> R, 1497 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 1498 This function z(x) is same as the merit function for the complementarity problem. 1499 */ 1500 1501 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1502 PetscReal minlambda,lambda,lambdatemp; 1503 #if defined(PETSC_USE_COMPLEX) 1504 PetscScalar cinitslope; 1505 #endif 1506 PetscErrorCode ierr; 1507 PetscInt count; 1508 SNES_VI *vi = (SNES_VI*)snes->data; 1509 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1510 MPI_Comm comm; 1511 1512 PetscFunctionBegin; 1513 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1514 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1515 *flag = PETSC_TRUE; 1516 1517 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1518 if (*ynorm == 0.0) { 1519 if (vi->lsmonitor) { 1520 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1521 } 1522 *gnorm = fnorm; 1523 ierr = VecCopy(x,w);CHKERRQ(ierr); 1524 ierr = VecCopy(f,g);CHKERRQ(ierr); 1525 *flag = PETSC_FALSE; 1526 goto theend1; 1527 } 1528 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1529 if (vi->lsmonitor) { 1530 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1531 } 1532 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1533 *ynorm = vi->maxstep; 1534 } 1535 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1536 minlambda = vi->minlambda/rellength; 1537 #if defined(PETSC_USE_COMPLEX) 1538 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1539 initslope = PetscRealPart(cinitslope); 1540 #else 1541 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1542 #endif 1543 if (initslope > 0.0) initslope = -initslope; 1544 if (initslope == 0.0) initslope = -1.0; 1545 1546 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1547 if (snes->nfuncs >= snes->max_funcs) { 1548 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1549 *flag = PETSC_FALSE; 1550 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1551 goto theend1; 1552 } 1553 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1554 if (snes->domainerror) { 1555 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1556 PetscFunctionReturn(0); 1557 } 1558 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1559 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1560 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1561 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1562 if (vi->lsmonitor) { 1563 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1564 } 1565 goto theend1; 1566 } 1567 1568 /* Fit points with quadratic */ 1569 lambda = 1.0; 1570 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1571 lambdaprev = lambda; 1572 gnormprev = *gnorm; 1573 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1574 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1575 else lambda = lambdatemp; 1576 1577 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1578 if (snes->nfuncs >= snes->max_funcs) { 1579 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1580 *flag = PETSC_FALSE; 1581 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1582 goto theend1; 1583 } 1584 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1585 if (snes->domainerror) { 1586 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1587 PetscFunctionReturn(0); 1588 } 1589 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1590 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1591 if (vi->lsmonitor) { 1592 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1593 } 1594 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1595 if (vi->lsmonitor) { 1596 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1597 } 1598 goto theend1; 1599 } 1600 1601 /* Fit points with cubic */ 1602 count = 1; 1603 while (PETSC_TRUE) { 1604 if (lambda <= minlambda) { 1605 if (vi->lsmonitor) { 1606 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1607 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); 1608 } 1609 *flag = PETSC_FALSE; 1610 break; 1611 } 1612 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1613 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1614 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1615 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1616 d = b*b - 3*a*initslope; 1617 if (d < 0.0) d = 0.0; 1618 if (a == 0.0) { 1619 lambdatemp = -initslope/(2.0*b); 1620 } else { 1621 lambdatemp = (-b + sqrt(d))/(3.0*a); 1622 } 1623 lambdaprev = lambda; 1624 gnormprev = *gnorm; 1625 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1626 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1627 else lambda = lambdatemp; 1628 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1629 if (snes->nfuncs >= snes->max_funcs) { 1630 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1631 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); 1632 *flag = PETSC_FALSE; 1633 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1634 break; 1635 } 1636 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1637 if (snes->domainerror) { 1638 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1639 PetscFunctionReturn(0); 1640 } 1641 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1642 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1643 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1644 if (vi->lsmonitor) { 1645 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1646 } 1647 break; 1648 } else { 1649 if (vi->lsmonitor) { 1650 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1651 } 1652 } 1653 count++; 1654 } 1655 theend1: 1656 /* Optional user-defined check for line search step validity */ 1657 if (vi->postcheckstep && *flag) { 1658 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1659 if (changed_y) { 1660 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1661 } 1662 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1663 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1664 if (snes->domainerror) { 1665 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 1669 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1670 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1671 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 1672 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1673 } 1674 } 1675 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 /* -------------------------------------------------------------------------- */ 1679 /* 1680 Projects vector X onto the variable bounds such that 1681 Xl[i] <= X[i] <= Xu[i] for i = 1..n 1682 */ 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "SNESVIProjectOntoBounds" 1686 PetscErrorCode SNESVIProjectOntoBounds(SNES snes,Vec X) 1687 { 1688 PetscErrorCode ierr; 1689 SNES_VI *vi = (SNES_VI*)snes->data; 1690 const PetscScalar *xl,*xu; 1691 PetscScalar *x; 1692 PetscInt i,n; 1693 1694 PetscFunctionBegin; 1695 1696 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr); 1697 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 1698 ierr = VecGetArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1699 ierr = VecGetArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1700 1701 for(i = 0;i<n;i++) { 1702 if (x[i] < xl[i]) x[i] = xl[i]; 1703 else if (x[i] > xu[i]) x[i] = xu[i]; 1704 } 1705 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 1706 ierr = VecRestoreArrayRead(vi->xl,&xl);CHKERRQ(ierr); 1707 ierr = VecRestoreArrayRead(vi->xu,&xu);CHKERRQ(ierr); 1708 1709 PetscFunctionReturn(0); 1710 } 1711 1712 /* -------------------------------------------------------------------------- */ 1713 #undef __FUNCT__ 1714 #define __FUNCT__ "SNESLineSearchBoundProjectedCubic_VI" 1715 /* 1716 This routine implements a bound projected cubic line search 1717 */ 1718 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) 1719 { 1720 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1721 PetscReal minlambda,lambda,lambdatemp; 1722 #if defined(PETSC_USE_COMPLEX) 1723 PetscScalar cinitslope; 1724 #endif 1725 PetscErrorCode ierr; 1726 PetscInt count; 1727 SNES_VI *vi = (SNES_VI*)snes->data; 1728 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1729 MPI_Comm comm; 1730 1731 PetscFunctionBegin; 1732 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1733 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1734 *flag = PETSC_TRUE; 1735 1736 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1737 if (*ynorm == 0.0) { 1738 if (vi->lsmonitor) { 1739 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1740 } 1741 *gnorm = fnorm; 1742 ierr = VecCopy(x,w);CHKERRQ(ierr); 1743 ierr = VecCopy(f,g);CHKERRQ(ierr); 1744 *flag = PETSC_FALSE; 1745 goto theend1; 1746 } 1747 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1748 if (vi->lsmonitor) { 1749 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1750 } 1751 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1752 *ynorm = vi->maxstep; 1753 } 1754 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1755 minlambda = vi->minlambda/rellength; 1756 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1757 #if defined(PETSC_USE_COMPLEX) 1758 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1759 initslope = PetscRealPart(cinitslope); 1760 #else 1761 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1762 #endif 1763 if (initslope > 0.0) initslope = -initslope; 1764 if (initslope == 0.0) initslope = -1.0; 1765 1766 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1767 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1768 if (snes->nfuncs >= snes->max_funcs) { 1769 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1770 *flag = PETSC_FALSE; 1771 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1772 goto theend1; 1773 } 1774 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1775 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1776 if (snes->domainerror) { 1777 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1778 PetscFunctionReturn(0); 1779 } 1780 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1781 ierr = PetscInfo2(snes,"Initial fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1782 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1783 if (vi->lsmonitor) { 1784 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1785 } 1786 goto theend1; 1787 } 1788 1789 /* Fit points with quadratic */ 1790 lambda = 1.0; 1791 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1792 lambdaprev = lambda; 1793 gnormprev = *gnorm; 1794 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1795 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1796 else lambda = lambdatemp; 1797 1798 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1799 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1800 if (snes->nfuncs >= snes->max_funcs) { 1801 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 1802 *flag = PETSC_FALSE; 1803 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1804 goto theend1; 1805 } 1806 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1807 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1808 if (snes->domainerror) { 1809 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1810 PetscFunctionReturn(0); 1811 } 1812 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1813 if (vi->lsmonitor) { 1814 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 1815 } 1816 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 1817 if (vi->lsmonitor) { 1818 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 1819 } 1820 goto theend1; 1821 } 1822 1823 /* Fit points with cubic */ 1824 count = 1; 1825 while (PETSC_TRUE) { 1826 if (lambda <= minlambda) { 1827 if (vi->lsmonitor) { 1828 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 1829 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); 1830 } 1831 *flag = PETSC_FALSE; 1832 break; 1833 } 1834 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 1835 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 1836 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1837 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 1838 d = b*b - 3*a*initslope; 1839 if (d < 0.0) d = 0.0; 1840 if (a == 0.0) { 1841 lambdatemp = -initslope/(2.0*b); 1842 } else { 1843 lambdatemp = (-b + sqrt(d))/(3.0*a); 1844 } 1845 lambdaprev = lambda; 1846 gnormprev = *gnorm; 1847 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1848 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1849 else lambda = lambdatemp; 1850 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1851 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1852 if (snes->nfuncs >= snes->max_funcs) { 1853 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1854 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); 1855 *flag = PETSC_FALSE; 1856 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1857 break; 1858 } 1859 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1860 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1861 if (snes->domainerror) { 1862 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1863 PetscFunctionReturn(0); 1864 } 1865 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1866 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* is reduction enough? */ 1867 if (vi->lsmonitor) { 1868 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1869 } 1870 break; 1871 } else { 1872 if (vi->lsmonitor) { 1873 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 1874 } 1875 } 1876 count++; 1877 } 1878 theend1: 1879 /* Optional user-defined check for line search step validity */ 1880 if (vi->postcheckstep && *flag) { 1881 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1882 if (changed_y) { 1883 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1884 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1885 } 1886 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 1887 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1888 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1889 if (snes->domainerror) { 1890 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1891 PetscFunctionReturn(0); 1892 } 1893 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1894 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 1895 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 1896 } 1897 } 1898 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1899 PetscFunctionReturn(0); 1900 } 1901 1902 #undef __FUNCT__ 1903 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 1904 /* 1905 This routine is a copy of SNESLineSearchQuadratic in snes/impls/ls/ls.c 1906 */ 1907 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) 1908 { 1909 /* 1910 Note that for line search purposes we work with with the related 1911 minimization problem: 1912 min z(x): R^n -> R, 1913 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 1914 z(x) is the same as the merit function for the complementarity problem 1915 */ 1916 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 1917 #if defined(PETSC_USE_COMPLEX) 1918 PetscScalar cinitslope; 1919 #endif 1920 PetscErrorCode ierr; 1921 PetscInt count; 1922 SNES_VI *vi = (SNES_VI*)snes->data; 1923 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1924 1925 PetscFunctionBegin; 1926 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1927 *flag = PETSC_TRUE; 1928 1929 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1930 if (*ynorm == 0.0) { 1931 if (vi->lsmonitor) { 1932 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 1933 } 1934 *gnorm = fnorm; 1935 ierr = VecCopy(x,w);CHKERRQ(ierr); 1936 ierr = VecCopy(f,g);CHKERRQ(ierr); 1937 *flag = PETSC_FALSE; 1938 goto theend2; 1939 } 1940 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1941 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1942 *ynorm = vi->maxstep; 1943 } 1944 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1945 minlambda = vi->minlambda/rellength; 1946 #if defined(PETSC_USE_COMPLEX) 1947 ierr = VecDot(vi->dpsi,y,&cinitslope);CHKERRQ(ierr); 1948 initslope = PetscRealPart(cinitslope); 1949 #else 1950 ierr = VecDot(vi->dpsi,y,&initslope);CHKERRQ(ierr); 1951 #endif 1952 if (initslope > 0.0) initslope = -initslope; 1953 if (initslope == 0.0) initslope = -1.0; 1954 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 1955 1956 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 1957 if (snes->nfuncs >= snes->max_funcs) { 1958 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 1959 *flag = PETSC_FALSE; 1960 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 1961 goto theend2; 1962 } 1963 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1964 if (snes->domainerror) { 1965 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1966 PetscFunctionReturn(0); 1967 } 1968 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 1969 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1970 if (.5*(*gnorm)*(*gnorm) <= .5*fnorm*fnorm + vi->alpha*initslope) { /* Sufficient reduction */ 1971 if (vi->lsmonitor) { 1972 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1973 } 1974 goto theend2; 1975 } 1976 1977 /* Fit points with quadratic */ 1978 lambda = 1.0; 1979 count = 1; 1980 while (PETSC_TRUE) { 1981 if (lambda <= minlambda) { /* bad luck; use full step */ 1982 if (vi->lsmonitor) { 1983 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 1984 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); 1985 } 1986 ierr = VecCopy(x,w);CHKERRQ(ierr); 1987 *flag = PETSC_FALSE; 1988 break; 1989 } 1990 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 1991 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 1992 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 1993 else lambda = lambdatemp; 1994 1995 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 1996 if (snes->nfuncs >= snes->max_funcs) { 1997 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 1998 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); 1999 *flag = PETSC_FALSE; 2000 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2001 break; 2002 } 2003 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2004 if (snes->domainerror) { 2005 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2006 PetscFunctionReturn(0); 2007 } 2008 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2009 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2010 if (.5*(*gnorm)*(*gnorm) < .5*fnorm*fnorm + lambda*vi->alpha*initslope) { /* sufficient reduction */ 2011 if (vi->lsmonitor) { 2012 ierr = PetscViewerASCIIMonitorPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 2013 } 2014 break; 2015 } 2016 count++; 2017 } 2018 theend2: 2019 /* Optional user-defined check for line search step validity */ 2020 if (vi->postcheckstep) { 2021 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2022 if (changed_y) { 2023 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2024 } 2025 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2026 ierr = SNESComputeFunction(snes,w,g); 2027 if (snes->domainerror) { 2028 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 ierr = VecNormBegin(g,NORM_2,gnorm);CHKERRQ(ierr); 2032 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2033 ierr = VecNormEnd(g,NORM_2,gnorm);CHKERRQ(ierr); 2034 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2035 if PetscIsInfOrNanReal(*gnorm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2036 } 2037 } 2038 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2039 PetscFunctionReturn(0); 2040 } 2041 2042 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*/ 2043 /* -------------------------------------------------------------------------- */ 2044 EXTERN_C_BEGIN 2045 #undef __FUNCT__ 2046 #define __FUNCT__ "SNESLineSearchSet_VI" 2047 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 2048 { 2049 PetscFunctionBegin; 2050 ((SNES_VI *)(snes->data))->LineSearch = func; 2051 ((SNES_VI *)(snes->data))->lsP = lsctx; 2052 PetscFunctionReturn(0); 2053 } 2054 EXTERN_C_END 2055 2056 /* -------------------------------------------------------------------------- */ 2057 EXTERN_C_BEGIN 2058 #undef __FUNCT__ 2059 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 2060 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 2061 { 2062 SNES_VI *vi = (SNES_VI*)snes->data; 2063 PetscErrorCode ierr; 2064 2065 PetscFunctionBegin; 2066 if (flg && !vi->lsmonitor) { 2067 ierr = PetscViewerASCIIMonitorCreate(((PetscObject)snes)->comm,"stdout",((PetscObject)snes)->tablevel,&vi->lsmonitor);CHKERRQ(ierr); 2068 } else if (!flg && vi->lsmonitor) { 2069 ierr = PetscViewerASCIIMonitorDestroy(vi->lsmonitor);CHKERRQ(ierr); 2070 } 2071 PetscFunctionReturn(0); 2072 } 2073 EXTERN_C_END 2074 2075 /* 2076 SNESView_VI - Prints info from the SNESVI data structure. 2077 2078 Input Parameters: 2079 . SNES - the SNES context 2080 . viewer - visualization context 2081 2082 Application Interface Routine: SNESView() 2083 */ 2084 #undef __FUNCT__ 2085 #define __FUNCT__ "SNESView_VI" 2086 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 2087 { 2088 SNES_VI *vi = (SNES_VI *)snes->data; 2089 const char *cstr,*tstr; 2090 PetscErrorCode ierr; 2091 PetscBool iascii; 2092 2093 PetscFunctionBegin; 2094 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2095 if (iascii) { 2096 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2097 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2098 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2099 else if (vi->LineSearch == SNESLineSearchBoundProjectedCubic_VI) cstr = "SNESLineSearchBoundProjectedCubic"; 2100 else cstr = "unknown"; 2101 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2102 else if (snes->ops->solve == SNESSolveVI_AS) tstr = "Active Set"; 2103 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2104 else tstr = "unknown"; 2105 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2106 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2107 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2108 } else { 2109 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ VI",((PetscObject)viewer)->type_name); 2110 } 2111 PetscFunctionReturn(0); 2112 } 2113 2114 /* 2115 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2116 2117 Input Parameters: 2118 . snes - the SNES context. 2119 . xl - lower bound. 2120 . xu - upper bound. 2121 2122 Notes: 2123 If this routine is not called then the lower and upper bounds are set to 2124 PETSC_VI_INF and PETSC_VI_NINF respectively during SNESSetUp(). 2125 2126 */ 2127 2128 #undef __FUNCT__ 2129 #define __FUNCT__ "SNESVISetVariableBounds" 2130 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2131 { 2132 SNES_VI *vi = (SNES_VI*)snes->data; 2133 2134 PetscFunctionBegin; 2135 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2136 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2137 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2138 2139 /* Check if SNESSetFunction is called */ 2140 if(!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2141 2142 /* Check if the vector sizes are compatible for lower and upper bounds */ 2143 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); 2144 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); 2145 vi->usersetxbounds = PETSC_TRUE; 2146 vi->xl = xl; 2147 vi->xu = xu; 2148 2149 PetscFunctionReturn(0); 2150 } 2151 /* -------------------------------------------------------------------------- */ 2152 /* 2153 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2154 2155 Input Parameter: 2156 . snes - the SNES context 2157 2158 Application Interface Routine: SNESSetFromOptions() 2159 */ 2160 #undef __FUNCT__ 2161 #define __FUNCT__ "SNESSetFromOptions_VI" 2162 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2163 { 2164 SNES_VI *vi = (SNES_VI *)snes->data; 2165 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2166 const char *vies[] = {"ss","as","rs"}; 2167 PetscErrorCode ierr; 2168 PetscInt indx; 2169 PetscBool flg,set,flg2; 2170 2171 PetscFunctionBegin; 2172 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2173 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2174 if (flg) { 2175 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2176 } 2177 ierr = PetscOptionsReal("-snes_vi_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2178 ierr = PetscOptionsReal("-snes_vi_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2179 ierr = PetscOptionsReal("-snes_vi_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2180 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2181 ierr = PetscOptionsBool("-snes_vi_lsmonitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2182 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2183 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2184 if (flg2) { 2185 switch (indx) { 2186 case 0: 2187 snes->ops->solve = SNESSolveVI_SS; 2188 break; 2189 case 1: 2190 snes->ops->solve = SNESSolveVI_AS; 2191 break; 2192 case 2: 2193 snes->ops->solve = SNESSolveVI_RS; 2194 break; 2195 } 2196 } 2197 if (snes->ops->solve == SNESSolveVI_RS) { 2198 ierr = SNESLineSearchSet(snes,SNESLineSearchBoundProjectedCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2199 } 2200 else { 2201 ierr = PetscOptionsEList("-snes_vi_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 2202 if (flg) { 2203 switch (indx) { 2204 case 0: 2205 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2206 break; 2207 case 1: 2208 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2209 break; 2210 case 2: 2211 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2212 break; 2213 case 3: 2214 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2215 break; 2216 } 2217 } 2218 } 2219 ierr = PetscOptionsTail();CHKERRQ(ierr); 2220 PetscFunctionReturn(0); 2221 } 2222 /* -------------------------------------------------------------------------- */ 2223 /*MC 2224 SNESVI - Semismooth newton method based nonlinear solver that uses a line search 2225 2226 Options Database: 2227 + -snes_vi [cubic,quadratic,basic,basicnonorms] - Selects line search 2228 . -snes_vi_alpha <alpha> - Sets alpha 2229 . -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) 2230 . -snes_vi_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 2231 - -snes_vi_monitor - print information about progress of line searches 2232 2233 2234 Level: beginner 2235 2236 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2237 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2238 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2239 2240 M*/ 2241 EXTERN_C_BEGIN 2242 #undef __FUNCT__ 2243 #define __FUNCT__ "SNESCreate_VI" 2244 PetscErrorCode SNESCreate_VI(SNES snes) 2245 { 2246 PetscErrorCode ierr; 2247 SNES_VI *vi; 2248 2249 PetscFunctionBegin; 2250 snes->ops->setup = SNESSetUp_VI; 2251 snes->ops->solve = SNESSolveVI_SS; 2252 snes->ops->destroy = SNESDestroy_VI; 2253 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2254 snes->ops->view = SNESView_VI; 2255 snes->ops->converged = SNESDefaultConverged_VI; 2256 2257 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2258 snes->data = (void*)vi; 2259 vi->alpha = 1.e-4; 2260 vi->maxstep = 1.e8; 2261 vi->minlambda = 1.e-12; 2262 vi->LineSearch = SNESLineSearchCubic_VI; 2263 vi->lsP = PETSC_NULL; 2264 vi->postcheckstep = PETSC_NULL; 2265 vi->postcheck = PETSC_NULL; 2266 vi->precheckstep = PETSC_NULL; 2267 vi->precheck = PETSC_NULL; 2268 vi->const_tol = 2.2204460492503131e-16; 2269 vi->computessfunction = ComputeFischerFunction; 2270 2271 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2272 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2273 2274 PetscFunctionReturn(0); 2275 } 2276 EXTERN_C_END 2277