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