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