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 /* remove later */ 1122 ierr = MatView(snes->jacobian,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1123 1124 /* Create active and inactive index sets */ 1125 1126 /*original 1127 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1128 */ 1129 ierr = SNESVIGetActiveSetIS(snes,X,F,&IS_act);CHKERRQ(ierr); 1130 1131 if (vi->checkredundancy) { 1132 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);CHKERRQ(ierr); 1133 if (IS_redact){ 1134 ierr = ISSort(IS_redact);CHKERRQ(ierr); 1135 ierr = ISComplement(IS_redact,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 1136 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 1137 } 1138 else { 1139 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 1140 } 1141 } else { 1142 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 1143 } 1144 1145 1146 /* Create inactive set submatrix */ 1147 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1148 1149 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 1150 if (0 && keptrows) { 1151 // if (keptrows) { 1152 PetscInt cnt,*nrows,k; 1153 const PetscInt *krows,*inact; 1154 PetscInt rstart=jac_inact_inact->rmap->rstart; 1155 1156 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 1157 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1158 1159 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 1160 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 1161 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 1162 ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); 1163 for (k=0; k<cnt; k++) { 1164 nrows[k] = inact[krows[k]-rstart]; 1165 } 1166 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 1167 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 1168 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 1169 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1170 1171 ierr = ISCreateGeneral(PETSC_COMM_WORLD,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 1172 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 1173 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 1174 } 1175 ierr = DMSetVI(snes->dm,IS_inact);CHKERRQ(ierr); 1176 /* remove later */ 1177 1178 /* 1179 ierr = VecView(vi->xu,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1180 ierr = VecView(vi->xl,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1181 ierr = VecView(X,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1182 ierr = VecView(F,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1183 ierr = ISView(IS_inact,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1184 */ 1185 1186 /* Get sizes of active and inactive sets */ 1187 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1188 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 1189 1190 /* Create active and inactive set vectors */ 1191 ierr = SNESVICreateSubVectors(snes,nis_inact,&F_inact);CHKERRQ(ierr); 1192 ierr = SNESVICreateSubVectors(snes,nis_act,&Y_act);CHKERRQ(ierr); 1193 ierr = SNESVICreateSubVectors(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 1194 1195 /* Create scatter contexts */ 1196 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 1197 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 1198 1199 /* Do a vec scatter to active and inactive set vectors */ 1200 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1201 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1202 1203 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1204 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1205 1206 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1207 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1208 1209 /* Active set direction = 0 */ 1210 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 1211 if (snes->jacobian != snes->jacobian_pre) { 1212 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 1213 } else prejac_inact_inact = jac_inact_inact; 1214 1215 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1216 if (!isequal) { 1217 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 1218 flg = DIFFERENT_NONZERO_PATTERN; 1219 } 1220 1221 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 1222 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 1223 /* ierr = MatView(snes->jacobian_pre,0); */ 1224 1225 /* remove later */ 1226 ierr = MatView(jac_inact_inact,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1227 ierr = ISView(IS_inact,PETSC_VIEWER_BINARY_(PETSC_COMM_WORLD));CHKERRQ(ierr); 1228 1229 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 1230 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1231 { 1232 PC pc; 1233 PetscBool flg; 1234 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1235 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1236 if (flg) { 1237 KSP *subksps; 1238 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1239 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1240 ierr = PetscFree(subksps);CHKERRQ(ierr); 1241 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1242 if (flg) { 1243 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1244 const PetscInt *ii; 1245 1246 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1247 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1248 for (j=0; j<n; j++) { 1249 if (ii[j] < N) cnts[0]++; 1250 else if (ii[j] < 2*N) cnts[1]++; 1251 else if (ii[j] < 3*N) cnts[2]++; 1252 } 1253 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1254 1255 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1256 } 1257 } 1258 } 1259 1260 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 1261 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1262 if (kspreason < 0) { 1263 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1264 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1265 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1266 break; 1267 } 1268 } 1269 1270 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1271 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1272 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1273 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1274 1275 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 1276 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 1277 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 1278 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 1279 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 1280 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1281 if (!isequal) { 1282 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1283 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1284 } 1285 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1286 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 1287 if (snes->jacobian != snes->jacobian_pre) { 1288 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 1289 } 1290 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1291 snes->linear_its += lits; 1292 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1293 /* 1294 if (vi->precheckstep) { 1295 PetscBool changed_y = PETSC_FALSE; 1296 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1297 } 1298 1299 if (PetscLogPrintInfo){ 1300 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1301 } 1302 */ 1303 /* Compute a (scaled) negative update in the line search routine: 1304 Y <- X - lambda*Y 1305 and evaluate G = function(Y) (depends on the line search). 1306 */ 1307 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1308 ynorm = 1; gnorm = fnorm; 1309 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1310 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1311 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1312 if (snes->domainerror) { 1313 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1314 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 1315 PetscFunctionReturn(0); 1316 } 1317 if (!lssucceed) { 1318 if (++snes->numFailures >= snes->maxFailures) { 1319 PetscBool ismin; 1320 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1321 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1322 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1323 break; 1324 } 1325 } 1326 /* Update function and solution vectors */ 1327 fnorm = gnorm; 1328 ierr = VecCopy(G,F);CHKERRQ(ierr); 1329 ierr = VecCopy(W,X);CHKERRQ(ierr); 1330 /* Monitor convergence */ 1331 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1332 snes->iter = i+1; 1333 snes->norm = fnorm; 1334 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1335 SNESLogConvHistory(snes,snes->norm,lits); 1336 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1337 /* Test for convergence, xnorm = || X || */ 1338 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1339 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1340 if (snes->reason) break; 1341 } 1342 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 1343 if (i == maxits) { 1344 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1345 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1346 } 1347 PetscFunctionReturn(0); 1348 } 1349 1350 #undef __FUNCT__ 1351 #define __FUNCT__ "SNESVISetRedundancyCheck" 1352 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 1353 { 1354 SNES_VI *vi = (SNES_VI*)snes->data; 1355 1356 PetscFunctionBegin; 1357 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1358 vi->checkredundancy = func; 1359 vi->ctxP = ctx; 1360 PetscFunctionReturn(0); 1361 } 1362 1363 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1364 #include <engine.h> 1365 #include <mex.h> 1366 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 1367 1368 #undef __FUNCT__ 1369 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 1370 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 1371 { 1372 PetscErrorCode ierr; 1373 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 1374 int nlhs = 1, nrhs = 5; 1375 mxArray *plhs[1], *prhs[5]; 1376 long long int l1 = 0, l2 = 0, ls = 0; 1377 PetscInt *indices=PETSC_NULL; 1378 1379 PetscFunctionBegin; 1380 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1381 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 1382 PetscValidPointer(is_redact,3); 1383 PetscCheckSameComm(snes,1,is_act,2); 1384 1385 /* Create IS for reduced active set of size 0, its size and indices will 1386 bet set by the Matlab function */ 1387 ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 1388 /* call Matlab function in ctx */ 1389 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 1390 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 1391 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 1392 prhs[0] = mxCreateDoubleScalar((double)ls); 1393 prhs[1] = mxCreateDoubleScalar((double)l1); 1394 prhs[2] = mxCreateDoubleScalar((double)l2); 1395 prhs[3] = mxCreateString(sctx->funcname); 1396 prhs[4] = sctx->ctx; 1397 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 1398 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 1399 mxDestroyArray(prhs[0]); 1400 mxDestroyArray(prhs[1]); 1401 mxDestroyArray(prhs[2]); 1402 mxDestroyArray(prhs[3]); 1403 mxDestroyArray(plhs[0]); 1404 PetscFunctionReturn(0); 1405 } 1406 1407 #undef __FUNCT__ 1408 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 1409 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 1410 { 1411 PetscErrorCode ierr; 1412 SNESMatlabContext *sctx; 1413 1414 PetscFunctionBegin; 1415 /* currently sctx is memory bleed */ 1416 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 1417 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 1418 sctx->ctx = mxDuplicateArray(ctx); 1419 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 #endif 1424 1425 1426 /* Variational Inequality solver using augmented space method. It does the opposite of the 1427 reduced space method i.e. it identifies the active set variables and instead of discarding 1428 them it augments the original system by introducing additional equality 1429 constraint equations for active set variables. The user can optionally provide an IS for 1430 a subset of 'redundant' active set variables which will treated as free variables. 1431 Specific implementation for Allen-Cahn problem 1432 */ 1433 #undef __FUNCT__ 1434 #define __FUNCT__ "SNESSolveVI_RSAUG" 1435 PetscErrorCode SNESSolveVI_RSAUG(SNES snes) 1436 { 1437 SNES_VI *vi = (SNES_VI*)snes->data; 1438 PetscErrorCode ierr; 1439 PetscInt maxits,i,lits; 1440 PetscBool lssucceed; 1441 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 1442 PetscReal fnorm,gnorm,xnorm=0,ynorm; 1443 Vec Y,X,F,G,W; 1444 KSPConvergedReason kspreason; 1445 1446 PetscFunctionBegin; 1447 snes->numFailures = 0; 1448 snes->numLinearSolveFailures = 0; 1449 snes->reason = SNES_CONVERGED_ITERATING; 1450 1451 maxits = snes->max_its; /* maximum number of iterations */ 1452 X = snes->vec_sol; /* solution vector */ 1453 F = snes->vec_func; /* residual vector */ 1454 Y = snes->work[0]; /* work vectors */ 1455 G = snes->work[1]; 1456 W = snes->work[2]; 1457 1458 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1459 snes->iter = 0; 1460 snes->norm = 0.0; 1461 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1462 1463 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 1464 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 1465 if (snes->domainerror) { 1466 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1467 PetscFunctionReturn(0); 1468 } 1469 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 1470 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 1471 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 1472 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1473 1474 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1475 snes->norm = fnorm; 1476 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1477 SNESLogConvHistory(snes,fnorm,0); 1478 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 1479 1480 /* set parameter for default relative tolerance convergence test */ 1481 snes->ttol = fnorm*snes->rtol; 1482 /* test convergence */ 1483 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1484 if (snes->reason) PetscFunctionReturn(0); 1485 1486 for (i=0; i<maxits; i++) { 1487 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 1488 IS IS_redact; /* redundant active set */ 1489 Mat J_aug,Jpre_aug; 1490 Vec F_aug,Y_aug; 1491 PetscInt nis_redact,nis_act; 1492 const PetscInt *idx_redact,*idx_act; 1493 PetscInt k; 1494 PetscInt *idx_actkept=PETSC_NULL,nkept=0; /* list of kept active set */ 1495 PetscScalar *f,*f2; 1496 PetscBool isequal; 1497 PetscInt i1=0,j1=0; 1498 1499 /* Call general purpose update function */ 1500 if (snes->ops->update) { 1501 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 1502 } 1503 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 1504 1505 /* Create active and inactive index sets */ 1506 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 1507 1508 /* Get local active set size */ 1509 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 1510 if (nis_act) { 1511 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1512 IS_redact = PETSC_NULL; 1513 if(vi->checkredundancy) { 1514 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP); 1515 } 1516 1517 if(!IS_redact) { 1518 /* User called checkredundancy function but didn't create IS_redact because 1519 there were no redundant active set variables */ 1520 /* Copy over all active set indices to the list */ 1521 ierr = PetscMalloc(nis_act*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1522 for(k=0;k < nis_act;k++) idx_actkept[k] = idx_act[k]; 1523 nkept = nis_act; 1524 } else { 1525 ierr = ISGetLocalSize(IS_redact,&nis_redact);CHKERRQ(ierr); 1526 ierr = PetscMalloc((nis_act-nis_redact)*sizeof(PetscInt),&idx_actkept);CHKERRQ(ierr); 1527 1528 /* Create reduced active set list */ 1529 ierr = ISGetIndices(IS_act,&idx_act);CHKERRQ(ierr); 1530 ierr = ISGetIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1531 j1 = 0;nkept = 0; 1532 for(k=0;k<nis_act;k++) { 1533 if(j1 < nis_redact && idx_act[k] == idx_redact[j1]) j1++; 1534 else idx_actkept[nkept++] = idx_act[k]; 1535 } 1536 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1537 ierr = ISRestoreIndices(IS_redact,&idx_redact);CHKERRQ(ierr); 1538 1539 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 1540 } 1541 1542 /* Create augmented F and Y */ 1543 ierr = VecCreate(((PetscObject)snes)->comm,&F_aug);CHKERRQ(ierr); 1544 ierr = VecSetSizes(F_aug,F->map->n+nkept,PETSC_DECIDE);CHKERRQ(ierr); 1545 ierr = VecSetFromOptions(F_aug);CHKERRQ(ierr); 1546 ierr = VecDuplicate(F_aug,&Y_aug);CHKERRQ(ierr); 1547 1548 /* Copy over F to F_aug in the first n locations */ 1549 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 1550 ierr = VecGetArray(F_aug,&f2);CHKERRQ(ierr); 1551 ierr = PetscMemcpy(f2,f,F->map->n*sizeof(PetscScalar));CHKERRQ(ierr); 1552 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 1553 ierr = VecRestoreArray(F_aug,&f2);CHKERRQ(ierr); 1554 1555 /* Create the augmented jacobian matrix */ 1556 ierr = MatCreate(((PetscObject)snes)->comm,&J_aug);CHKERRQ(ierr); 1557 ierr = MatSetSizes(J_aug,X->map->n+nkept,X->map->n+nkept,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1558 ierr = MatSetFromOptions(J_aug);CHKERRQ(ierr); 1559 1560 1561 { /* local vars */ 1562 /* Preallocate augmented matrix and set values in it...Doing only sequential case first*/ 1563 PetscInt ncols; 1564 const PetscInt *cols; 1565 const PetscScalar *vals; 1566 PetscScalar value[2]; 1567 PetscInt row,col[2]; 1568 PetscInt *d_nnz; 1569 value[0] = 1.0; value[1] = 0.0; 1570 ierr = PetscMalloc((X->map->n+nkept)*sizeof(PetscInt),&d_nnz);CHKERRQ(ierr); 1571 ierr = PetscMemzero(d_nnz,(X->map->n+nkept)*sizeof(PetscInt));CHKERRQ(ierr); 1572 for(row=0;row<snes->jacobian->rmap->n;row++) { 1573 ierr = MatGetRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1574 d_nnz[row] += ncols; 1575 ierr = MatRestoreRow(snes->jacobian,row,&ncols,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 1576 } 1577 1578 for(k=0;k<nkept;k++) { 1579 d_nnz[idx_actkept[k]] += 1; 1580 d_nnz[snes->jacobian->rmap->n+k] += 2; 1581 } 1582 ierr = MatSeqAIJSetPreallocation(J_aug,PETSC_NULL,d_nnz);CHKERRQ(ierr); 1583 1584 ierr = PetscFree(d_nnz);CHKERRQ(ierr); 1585 1586 /* Set the original jacobian matrix in J_aug */ 1587 for(row=0;row<snes->jacobian->rmap->n;row++) { 1588 ierr = MatGetRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1589 ierr = MatSetValues(J_aug,1,&row,ncols,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 1590 ierr = MatRestoreRow(snes->jacobian,row,&ncols,&cols,&vals);CHKERRQ(ierr); 1591 } 1592 /* Add the augmented part */ 1593 for(k=0;k<nkept;k++) { 1594 row = snes->jacobian->rmap->n+k; 1595 col[0] = idx_actkept[k]; col[1] = snes->jacobian->rmap->n+k; 1596 ierr = MatSetValues(J_aug,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 1597 ierr = MatSetValues(J_aug,1,&col[0],1,&row,&value[0],INSERT_VALUES);CHKERRQ(ierr); 1598 } 1599 ierr = MatAssemblyBegin(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1600 ierr = MatAssemblyEnd(J_aug,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1601 /* Only considering prejac = jac for now */ 1602 Jpre_aug = J_aug; 1603 } /* local vars*/ 1604 ierr = ISRestoreIndices(IS_act,&idx_act);CHKERRQ(ierr); 1605 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 1606 } else { 1607 F_aug = F; J_aug = snes->jacobian; Y_aug = Y; Jpre_aug = snes->jacobian_pre; 1608 } 1609 1610 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 1611 if (!isequal) { 1612 ierr = SNESVIResetPCandKSP(snes,J_aug,Jpre_aug);CHKERRQ(ierr); 1613 } 1614 ierr = KSPSetOperators(snes->ksp,J_aug,Jpre_aug,flg);CHKERRQ(ierr); 1615 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 1616 /* { 1617 PC pc; 1618 PetscBool flg; 1619 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 1620 ierr = PetscTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 1621 if (flg) { 1622 KSP *subksps; 1623 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 1624 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 1625 ierr = PetscFree(subksps);CHKERRQ(ierr); 1626 ierr = PetscTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 1627 if (flg) { 1628 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 1629 const PetscInt *ii; 1630 1631 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 1632 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 1633 for (j=0; j<n; j++) { 1634 if (ii[j] < N) cnts[0]++; 1635 else if (ii[j] < 2*N) cnts[1]++; 1636 else if (ii[j] < 3*N) cnts[2]++; 1637 } 1638 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 1639 1640 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 1641 } 1642 } 1643 } 1644 */ 1645 ierr = SNES_KSPSolve(snes,snes->ksp,F_aug,Y_aug);CHKERRQ(ierr); 1646 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 1647 if (kspreason < 0) { 1648 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 1649 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 1650 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 1651 break; 1652 } 1653 } 1654 1655 if(nis_act) { 1656 PetscScalar *y1,*y2; 1657 ierr = VecGetArray(Y,&y1);CHKERRQ(ierr); 1658 ierr = VecGetArray(Y_aug,&y2);CHKERRQ(ierr); 1659 /* Copy over inactive Y_aug to Y */ 1660 j1 = 0; 1661 for(i1=Y->map->rstart;i1 < Y->map->rend;i1++) { 1662 if(j1 < nkept && idx_actkept[j1] == i1) j1++; 1663 else y1[i1-Y->map->rstart] = y2[i1-Y->map->rstart]; 1664 } 1665 ierr = VecRestoreArray(Y,&y1);CHKERRQ(ierr); 1666 ierr = VecRestoreArray(Y_aug,&y2);CHKERRQ(ierr); 1667 1668 ierr = VecDestroy(&F_aug);CHKERRQ(ierr); 1669 ierr = VecDestroy(&Y_aug);CHKERRQ(ierr); 1670 ierr = MatDestroy(&J_aug);CHKERRQ(ierr); 1671 ierr = PetscFree(idx_actkept);CHKERRQ(ierr); 1672 } 1673 1674 if (!isequal) { 1675 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1676 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 1677 } 1678 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 1679 1680 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 1681 snes->linear_its += lits; 1682 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 1683 /* 1684 if (vi->precheckstep) { 1685 PetscBool changed_y = PETSC_FALSE; 1686 ierr = (*vi->precheckstep)(snes,X,Y,vi->precheck,&changed_y);CHKERRQ(ierr); 1687 } 1688 1689 if (PetscLogPrintInfo){ 1690 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 1691 } 1692 */ 1693 /* Compute a (scaled) negative update in the line search routine: 1694 Y <- X - lambda*Y 1695 and evaluate G = function(Y) (depends on the line search). 1696 */ 1697 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 1698 ynorm = 1; gnorm = fnorm; 1699 ierr = (*vi->LineSearch)(snes,vi->lsP,X,F,G,Y,W,fnorm,xnorm,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr); 1700 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",fnorm,gnorm,ynorm,(int)lssucceed);CHKERRQ(ierr); 1701 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 1702 if (snes->domainerror) { 1703 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 1704 PetscFunctionReturn(0); 1705 } 1706 if (!lssucceed) { 1707 if (++snes->numFailures >= snes->maxFailures) { 1708 PetscBool ismin; 1709 snes->reason = SNES_DIVERGED_LINE_SEARCH; 1710 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 1711 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 1712 break; 1713 } 1714 } 1715 /* Update function and solution vectors */ 1716 fnorm = gnorm; 1717 ierr = VecCopy(G,F);CHKERRQ(ierr); 1718 ierr = VecCopy(W,X);CHKERRQ(ierr); 1719 /* Monitor convergence */ 1720 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 1721 snes->iter = i+1; 1722 snes->norm = fnorm; 1723 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 1724 SNESLogConvHistory(snes,snes->norm,lits); 1725 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 1726 /* Test for convergence, xnorm = || X || */ 1727 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 1728 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 1729 if (snes->reason) break; 1730 } 1731 if (i == maxits) { 1732 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 1733 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 1734 } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 /* -------------------------------------------------------------------------- */ 1739 /* 1740 SNESSetUp_VI - Sets up the internal data structures for the later use 1741 of the SNESVI nonlinear solver. 1742 1743 Input Parameter: 1744 . snes - the SNES context 1745 . x - the solution vector 1746 1747 Application Interface Routine: SNESSetUp() 1748 1749 Notes: 1750 For basic use of the SNES solvers, the user need not explicitly call 1751 SNESSetUp(), since these actions will automatically occur during 1752 the call to SNESSolve(). 1753 */ 1754 #undef __FUNCT__ 1755 #define __FUNCT__ "SNESSetUp_VI" 1756 PetscErrorCode SNESSetUp_VI(SNES snes) 1757 { 1758 PetscErrorCode ierr; 1759 SNES_VI *vi = (SNES_VI*) snes->data; 1760 PetscInt i_start[3],i_end[3]; 1761 1762 PetscFunctionBegin; 1763 1764 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 1765 1766 if (vi->computevariablebounds) { 1767 if (!vi->xl) {ierr = VecDuplicate(snes->vec_sol,&vi->xl);CHKERRQ(ierr);} 1768 if (!vi->xu) {ierr = VecDuplicate(snes->vec_sol,&vi->xu);CHKERRQ(ierr);} 1769 ierr = (*vi->computevariablebounds)(snes,vi->xl,vi->xu);CHKERRQ(ierr); 1770 } else if (!vi->xl && !vi->xu) { 1771 /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */ 1772 ierr = VecDuplicate(snes->vec_sol, &vi->xl);CHKERRQ(ierr); 1773 ierr = VecSet(vi->xl,SNES_VI_NINF);CHKERRQ(ierr); 1774 ierr = VecDuplicate(snes->vec_sol, &vi->xu);CHKERRQ(ierr); 1775 ierr = VecSet(vi->xu,SNES_VI_INF);CHKERRQ(ierr); 1776 } else { 1777 /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */ 1778 ierr = VecGetOwnershipRange(snes->vec_sol,i_start,i_end);CHKERRQ(ierr); 1779 ierr = VecGetOwnershipRange(vi->xl,i_start+1,i_end+1);CHKERRQ(ierr); 1780 ierr = VecGetOwnershipRange(vi->xu,i_start+2,i_end+2);CHKERRQ(ierr); 1781 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])) 1782 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."); 1783 } 1784 if (snes->ops->solve != SNESSolveVI_SS) { 1785 /* Set up previous active index set for the first snes solve 1786 vi->IS_inact_prev = 0,1,2,....N */ 1787 PetscInt *indices; 1788 PetscInt i,n,rstart,rend; 1789 1790 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 1791 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 1792 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 1793 for(i=0;i < n; i++) indices[i] = rstart + i; 1794 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev); 1795 } 1796 1797 if (snes->ops->solve == SNESSolveVI_SS) { 1798 ierr = VecDuplicate(snes->vec_sol, &vi->dpsi);CHKERRQ(ierr); 1799 ierr = VecDuplicate(snes->vec_sol, &vi->phi);CHKERRQ(ierr); 1800 ierr = VecDuplicate(snes->vec_sol, &vi->Da);CHKERRQ(ierr); 1801 ierr = VecDuplicate(snes->vec_sol, &vi->Db);CHKERRQ(ierr); 1802 ierr = VecDuplicate(snes->vec_sol, &vi->z);CHKERRQ(ierr); 1803 ierr = VecDuplicate(snes->vec_sol, &vi->t);CHKERRQ(ierr); 1804 } 1805 PetscFunctionReturn(0); 1806 } 1807 /* -------------------------------------------------------------------------- */ 1808 #undef __FUNCT__ 1809 #define __FUNCT__ "SNESReset_VI" 1810 PetscErrorCode SNESReset_VI(SNES snes) 1811 { 1812 SNES_VI *vi = (SNES_VI*) snes->data; 1813 PetscErrorCode ierr; 1814 1815 PetscFunctionBegin; 1816 ierr = VecDestroy(&vi->xl);CHKERRQ(ierr); 1817 ierr = VecDestroy(&vi->xu);CHKERRQ(ierr); 1818 if (snes->ops->solve != SNESSolveVI_SS) { 1819 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 1820 } 1821 PetscFunctionReturn(0); 1822 } 1823 1824 /* 1825 SNESDestroy_VI - Destroys the private SNES_VI context that was created 1826 with SNESCreate_VI(). 1827 1828 Input Parameter: 1829 . snes - the SNES context 1830 1831 Application Interface Routine: SNESDestroy() 1832 */ 1833 #undef __FUNCT__ 1834 #define __FUNCT__ "SNESDestroy_VI" 1835 PetscErrorCode SNESDestroy_VI(SNES snes) 1836 { 1837 SNES_VI *vi = (SNES_VI*) snes->data; 1838 PetscErrorCode ierr; 1839 1840 PetscFunctionBegin; 1841 1842 if (snes->ops->solve == SNESSolveVI_SS) { 1843 /* clear vectors */ 1844 ierr = VecDestroy(&vi->dpsi);CHKERRQ(ierr); 1845 ierr = VecDestroy(&vi->phi);CHKERRQ(ierr); 1846 ierr = VecDestroy(&vi->Da);CHKERRQ(ierr); 1847 ierr = VecDestroy(&vi->Db);CHKERRQ(ierr); 1848 ierr = VecDestroy(&vi->z);CHKERRQ(ierr); 1849 ierr = VecDestroy(&vi->t);CHKERRQ(ierr); 1850 } 1851 1852 ierr = PetscViewerDestroy(&vi->lsmonitor);CHKERRQ(ierr); 1853 ierr = PetscFree(snes->data);CHKERRQ(ierr); 1854 1855 /* clear composed functions */ 1856 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","",PETSC_NULL);CHKERRQ(ierr); 1857 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","",PETSC_NULL);CHKERRQ(ierr); 1858 PetscFunctionReturn(0); 1859 } 1860 1861 /* -------------------------------------------------------------------------- */ 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "SNESLineSearchNo_VI" 1864 1865 /* 1866 This routine does not actually do a line search but it takes a full newton 1867 step while ensuring that the new iterates remain within the constraints. 1868 1869 */ 1870 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) 1871 { 1872 PetscErrorCode ierr; 1873 SNES_VI *vi = (SNES_VI*)snes->data; 1874 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1875 1876 PetscFunctionBegin; 1877 *flag = PETSC_TRUE; 1878 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1879 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); /* ynorm = || y || */ 1880 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1881 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1882 if (vi->postcheckstep) { 1883 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1884 } 1885 if (changed_y) { 1886 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1887 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1888 } 1889 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1890 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1891 if (!snes->domainerror) { 1892 if (snes->ops->solve != SNESSolveVI_SS) { 1893 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 1894 } else { 1895 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); /* gnorm = || g || */ 1896 } 1897 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 1898 } 1899 if (vi->lsmonitor) { 1900 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1901 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 1902 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1903 } 1904 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1905 PetscFunctionReturn(0); 1906 } 1907 1908 /* -------------------------------------------------------------------------- */ 1909 #undef __FUNCT__ 1910 #define __FUNCT__ "SNESLineSearchNoNorms_VI" 1911 1912 /* 1913 This routine is a copy of SNESLineSearchNoNorms in snes/impls/ls/ls.c 1914 */ 1915 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) 1916 { 1917 PetscErrorCode ierr; 1918 SNES_VI *vi = (SNES_VI*)snes->data; 1919 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1920 1921 PetscFunctionBegin; 1922 *flag = PETSC_TRUE; 1923 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1924 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1925 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1926 if (vi->postcheckstep) { 1927 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 1928 } 1929 if (changed_y) { 1930 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); /* w <- x - y */ 1931 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 1932 } 1933 1934 /* don't evaluate function the last time through */ 1935 if (snes->iter < snes->max_its-1) { 1936 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 1937 } 1938 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1939 PetscFunctionReturn(0); 1940 } 1941 1942 /* -------------------------------------------------------------------------- */ 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "SNESLineSearchCubic_VI" 1945 /* 1946 This routine implements a cubic line search while doing a projection on the variable bounds 1947 */ 1948 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) 1949 { 1950 PetscReal initslope,lambdaprev,gnormprev,a,b,d,t1,t2,rellength; 1951 PetscReal minlambda,lambda,lambdatemp; 1952 #if defined(PETSC_USE_COMPLEX) 1953 PetscScalar cinitslope; 1954 #endif 1955 PetscErrorCode ierr; 1956 PetscInt count; 1957 SNES_VI *vi = (SNES_VI*)snes->data; 1958 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 1959 MPI_Comm comm; 1960 1961 PetscFunctionBegin; 1962 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 1963 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 1964 *flag = PETSC_TRUE; 1965 1966 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 1967 if (*ynorm == 0.0) { 1968 if (vi->lsmonitor) { 1969 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1970 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Initial direction and size is 0\n");CHKERRQ(ierr); 1971 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1972 } 1973 *gnorm = fnorm; 1974 ierr = VecCopy(x,w);CHKERRQ(ierr); 1975 ierr = VecCopy(f,g);CHKERRQ(ierr); 1976 *flag = PETSC_FALSE; 1977 goto theend1; 1978 } 1979 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 1980 if (vi->lsmonitor) { 1981 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1982 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Scaling step by %G old ynorm %G\n",vi->maxstep/(*ynorm),*ynorm);CHKERRQ(ierr); 1983 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 1984 } 1985 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 1986 *ynorm = vi->maxstep; 1987 } 1988 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 1989 minlambda = vi->minlambda/rellength; 1990 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 1991 #if defined(PETSC_USE_COMPLEX) 1992 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 1993 initslope = PetscRealPart(cinitslope); 1994 #else 1995 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 1996 #endif 1997 if (initslope > 0.0) initslope = -initslope; 1998 if (initslope == 0.0) initslope = -1.0; 1999 2000 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2001 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2002 if (snes->nfuncs >= snes->max_funcs) { 2003 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 2004 *flag = PETSC_FALSE; 2005 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2006 goto theend1; 2007 } 2008 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2009 if (snes->ops->solve != SNESSolveVI_SS) { 2010 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2011 } else { 2012 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2013 } 2014 if (snes->domainerror) { 2015 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2016 PetscFunctionReturn(0); 2017 } 2018 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2019 ierr = PetscInfo4(snes,"Initial fnorm %G gnorm %G alpha %G initslope %G\n",fnorm,*gnorm,vi->alpha,initslope);CHKERRQ(ierr); 2020 if ((*gnorm)*(*gnorm) <= (1.0 - vi->alpha)*fnorm*fnorm ) { /* Sufficient reduction */ 2021 if (vi->lsmonitor) { 2022 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2023 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 2024 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2025 } 2026 goto theend1; 2027 } 2028 2029 /* Fit points with quadratic */ 2030 lambda = 1.0; 2031 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 2032 lambdaprev = lambda; 2033 gnormprev = *gnorm; 2034 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2035 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2036 else lambda = lambdatemp; 2037 2038 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2039 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2040 if (snes->nfuncs >= snes->max_funcs) { 2041 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while attempting quadratic backtracking! %D \n",snes->nfuncs);CHKERRQ(ierr); 2042 *flag = PETSC_FALSE; 2043 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2044 goto theend1; 2045 } 2046 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2047 if (snes->ops->solve != SNESSolveVI_SS) { 2048 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2049 } else { 2050 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2051 } 2052 if (snes->domainerror) { 2053 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2057 if (vi->lsmonitor) { 2058 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2059 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: gnorm after quadratic fit %G\n",*gnorm);CHKERRQ(ierr); 2060 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2061 } 2062 if ((*gnorm)*(*gnorm) < (1.0 - vi->alpha)*fnorm*fnorm ) { /* sufficient reduction */ 2063 if (vi->lsmonitor) { 2064 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2065 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Quadratically determined step, lambda=%18.16e\n",lambda);CHKERRQ(ierr); 2066 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2067 } 2068 goto theend1; 2069 } 2070 2071 /* Fit points with cubic */ 2072 count = 1; 2073 while (PETSC_TRUE) { 2074 if (lambda <= minlambda) { 2075 if (vi->lsmonitor) { 2076 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2077 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: unable to find good step length! After %D tries \n",count);CHKERRQ(ierr); 2078 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); 2079 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2080 } 2081 *flag = PETSC_FALSE; 2082 break; 2083 } 2084 t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope; 2085 t2 = .5*(gnormprev*gnormprev - fnorm*fnorm) - lambdaprev*initslope; 2086 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 2087 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 2088 d = b*b - 3*a*initslope; 2089 if (d < 0.0) d = 0.0; 2090 if (a == 0.0) { 2091 lambdatemp = -initslope/(2.0*b); 2092 } else { 2093 lambdatemp = (-b + sqrt(d))/(3.0*a); 2094 } 2095 lambdaprev = lambda; 2096 gnormprev = *gnorm; 2097 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2098 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2099 else lambda = lambdatemp; 2100 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2101 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2102 if (snes->nfuncs >= snes->max_funcs) { 2103 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 2104 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); 2105 *flag = PETSC_FALSE; 2106 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2107 break; 2108 } 2109 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2110 if (snes->ops->solve != SNESSolveVI_SS) { 2111 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2112 } else { 2113 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2114 } 2115 if (snes->domainerror) { 2116 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2120 if ((*gnorm)*(*gnorm) < (1.0 - vi->alpha)*fnorm*fnorm) { /* is reduction enough? */ 2121 if (vi->lsmonitor) { 2122 ierr = PetscPrintf(comm," Line search: Cubically determined step, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 2123 } 2124 break; 2125 } else { 2126 if (vi->lsmonitor) { 2127 ierr = PetscPrintf(comm," Line search: Cubic step no good, shrinking lambda, current gnorm %G lambda=%18.16e\n",*gnorm,lambda);CHKERRQ(ierr); 2128 } 2129 } 2130 count++; 2131 } 2132 theend1: 2133 /* Optional user-defined check for line search step validity */ 2134 if (vi->postcheckstep && *flag) { 2135 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2136 if (changed_y) { 2137 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2138 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2139 } 2140 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2141 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2142 if (snes->ops->solve != SNESSolveVI_SS) { 2143 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2144 } else { 2145 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2146 } 2147 if (snes->domainerror) { 2148 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2149 PetscFunctionReturn(0); 2150 } 2151 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2152 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2153 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2154 } 2155 } 2156 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2157 PetscFunctionReturn(0); 2158 } 2159 2160 #undef __FUNCT__ 2161 #define __FUNCT__ "SNESLineSearchQuadratic_VI" 2162 /* 2163 This routine does a quadratic line search while keeping the iterates within the variable bounds 2164 */ 2165 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) 2166 { 2167 /* 2168 Note that for line search purposes we work with with the related 2169 minimization problem: 2170 min z(x): R^n -> R, 2171 where z(x) = .5 * fnorm*fnorm,and fnorm = || f ||_2. 2172 */ 2173 PetscReal initslope,minlambda,lambda,lambdatemp,rellength; 2174 #if defined(PETSC_USE_COMPLEX) 2175 PetscScalar cinitslope; 2176 #endif 2177 PetscErrorCode ierr; 2178 PetscInt count; 2179 SNES_VI *vi = (SNES_VI*)snes->data; 2180 PetscBool changed_w = PETSC_FALSE,changed_y = PETSC_FALSE; 2181 2182 PetscFunctionBegin; 2183 ierr = PetscLogEventBegin(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2184 *flag = PETSC_TRUE; 2185 2186 ierr = VecNorm(y,NORM_2,ynorm);CHKERRQ(ierr); 2187 if (*ynorm == 0.0) { 2188 if (vi->lsmonitor) { 2189 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2190 ierr = PetscViewerASCIIPrintf(vi->lsmonitor,"Line search: Direction and size is 0\n");CHKERRQ(ierr); 2191 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2192 } 2193 *gnorm = fnorm; 2194 ierr = VecCopy(x,w);CHKERRQ(ierr); 2195 ierr = VecCopy(f,g);CHKERRQ(ierr); 2196 *flag = PETSC_FALSE; 2197 goto theend2; 2198 } 2199 if (*ynorm > vi->maxstep) { /* Step too big, so scale back */ 2200 ierr = VecScale(y,vi->maxstep/(*ynorm));CHKERRQ(ierr); 2201 *ynorm = vi->maxstep; 2202 } 2203 ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr); 2204 minlambda = vi->minlambda/rellength; 2205 ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr); 2206 #if defined(PETSC_USE_COMPLEX) 2207 ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr); 2208 initslope = PetscRealPart(cinitslope); 2209 #else 2210 ierr = VecDot(f,w,&initslope);CHKERRQ(ierr); 2211 #endif 2212 if (initslope > 0.0) initslope = -initslope; 2213 if (initslope == 0.0) initslope = -1.0; 2214 ierr = PetscInfo1(snes,"Initslope %G \n",initslope);CHKERRQ(ierr); 2215 2216 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2217 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2218 if (snes->nfuncs >= snes->max_funcs) { 2219 ierr = PetscInfo(snes,"Exceeded maximum function evaluations, while checking full step length!\n");CHKERRQ(ierr); 2220 *flag = PETSC_FALSE; 2221 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2222 goto theend2; 2223 } 2224 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2225 if (snes->ops->solve != SNESSolveVI_SS) { 2226 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2227 } else { 2228 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2229 } 2230 if (snes->domainerror) { 2231 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2232 PetscFunctionReturn(0); 2233 } 2234 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2235 if ((*gnorm)*(*gnorm) <= (1.0 - vi->alpha)*fnorm*fnorm) { /* Sufficient reduction */ 2236 if (vi->lsmonitor) { 2237 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2238 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line search: Using full step: fnorm %G gnorm %G\n",fnorm,*gnorm);CHKERRQ(ierr); 2239 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2240 } 2241 goto theend2; 2242 } 2243 2244 /* Fit points with quadratic */ 2245 lambda = 1.0; 2246 count = 1; 2247 while (PETSC_TRUE) { 2248 if (lambda <= minlambda) { /* bad luck; use full step */ 2249 if (vi->lsmonitor) { 2250 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2251 ierr = PetscViewerASCIIPrintf(vi->lsmonitor,"Line search: Unable to find good step length! %D \n",count);CHKERRQ(ierr); 2252 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); 2253 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2254 } 2255 ierr = VecCopy(x,w);CHKERRQ(ierr); 2256 *flag = PETSC_FALSE; 2257 break; 2258 } 2259 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 2260 if (lambdatemp > .5*lambda) lambdatemp = .5*lambda; 2261 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 2262 else lambda = lambdatemp; 2263 2264 ierr = VecWAXPY(w,-lambda,y,x);CHKERRQ(ierr); 2265 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2266 if (snes->nfuncs >= snes->max_funcs) { 2267 ierr = PetscInfo1(snes,"Exceeded maximum function evaluations, while looking for good step length! %D \n",count);CHKERRQ(ierr); 2268 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); 2269 *flag = PETSC_FALSE; 2270 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 2271 break; 2272 } 2273 ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr); 2274 if (snes->domainerror) { 2275 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2276 PetscFunctionReturn(0); 2277 } 2278 if (snes->ops->solve != SNESSolveVI_SS) { 2279 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2280 } else { 2281 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2282 } 2283 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2284 if ((*gnorm)*(*gnorm) < (1.0 - vi->alpha)*fnorm*fnorm) { /* sufficient reduction */ 2285 if (vi->lsmonitor) { 2286 ierr = PetscViewerASCIIAddTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2287 ierr = PetscViewerASCIIPrintf(vi->lsmonitor," Line Search: Quadratically determined step, lambda=%G\n",lambda);CHKERRQ(ierr); 2288 ierr = PetscViewerASCIISubtractTab(vi->lsmonitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr); 2289 } 2290 break; 2291 } 2292 count++; 2293 } 2294 theend2: 2295 /* Optional user-defined check for line search step validity */ 2296 if (vi->postcheckstep) { 2297 ierr = (*vi->postcheckstep)(snes,x,y,w,vi->postcheck,&changed_y,&changed_w);CHKERRQ(ierr); 2298 if (changed_y) { 2299 ierr = VecWAXPY(w,-1.0,y,x);CHKERRQ(ierr); 2300 ierr = SNESVIProjectOntoBounds(snes,w);CHKERRQ(ierr); 2301 } 2302 if (changed_y || changed_w) { /* recompute the function if the step has changed */ 2303 ierr = SNESComputeFunction(snes,w,g); 2304 if (snes->domainerror) { 2305 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2306 PetscFunctionReturn(0); 2307 } 2308 if (snes->ops->solve != SNESSolveVI_SS) { 2309 ierr = SNESVIComputeInactiveSetFnorm(snes,g,w,gnorm);CHKERRQ(ierr); 2310 } else { 2311 ierr = VecNorm(g,NORM_2,gnorm);CHKERRQ(ierr); 2312 } 2313 2314 ierr = VecNormBegin(y,NORM_2,ynorm);CHKERRQ(ierr); 2315 ierr = VecNormEnd(y,NORM_2,ynorm);CHKERRQ(ierr); 2316 if (PetscIsInfOrNanReal(*gnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 2317 } 2318 } 2319 ierr = PetscLogEventEnd(SNES_LineSearch,snes,x,f,g);CHKERRQ(ierr); 2320 PetscFunctionReturn(0); 2321 } 2322 2323 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*/ 2324 /* -------------------------------------------------------------------------- */ 2325 EXTERN_C_BEGIN 2326 #undef __FUNCT__ 2327 #define __FUNCT__ "SNESLineSearchSet_VI" 2328 PetscErrorCode SNESLineSearchSet_VI(SNES snes,FCN2 func,void *lsctx) 2329 { 2330 PetscFunctionBegin; 2331 ((SNES_VI *)(snes->data))->LineSearch = func; 2332 ((SNES_VI *)(snes->data))->lsP = lsctx; 2333 PetscFunctionReturn(0); 2334 } 2335 EXTERN_C_END 2336 2337 /* -------------------------------------------------------------------------- */ 2338 EXTERN_C_BEGIN 2339 #undef __FUNCT__ 2340 #define __FUNCT__ "SNESLineSearchSetMonitor_VI" 2341 PetscErrorCode SNESLineSearchSetMonitor_VI(SNES snes,PetscBool flg) 2342 { 2343 SNES_VI *vi = (SNES_VI*)snes->data; 2344 PetscErrorCode ierr; 2345 2346 PetscFunctionBegin; 2347 if (flg && !vi->lsmonitor) { 2348 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,"stdout",&vi->lsmonitor);CHKERRQ(ierr); 2349 } else if (!flg && vi->lsmonitor) { 2350 ierr = PetscViewerDestroy(&vi->lsmonitor);CHKERRQ(ierr); 2351 } 2352 PetscFunctionReturn(0); 2353 } 2354 EXTERN_C_END 2355 2356 /* 2357 SNESView_VI - Prints info from the SNESVI data structure. 2358 2359 Input Parameters: 2360 . SNES - the SNES context 2361 . viewer - visualization context 2362 2363 Application Interface Routine: SNESView() 2364 */ 2365 #undef __FUNCT__ 2366 #define __FUNCT__ "SNESView_VI" 2367 static PetscErrorCode SNESView_VI(SNES snes,PetscViewer viewer) 2368 { 2369 SNES_VI *vi = (SNES_VI *)snes->data; 2370 const char *cstr,*tstr; 2371 PetscErrorCode ierr; 2372 PetscBool iascii; 2373 2374 PetscFunctionBegin; 2375 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2376 if (iascii) { 2377 if (vi->LineSearch == SNESLineSearchNo_VI) cstr = "SNESLineSearchNo"; 2378 else if (vi->LineSearch == SNESLineSearchQuadratic_VI) cstr = "SNESLineSearchQuadratic"; 2379 else if (vi->LineSearch == SNESLineSearchCubic_VI) cstr = "SNESLineSearchCubic"; 2380 else cstr = "unknown"; 2381 if (snes->ops->solve == SNESSolveVI_SS) tstr = "Semismooth"; 2382 else if (snes->ops->solve == SNESSolveVI_RS) tstr = "Reduced Space"; 2383 else if (snes->ops->solve == SNESSolveVI_RSAUG) tstr = "Reduced space with augmented variables"; 2384 else tstr = "unknown"; 2385 ierr = PetscViewerASCIIPrintf(viewer," VI algorithm: %s\n",tstr);CHKERRQ(ierr); 2386 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 2387 ierr = PetscViewerASCIIPrintf(viewer," alpha=%G, maxstep=%G, minlambda=%G\n",vi->alpha,vi->maxstep,vi->minlambda);CHKERRQ(ierr); 2388 } 2389 PetscFunctionReturn(0); 2390 } 2391 2392 #undef __FUNCT__ 2393 #define __FUNCT__ "SNESVISetVariableBounds" 2394 /*@ 2395 SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. xl <= x <= xu. 2396 2397 Input Parameters: 2398 . snes - the SNES context. 2399 . xl - lower bound. 2400 . xu - upper bound. 2401 2402 Notes: 2403 If this routine is not called then the lower and upper bounds are set to 2404 SNES_VI_INF and SNES_VI_NINF respectively during SNESSetUp(). 2405 2406 @*/ 2407 PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu) 2408 { 2409 SNES_VI *vi; 2410 PetscErrorCode ierr; 2411 2412 PetscFunctionBegin; 2413 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2414 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 2415 PetscValidHeaderSpecific(xu,VEC_CLASSID,3); 2416 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call SNESSetFunction() first"); 2417 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); 2418 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); 2419 ierr = SNESSetType(snes,SNESVI);CHKERRQ(ierr); 2420 vi = (SNES_VI*)snes->data; 2421 ierr = PetscObjectReference((PetscObject)xl);CHKERRQ(ierr); 2422 ierr = PetscObjectReference((PetscObject)xu);CHKERRQ(ierr); 2423 ierr = VecDestroy(&vi->xl);CHKERRQ(ierr); 2424 ierr = VecDestroy(&vi->xu);CHKERRQ(ierr); 2425 vi->xl = xl; 2426 vi->xu = xu; 2427 PetscFunctionReturn(0); 2428 } 2429 2430 /* -------------------------------------------------------------------------- */ 2431 /* 2432 SNESSetFromOptions_VI - Sets various parameters for the SNESVI method. 2433 2434 Input Parameter: 2435 . snes - the SNES context 2436 2437 Application Interface Routine: SNESSetFromOptions() 2438 */ 2439 #undef __FUNCT__ 2440 #define __FUNCT__ "SNESSetFromOptions_VI" 2441 static PetscErrorCode SNESSetFromOptions_VI(SNES snes) 2442 { 2443 SNES_VI *vi = (SNES_VI *)snes->data; 2444 const char *lses[] = {"basic","basicnonorms","quadratic","cubic"}; 2445 const char *vies[] = {"ss","rs","rsaug"}; 2446 PetscErrorCode ierr; 2447 PetscInt indx; 2448 PetscBool flg,set,flg2; 2449 2450 PetscFunctionBegin; 2451 ierr = PetscOptionsHead("SNES semismooth method options");CHKERRQ(ierr); 2452 ierr = PetscOptionsBool("-snes_vi_monitor","Monitor all non-active variables","None",PETSC_FALSE,&flg,0);CHKERRQ(ierr); 2453 if (flg) { 2454 ierr = SNESMonitorSet(snes,SNESMonitorVI,0,0);CHKERRQ(ierr); 2455 } 2456 ierr = PetscOptionsReal("-snes_ls_alpha","Function norm must decrease by","None",vi->alpha,&vi->alpha,0);CHKERRQ(ierr); 2457 ierr = PetscOptionsReal("-snes_ls_maxstep","Step must be less than","None",vi->maxstep,&vi->maxstep,0);CHKERRQ(ierr); 2458 ierr = PetscOptionsReal("-snes_ls_minlambda","Minimum lambda allowed","None",vi->minlambda,&vi->minlambda,0);CHKERRQ(ierr); 2459 ierr = PetscOptionsReal("-snes_vi_const_tol","constraint tolerance","None",vi->const_tol,&vi->const_tol,0);CHKERRQ(ierr); 2460 ierr = PetscOptionsBool("-snes_ls_monitor","Print progress of line searches","SNESLineSearchSetMonitor",vi->lsmonitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); 2461 if (set) {ierr = SNESLineSearchSetMonitor(snes,flg);CHKERRQ(ierr);} 2462 ierr = PetscOptionsEList("-snes_vi_type","Semismooth algorithm used","",vies,3,"ss",&indx,&flg2);CHKERRQ(ierr); 2463 if (flg2) { 2464 switch (indx) { 2465 case 0: 2466 snes->ops->solve = SNESSolveVI_SS; 2467 break; 2468 case 1: 2469 snes->ops->solve = SNESSolveVI_RS; 2470 break; 2471 case 2: 2472 snes->ops->solve = SNESSolveVI_RSAUG; 2473 } 2474 } 2475 ierr = PetscOptionsEList("-snes_ls","Line search used","SNESLineSearchSet",lses,4,"cubic",&indx,&flg);CHKERRQ(ierr); 2476 if (flg) { 2477 switch (indx) { 2478 case 0: 2479 ierr = SNESLineSearchSet(snes,SNESLineSearchNo_VI,PETSC_NULL);CHKERRQ(ierr); 2480 break; 2481 case 1: 2482 ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms_VI,PETSC_NULL);CHKERRQ(ierr); 2483 break; 2484 case 2: 2485 ierr = SNESLineSearchSet(snes,SNESLineSearchQuadratic_VI,PETSC_NULL);CHKERRQ(ierr); 2486 break; 2487 case 3: 2488 ierr = SNESLineSearchSet(snes,SNESLineSearchCubic_VI,PETSC_NULL);CHKERRQ(ierr); 2489 break; 2490 } 2491 } 2492 ierr = PetscOptionsTail();CHKERRQ(ierr); 2493 PetscFunctionReturn(0); 2494 } 2495 /* -------------------------------------------------------------------------- */ 2496 /*MC 2497 SNESVI - Various solvers for variational inequalities based on Newton's method 2498 2499 Options Database: 2500 + -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 2501 additional variables that enforce the constraints 2502 - -snes_vi_monitor - prints the number of active constraints at each iteration. 2503 2504 2505 Level: beginner 2506 2507 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 2508 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 2509 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 2510 2511 M*/ 2512 EXTERN_C_BEGIN 2513 #undef __FUNCT__ 2514 #define __FUNCT__ "SNESCreate_VI" 2515 PetscErrorCode SNESCreate_VI(SNES snes) 2516 { 2517 PetscErrorCode ierr; 2518 SNES_VI *vi; 2519 2520 PetscFunctionBegin; 2521 snes->ops->reset = SNESReset_VI; 2522 snes->ops->setup = SNESSetUp_VI; 2523 snes->ops->solve = SNESSolveVI_RS; 2524 snes->ops->destroy = SNESDestroy_VI; 2525 snes->ops->setfromoptions = SNESSetFromOptions_VI; 2526 snes->ops->view = SNESView_VI; 2527 snes->ops->converged = SNESDefaultConverged_VI; 2528 2529 ierr = PetscNewLog(snes,SNES_VI,&vi);CHKERRQ(ierr); 2530 snes->data = (void*)vi; 2531 vi->alpha = 1.e-4; 2532 vi->maxstep = 1.e8; 2533 vi->minlambda = 1.e-12; 2534 vi->LineSearch = SNESLineSearchCubic_VI; 2535 vi->lsP = PETSC_NULL; 2536 vi->postcheckstep = PETSC_NULL; 2537 vi->postcheck = PETSC_NULL; 2538 vi->precheckstep = PETSC_NULL; 2539 vi->precheck = PETSC_NULL; 2540 vi->const_tol = 2.2204460492503131e-16; 2541 vi->checkredundancy = PETSC_NULL; 2542 2543 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetMonitor_C","SNESLineSearchSetMonitor_VI",SNESLineSearchSetMonitor_VI);CHKERRQ(ierr); 2544 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSet_C","SNESLineSearchSet_VI",SNESLineSearchSet_VI);CHKERRQ(ierr); 2545 2546 PetscFunctionReturn(0); 2547 } 2548 EXTERN_C_END 2549 2550 #undef __FUNCT__ 2551 #define __FUNCT__ "SNESVIGetInactiveSet" 2552 /* 2553 SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear 2554 system is solved on) 2555 2556 Input parameter 2557 . snes - the SNES context 2558 2559 Output parameter 2560 . ISact - active set index set 2561 2562 */ 2563 PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS* inact) 2564 { 2565 SNES_VI *vi = (SNES_VI*)snes->data; 2566 PetscFunctionBegin; 2567 *inact = vi->IS_inact_prev; 2568 PetscFunctionReturn(0); 2569 } 2570