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