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