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