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