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