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