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