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