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