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