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