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