1 2 #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/ 3 #include <petsc-private/kspimpl.h> 4 #include <petsc-private/matimpl.h> 5 #include <petsc-private/dmimpl.h> 6 #include <petsc-private/vecimpl.h> 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "SNESVIGetInactiveSet" 10 /* 11 SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear 12 system is solved on) 13 14 Input parameter 15 . snes - the SNES context 16 17 Output parameter 18 . ISact - active set index set 19 20 */ 21 PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact) 22 { 23 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 24 25 PetscFunctionBegin; 26 *inact = vi->IS_inact_prev; 27 PetscFunctionReturn(0); 28 } 29 30 /* 31 Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors 32 defined by the reduced space method. 33 34 Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints. 35 36 <*/ 37 typedef struct { 38 PetscInt n; /* size of vectors in the reduced DM space */ 39 IS inactive; 40 41 PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */ 42 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*); 43 PetscErrorCode (*createglobalvector)(DM,Vec*); 44 45 DM dm; /* when destroying this object we need to reset the above function into the base DM */ 46 } DM_SNESVI; 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMCreateGlobalVector_SNESVI" 50 /* 51 DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space 52 53 */ 54 PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec) 55 { 56 PetscErrorCode ierr; 57 PetscContainer isnes; 58 DM_SNESVI *dmsnesvi; 59 60 PetscFunctionBegin; 61 ierr = PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);CHKERRQ(ierr); 62 if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing"); 63 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi);CHKERRQ(ierr); 64 ierr = VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "DMCreateInterpolation_SNESVI" 70 /* 71 DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints. 72 73 */ 74 PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec) 75 { 76 PetscErrorCode ierr; 77 PetscContainer isnes; 78 DM_SNESVI *dmsnesvi1,*dmsnesvi2; 79 Mat interp; 80 81 PetscFunctionBegin; 82 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);CHKERRQ(ierr); 83 if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 84 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 85 ierr = PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);CHKERRQ(ierr); 86 if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 87 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);CHKERRQ(ierr); 88 89 ierr = (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);CHKERRQ(ierr); 90 ierr = MatGetSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 91 ierr = MatDestroy(&interp);CHKERRQ(ierr); 92 *vec = 0; 93 PetscFunctionReturn(0); 94 } 95 96 extern PetscErrorCode DMSetVI(DM,IS); 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "DMCoarsen_SNESVI" 100 /* 101 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 102 103 */ 104 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 105 { 106 PetscErrorCode ierr; 107 PetscContainer isnes; 108 DM_SNESVI *dmsnesvi1; 109 Vec finemarked,coarsemarked; 110 IS inactive; 111 Mat inject; 112 VecScatter vscat; 113 const PetscInt *index; 114 PetscInt n,k,cnt = 0,rstart,*coarseindex; 115 PetscScalar *marked; 116 117 PetscFunctionBegin; 118 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);CHKERRQ(ierr); 119 if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing"); 120 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 121 122 /* get the original coarsen */ 123 ierr = (*dmsnesvi1->coarsen)(dm1,comm,dm2);CHKERRQ(ierr); 124 125 /* not sure why this extra reference is needed, but without the dm2 disappears too early */ 126 /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */ 127 /* ierr = PetscObjectReference((PetscObject)*dm2);CHKERRQ(ierr);*/ 128 129 /* need to set back global vectors in order to use the original injection */ 130 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 131 132 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector; 133 134 ierr = DMCreateGlobalVector(dm1,&finemarked);CHKERRQ(ierr); 135 ierr = DMCreateGlobalVector(*dm2,&coarsemarked);CHKERRQ(ierr); 136 137 /* 138 fill finemarked with locations of inactive points 139 */ 140 ierr = ISGetIndices(dmsnesvi1->inactive,&index);CHKERRQ(ierr); 141 ierr = ISGetLocalSize(dmsnesvi1->inactive,&n);CHKERRQ(ierr); 142 ierr = VecSet(finemarked,0.0);CHKERRQ(ierr); 143 for (k=0; k<n; k++) { 144 ierr = VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);CHKERRQ(ierr); 145 } 146 ierr = VecAssemblyBegin(finemarked);CHKERRQ(ierr); 147 ierr = VecAssemblyEnd(finemarked);CHKERRQ(ierr); 148 149 ierr = DMCreateInjection(*dm2,dm1,&inject);CHKERRQ(ierr); 150 ierr = MatScatterGetVecScatter(inject,&vscat);CHKERRQ(ierr); 151 ierr = VecScatterBegin(vscat,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 152 ierr = VecScatterEnd(vscat,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 153 ierr = MatDestroy(&inject);CHKERRQ(ierr); 154 155 /* 156 create index set list of coarse inactive points from coarsemarked 157 */ 158 ierr = VecGetLocalSize(coarsemarked,&n);CHKERRQ(ierr); 159 ierr = VecGetOwnershipRange(coarsemarked,&rstart,NULL);CHKERRQ(ierr); 160 ierr = VecGetArray(coarsemarked,&marked);CHKERRQ(ierr); 161 for (k=0; k<n; k++) { 162 if (marked[k] != 0.0) cnt++; 163 } 164 ierr = PetscMalloc1(cnt,&coarseindex);CHKERRQ(ierr); 165 cnt = 0; 166 for (k=0; k<n; k++) { 167 if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart; 168 } 169 ierr = VecRestoreArray(coarsemarked,&marked);CHKERRQ(ierr); 170 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);CHKERRQ(ierr); 171 172 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 173 174 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 175 176 ierr = DMSetVI(*dm2,inactive);CHKERRQ(ierr); 177 178 ierr = VecDestroy(&finemarked);CHKERRQ(ierr); 179 ierr = VecDestroy(&coarsemarked);CHKERRQ(ierr); 180 ierr = ISDestroy(&inactive);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "DMDestroy_SNESVI" 186 PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi) 187 { 188 PetscErrorCode ierr; 189 190 PetscFunctionBegin; 191 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */ 192 dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation; 193 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen; 194 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector; 195 /* need to clear out this vectors because some of them may not have a reference to the DM 196 but they are counted as having references to the DM in DMDestroy() */ 197 ierr = DMClearGlobalVectors(dmsnesvi->dm);CHKERRQ(ierr); 198 199 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 200 ierr = PetscFree(dmsnesvi);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "DMSetVI" 206 /* 207 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 208 be restricted to only those variables NOT associated with active constraints. 209 210 */ 211 PetscErrorCode DMSetVI(DM dm,IS inactive) 212 { 213 PetscErrorCode ierr; 214 PetscContainer isnes; 215 DM_SNESVI *dmsnesvi; 216 217 PetscFunctionBegin; 218 if (!dm) PetscFunctionReturn(0); 219 220 ierr = PetscObjectReference((PetscObject)inactive);CHKERRQ(ierr); 221 222 ierr = PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);CHKERRQ(ierr); 223 if (!isnes) { 224 ierr = PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);CHKERRQ(ierr); 225 ierr = PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);CHKERRQ(ierr); 226 ierr = PetscNew(&dmsnesvi);CHKERRQ(ierr); 227 ierr = PetscContainerSetPointer(isnes,(void*)dmsnesvi);CHKERRQ(ierr); 228 ierr = PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);CHKERRQ(ierr); 229 ierr = PetscContainerDestroy(&isnes);CHKERRQ(ierr); 230 231 dmsnesvi->createinterpolation = dm->ops->createinterpolation; 232 dm->ops->createinterpolation = DMCreateInterpolation_SNESVI; 233 dmsnesvi->coarsen = dm->ops->coarsen; 234 dm->ops->coarsen = DMCoarsen_SNESVI; 235 dmsnesvi->createglobalvector = dm->ops->createglobalvector; 236 dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 237 } else { 238 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi);CHKERRQ(ierr); 239 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 240 } 241 ierr = DMClearGlobalVectors(dm);CHKERRQ(ierr); 242 ierr = ISGetLocalSize(inactive,&dmsnesvi->n);CHKERRQ(ierr); 243 244 dmsnesvi->inactive = inactive; 245 dmsnesvi->dm = dm; 246 PetscFunctionReturn(0); 247 } 248 249 #undef __FUNCT__ 250 #define __FUNCT__ "DMDestroyVI" 251 /* 252 DMDestroyVI - Frees the DM_SNESVI object contained in the DM 253 - also resets the function pointers in the DM for createinterpolation() etc to use the original DM 254 */ 255 PetscErrorCode DMDestroyVI(DM dm) 256 { 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 if (!dm) PetscFunctionReturn(0); 261 ierr = PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);CHKERRQ(ierr); 262 PetscFunctionReturn(0); 263 } 264 265 /* --------------------------------------------------------------------------------------------------------*/ 266 267 268 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "SNESCreateIndexSets_VINEWTONRSLS" 272 PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact) 273 { 274 PetscErrorCode ierr; 275 276 PetscFunctionBegin; 277 ierr = SNESVIGetActiveSetIS(snes,X,F,ISact);CHKERRQ(ierr); 278 ierr = ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 283 #undef __FUNCT__ 284 #define __FUNCT__ "SNESCreateSubVectors_VINEWTONRSLS" 285 PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv) 286 { 287 PetscErrorCode ierr; 288 Vec v; 289 290 PetscFunctionBegin; 291 ierr = VecCreate(PetscObjectComm((PetscObject)snes),&v);CHKERRQ(ierr); 292 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 293 ierr = VecSetType(v,VECSTANDARD);CHKERRQ(ierr); 294 *newv = v; 295 PetscFunctionReturn(0); 296 } 297 298 /* Resets the snes PC and KSP when the active set sizes change */ 299 #undef __FUNCT__ 300 #define __FUNCT__ "SNESVIResetPCandKSP" 301 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 302 { 303 PetscErrorCode ierr; 304 KSP snesksp; 305 306 PetscFunctionBegin; 307 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 308 ierr = KSPReset(snesksp);CHKERRQ(ierr); 309 310 /* 311 KSP kspnew; 312 PC pcnew; 313 const MatSolverPackage stype; 314 315 316 ierr = KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);CHKERRQ(ierr); 317 kspnew->pc_side = snesksp->pc_side; 318 kspnew->rtol = snesksp->rtol; 319 kspnew->abstol = snesksp->abstol; 320 kspnew->max_it = snesksp->max_it; 321 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 322 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 323 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 324 ierr = PCSetOperators(kspnew->pc,Amat,Pmat);CHKERRQ(ierr); 325 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 326 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 327 ierr = KSPDestroy(&snesksp);CHKERRQ(ierr); 328 snes->ksp = kspnew; 329 ierr = PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);CHKERRQ(ierr); 330 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr);*/ 331 PetscFunctionReturn(0); 332 } 333 334 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 335 implemented in this algorithm. It basically identifies the active constraints and does 336 a linear solve on the other variables (those not associated with the active constraints). */ 337 #undef __FUNCT__ 338 #define __FUNCT__ "SNESSolve_VINEWTONRSLS" 339 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 340 { 341 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 342 PetscErrorCode ierr; 343 PetscInt maxits,i,lits; 344 PetscBool lssucceed; 345 PetscReal fnorm,gnorm,xnorm=0,ynorm; 346 Vec Y,X,F; 347 KSPConvergedReason kspreason; 348 349 PetscFunctionBegin; 350 snes->numFailures = 0; 351 snes->numLinearSolveFailures = 0; 352 snes->reason = SNES_CONVERGED_ITERATING; 353 354 maxits = snes->max_its; /* maximum number of iterations */ 355 X = snes->vec_sol; /* solution vector */ 356 F = snes->vec_func; /* residual vector */ 357 Y = snes->work[0]; /* work vectors */ 358 359 ierr = SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);CHKERRQ(ierr); 360 ierr = SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 361 ierr = SNESLineSearchSetUp(snes->linesearch);CHKERRQ(ierr); 362 363 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 364 snes->iter = 0; 365 snes->norm = 0.0; 366 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 367 368 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 369 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 370 if (snes->domainerror) { 371 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 372 PetscFunctionReturn(0); 373 } 374 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 375 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 376 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 377 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 378 379 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 380 snes->norm = fnorm; 381 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 382 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 383 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 384 385 /* test convergence */ 386 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 387 if (snes->reason) PetscFunctionReturn(0); 388 389 390 for (i=0; i<maxits; i++) { 391 392 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 393 IS IS_redact; /* redundant active set */ 394 VecScatter scat_act,scat_inact; 395 PetscInt nis_act,nis_inact; 396 Vec Y_act,Y_inact,F_inact; 397 Mat jac_inact_inact,prejac_inact_inact; 398 PetscBool isequal; 399 400 /* Call general purpose update function */ 401 if (snes->ops->update) { 402 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 403 } 404 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 405 406 407 /* Create active and inactive index sets */ 408 409 /*original 410 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 411 */ 412 ierr = SNESVIGetActiveSetIS(snes,X,F,&IS_act);CHKERRQ(ierr); 413 414 if (vi->checkredundancy) { 415 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);CHKERRQ(ierr); 416 if (IS_redact) { 417 ierr = ISSort(IS_redact);CHKERRQ(ierr); 418 ierr = ISComplement(IS_redact,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 419 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 420 } else { 421 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 422 } 423 } else { 424 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 425 } 426 427 428 /* Create inactive set submatrix */ 429 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 430 431 if (0) { /* Dead code (temporary developer hack) */ 432 IS keptrows; 433 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 434 if (keptrows) { 435 PetscInt cnt,*nrows,k; 436 const PetscInt *krows,*inact; 437 PetscInt rstart=jac_inact_inact->rmap->rstart; 438 439 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 440 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 441 442 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 443 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 444 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 445 ierr = PetscMalloc1(cnt,&nrows);CHKERRQ(ierr); 446 for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart]; 447 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 448 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 449 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 450 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 451 452 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 453 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 454 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 455 } 456 } 457 ierr = DMSetVI(snes->dm,IS_inact);CHKERRQ(ierr); 458 /* remove later */ 459 460 /* 461 ierr = VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));CHKERRQ(ierr); 462 ierr = VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));CHKERRQ(ierr); 463 ierr = VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));CHKERRQ(ierr); 464 ierr = VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));CHKERRQ(ierr); 465 ierr = ISView(IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)IS_inact)));CHKERRQ(ierr); 466 */ 467 468 /* Get sizes of active and inactive sets */ 469 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 470 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 471 472 /* Create active and inactive set vectors */ 473 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);CHKERRQ(ierr); 474 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);CHKERRQ(ierr); 475 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 476 477 /* Create scatter contexts */ 478 ierr = VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);CHKERRQ(ierr); 479 ierr = VecScatterCreate(Y,IS_inact,Y_inact,NULL,&scat_inact);CHKERRQ(ierr); 480 481 /* Do a vec scatter to active and inactive set vectors */ 482 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 483 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 484 485 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 486 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 487 488 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 489 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 490 491 /* Active set direction = 0 */ 492 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 493 if (snes->jacobian != snes->jacobian_pre) { 494 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 495 } else prejac_inact_inact = jac_inact_inact; 496 497 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 498 if (!isequal) { 499 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 500 } 501 502 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 503 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 504 /* ierr = MatView(snes->jacobian_pre,0); */ 505 506 507 508 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 509 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 510 { 511 PC pc; 512 PetscBool flg; 513 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 514 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 515 if (flg) { 516 KSP *subksps; 517 ierr = PCFieldSplitGetSubKSP(pc,NULL,&subksps);CHKERRQ(ierr); 518 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 519 ierr = PetscFree(subksps);CHKERRQ(ierr); 520 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 521 if (flg) { 522 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 523 const PetscInt *ii; 524 525 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 526 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 527 for (j=0; j<n; j++) { 528 if (ii[j] < N) cnts[0]++; 529 else if (ii[j] < 2*N) cnts[1]++; 530 else if (ii[j] < 3*N) cnts[2]++; 531 } 532 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 533 534 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 535 } 536 } 537 } 538 539 ierr = KSPSolve(snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 540 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 541 if (kspreason < 0) { 542 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 543 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 544 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 545 break; 546 } 547 } 548 549 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 550 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 551 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 552 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 553 554 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 555 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 556 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 557 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 558 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 559 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 560 if (!isequal) { 561 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 562 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 563 } 564 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 565 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 566 if (snes->jacobian != snes->jacobian_pre) { 567 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 568 } 569 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 570 snes->linear_its += lits; 571 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 572 /* 573 if (snes->ops->precheck) { 574 PetscBool changed_y = PETSC_FALSE; 575 ierr = (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr); 576 } 577 578 if (PetscLogPrintInfo) { 579 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 580 } 581 */ 582 /* Compute a (scaled) negative update in the line search routine: 583 Y <- X - lambda*Y 584 and evaluate G = function(Y) (depends on the line search). 585 */ 586 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 587 ynorm = 1; gnorm = fnorm; 588 ierr = SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);CHKERRQ(ierr); 589 ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr); 590 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 591 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);CHKERRQ(ierr); 592 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 593 if (snes->domainerror) { 594 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 595 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 596 PetscFunctionReturn(0); 597 } 598 if (!lssucceed) { 599 if (++snes->numFailures >= snes->maxFailures) { 600 PetscBool ismin; 601 snes->reason = SNES_DIVERGED_LINE_SEARCH; 602 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);CHKERRQ(ierr); 603 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 604 break; 605 } 606 } 607 /* Update function and solution vectors */ 608 fnorm = gnorm; 609 /* Monitor convergence */ 610 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 611 snes->iter = i+1; 612 snes->norm = fnorm; 613 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 614 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 615 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 616 /* Test for convergence, xnorm = || X || */ 617 if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 618 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 619 if (snes->reason) break; 620 } 621 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 622 if (i == maxits) { 623 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 624 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 625 } 626 PetscFunctionReturn(0); 627 } 628 629 #undef __FUNCT__ 630 #define __FUNCT__ "SNESVISetRedundancyCheck" 631 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 632 { 633 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 634 635 PetscFunctionBegin; 636 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 637 vi->checkredundancy = func; 638 vi->ctxP = ctx; 639 PetscFunctionReturn(0); 640 } 641 642 #if defined(PETSC_HAVE_MATLAB_ENGINE) 643 #include <engine.h> 644 #include <mex.h> 645 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 646 647 #undef __FUNCT__ 648 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 649 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx) 650 { 651 PetscErrorCode ierr; 652 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 653 int nlhs = 1, nrhs = 5; 654 mxArray *plhs[1], *prhs[5]; 655 long long int l1 = 0, l2 = 0, ls = 0; 656 PetscInt *indices=NULL; 657 658 PetscFunctionBegin; 659 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 660 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 661 PetscValidPointer(is_redact,3); 662 PetscCheckSameComm(snes,1,is_act,2); 663 664 /* Create IS for reduced active set of size 0, its size and indices will 665 bet set by the Matlab function */ 666 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 667 /* call Matlab function in ctx */ 668 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 669 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 670 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 671 prhs[0] = mxCreateDoubleScalar((double)ls); 672 prhs[1] = mxCreateDoubleScalar((double)l1); 673 prhs[2] = mxCreateDoubleScalar((double)l2); 674 prhs[3] = mxCreateString(sctx->funcname); 675 prhs[4] = sctx->ctx; 676 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 677 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 678 mxDestroyArray(prhs[0]); 679 mxDestroyArray(prhs[1]); 680 mxDestroyArray(prhs[2]); 681 mxDestroyArray(prhs[3]); 682 mxDestroyArray(plhs[0]); 683 PetscFunctionReturn(0); 684 } 685 686 #undef __FUNCT__ 687 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 688 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx) 689 { 690 PetscErrorCode ierr; 691 SNESMatlabContext *sctx; 692 693 PetscFunctionBegin; 694 /* currently sctx is memory bleed */ 695 ierr = PetscNew(&sctx);CHKERRQ(ierr); 696 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 697 sctx->ctx = mxDuplicateArray(ctx); 698 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 699 PetscFunctionReturn(0); 700 } 701 702 #endif 703 704 /* -------------------------------------------------------------------------- */ 705 /* 706 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 707 of the SNESVI nonlinear solver. 708 709 Input Parameter: 710 . snes - the SNES context 711 . x - the solution vector 712 713 Application Interface Routine: SNESSetUp() 714 715 Notes: 716 For basic use of the SNES solvers, the user need not explicitly call 717 SNESSetUp(), since these actions will automatically occur during 718 the call to SNESSolve(). 719 */ 720 #undef __FUNCT__ 721 #define __FUNCT__ "SNESSetUp_VINEWTONRSLS" 722 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 723 { 724 PetscErrorCode ierr; 725 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 726 PetscInt *indices; 727 PetscInt i,n,rstart,rend; 728 SNESLineSearch linesearch; 729 730 PetscFunctionBegin; 731 ierr = SNESSetUp_VI(snes);CHKERRQ(ierr); 732 733 /* Set up previous active index set for the first snes solve 734 vi->IS_inact_prev = 0,1,2,....N */ 735 736 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 737 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 738 ierr = PetscMalloc1(n,&indices);CHKERRQ(ierr); 739 for (i=0; i < n; i++) indices[i] = rstart + i; 740 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);CHKERRQ(ierr); 741 742 /* set the line search functions */ 743 if (!snes->linesearch) { 744 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 745 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 746 } 747 PetscFunctionReturn(0); 748 } 749 /* -------------------------------------------------------------------------- */ 750 #undef __FUNCT__ 751 #define __FUNCT__ "SNESReset_VINEWTONRSLS" 752 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 753 { 754 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 755 PetscErrorCode ierr; 756 757 PetscFunctionBegin; 758 ierr = SNESReset_VI(snes);CHKERRQ(ierr); 759 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 760 PetscFunctionReturn(0); 761 } 762 763 /* -------------------------------------------------------------------------- */ 764 /*MC 765 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 766 767 Options Database: 768 + -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 additional variables that enforce the constraints 769 - -snes_vi_monitor - prints the number of active constraints at each iteration. 770 771 Level: beginner 772 773 References: 774 - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large-Scale 775 Applications, Optimization Methods and Software, 21 (2006). 776 777 .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(), 778 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 779 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 780 781 M*/ 782 #undef __FUNCT__ 783 #define __FUNCT__ "SNESCreate_VINEWTONRSLS" 784 PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 785 { 786 PetscErrorCode ierr; 787 SNES_VINEWTONRSLS *vi; 788 789 PetscFunctionBegin; 790 snes->ops->reset = SNESReset_VINEWTONRSLS; 791 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 792 snes->ops->solve = SNESSolve_VINEWTONRSLS; 793 snes->ops->destroy = SNESDestroy_VI; 794 snes->ops->setfromoptions = SNESSetFromOptions_VI; 795 snes->ops->view = NULL; 796 snes->ops->converged = SNESConvergedDefault_VI; 797 798 snes->usesksp = PETSC_TRUE; 799 snes->usespc = PETSC_FALSE; 800 801 ierr = PetscNewLog(snes,&vi);CHKERRQ(ierr); 802 snes->data = (void*)vi; 803 vi->checkredundancy = NULL; 804 805 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);CHKERRQ(ierr); 806 ierr = PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810