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