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