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