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