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