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