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