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