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