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