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