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