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