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