1 2 #include <../src/snes/impls/vi/rs/virsimpl.h> /*I "petscsnes.h" I*/ 3 #include <../include/petsc-private/kspimpl.h> 4 #include <../include/petsc-private/matimpl.h> 5 #include <../include/petsc-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_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 23 24 PetscFunctionBegin; 25 *inact = vi->IS_inact_prev; 26 PetscFunctionReturn(0); 27 } 28 29 /* 30 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 31 defined by the reduced space method. 32 33 Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints. 34 35 <*/ 36 typedef struct { 37 PetscInt n; /* size of vectors in the reduced DM space */ 38 IS inactive; 39 PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */ 40 PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*); 41 PetscErrorCode (*createglobalvector)(DM,Vec*); 42 DM dm; /* when destroying this object we need to reset the above function into the base DM */ 43 } DM_SNESVI; 44 45 #undef __FUNCT__ 46 #define __FUNCT__ "DMCreateGlobalVector_SNESVI" 47 /* 48 DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space 49 50 */ 51 PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec) 52 { 53 PetscErrorCode ierr; 54 PetscContainer isnes; 55 DM_SNESVI *dmsnesvi; 56 57 PetscFunctionBegin; 58 ierr = PetscObjectQuery((PetscObject)dm,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 59 if (!isnes) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_PLIB,"Composed SNES is missing"); 60 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi);CHKERRQ(ierr); 61 ierr = VecCreateMPI(((PetscObject)dm)->comm,dmsnesvi->n,PETSC_DETERMINE,vec);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "DMCreateInterpolation_SNESVI" 67 /* 68 DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints. 69 70 */ 71 PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec) 72 { 73 PetscErrorCode ierr; 74 PetscContainer isnes; 75 DM_SNESVI *dmsnesvi1,*dmsnesvi2; 76 Mat interp; 77 78 PetscFunctionBegin; 79 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 80 if (!isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed VI data structure is missing"); 81 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 82 ierr = PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 83 if (!isnes) SETERRQ(((PetscObject)dm2)->comm,PETSC_ERR_PLIB,"Composed VI data structure is missing"); 84 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);CHKERRQ(ierr); 85 86 ierr = (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,PETSC_NULL);CHKERRQ(ierr); 87 ierr = MatGetSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 88 ierr = MatDestroy(&interp);CHKERRQ(ierr); 89 *vec = 0; 90 PetscFunctionReturn(0); 91 } 92 93 extern PetscErrorCode DMSetVI(DM,IS); 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "DMCoarsen_SNESVI" 97 /* 98 DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set 99 100 */ 101 PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2) 102 { 103 PetscErrorCode ierr; 104 PetscContainer isnes; 105 DM_SNESVI *dmsnesvi1; 106 Vec finemarked,coarsemarked; 107 IS inactive; 108 VecScatter inject; 109 const PetscInt *index; 110 PetscInt n,k,cnt = 0,rstart,*coarseindex; 111 PetscScalar *marked; 112 113 PetscFunctionBegin; 114 ierr = PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 115 if (!isnes) SETERRQ(((PetscObject)dm1)->comm,PETSC_ERR_PLIB,"Composed VI data structure is missing"); 116 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);CHKERRQ(ierr); 117 118 /* get the original coarsen */ 119 ierr = (*dmsnesvi1->coarsen)(dm1,comm,dm2);CHKERRQ(ierr); 120 121 /* not sure why this extra reference is needed, but without the dm2 disappears too early */ 122 ierr = PetscObjectReference((PetscObject)*dm2);CHKERRQ(ierr); 123 124 /* need to set back global vectors in order to use the original injection */ 125 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 126 dm1->ops->createglobalvector = dmsnesvi1->createglobalvector; 127 ierr = DMCreateGlobalVector(dm1,&finemarked);CHKERRQ(ierr); 128 ierr = DMCreateGlobalVector(*dm2,&coarsemarked);CHKERRQ(ierr); 129 130 /* 131 fill finemarked with locations of inactive points 132 */ 133 ierr = ISGetIndices(dmsnesvi1->inactive,&index);CHKERRQ(ierr); 134 ierr = ISGetLocalSize(dmsnesvi1->inactive,&n);CHKERRQ(ierr); 135 ierr = VecSet(finemarked,0.0);CHKERRQ(ierr); 136 for (k=0;k<n;k++) { 137 ierr = VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);CHKERRQ(ierr); 138 } 139 ierr = VecAssemblyBegin(finemarked);CHKERRQ(ierr); 140 ierr = VecAssemblyEnd(finemarked);CHKERRQ(ierr); 141 142 ierr = DMCreateInjection(*dm2,dm1,&inject);CHKERRQ(ierr); 143 ierr = VecScatterBegin(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 144 ierr = VecScatterEnd(inject,finemarked,coarsemarked,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 145 ierr = VecScatterDestroy(&inject);CHKERRQ(ierr); 146 147 /* 148 create index set list of coarse inactive points from coarsemarked 149 */ 150 ierr = VecGetLocalSize(coarsemarked,&n);CHKERRQ(ierr); 151 ierr = VecGetOwnershipRange(coarsemarked,&rstart,PETSC_NULL);CHKERRQ(ierr); 152 ierr = VecGetArray(coarsemarked,&marked);CHKERRQ(ierr); 153 for (k=0; k<n; k++) { 154 if (marked[k] != 0.0) cnt++; 155 } 156 ierr = PetscMalloc(cnt*sizeof(PetscInt),&coarseindex);CHKERRQ(ierr); 157 cnt = 0; 158 for (k=0; k<n; k++) { 159 if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart; 160 } 161 ierr = VecRestoreArray(coarsemarked,&marked);CHKERRQ(ierr); 162 ierr = ISCreateGeneral(((PetscObject)coarsemarked)->comm,cnt,coarseindex,PETSC_OWN_POINTER,&inactive);CHKERRQ(ierr); 163 164 ierr = DMClearGlobalVectors(dm1);CHKERRQ(ierr); 165 dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 166 ierr = DMSetVI(*dm2,inactive);CHKERRQ(ierr); 167 168 ierr = VecDestroy(&finemarked);CHKERRQ(ierr); 169 ierr = VecDestroy(&coarsemarked);CHKERRQ(ierr); 170 ierr = ISDestroy(&inactive);CHKERRQ(ierr); 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "DMDestroy_SNESVI" 176 PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi) 177 { 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */ 182 dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation; 183 dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen; 184 dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector; 185 /* need to clear out this vectors because some of them may not have a reference to the DM 186 but they are counted as having references to the DM in DMDestroy() */ 187 ierr = DMClearGlobalVectors(dmsnesvi->dm);CHKERRQ(ierr); 188 189 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 190 ierr = PetscFree(dmsnesvi);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 #undef __FUNCT__ 195 #define __FUNCT__ "DMSetVI" 196 /* 197 DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to 198 be restricted to only those variables NOT associated with active constraints. 199 200 */ 201 PetscErrorCode DMSetVI(DM dm,IS inactive) 202 { 203 PetscErrorCode ierr; 204 PetscContainer isnes; 205 DM_SNESVI *dmsnesvi; 206 207 PetscFunctionBegin; 208 if (!dm) PetscFunctionReturn(0); 209 210 ierr = PetscObjectReference((PetscObject)inactive);CHKERRQ(ierr); 211 212 ierr = PetscObjectQuery((PetscObject)dm,"VI",(PetscObject *)&isnes);CHKERRQ(ierr); 213 if (!isnes) { 214 ierr = PetscContainerCreate(((PetscObject)dm)->comm,&isnes);CHKERRQ(ierr); 215 ierr = PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);CHKERRQ(ierr); 216 ierr = PetscNew(DM_SNESVI,&dmsnesvi);CHKERRQ(ierr); 217 ierr = PetscContainerSetPointer(isnes,(void*)dmsnesvi);CHKERRQ(ierr); 218 ierr = PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);CHKERRQ(ierr); 219 ierr = PetscContainerDestroy(&isnes);CHKERRQ(ierr); 220 dmsnesvi->createinterpolation = dm->ops->createinterpolation; 221 dm->ops->createinterpolation = DMCreateInterpolation_SNESVI; 222 dmsnesvi->coarsen = dm->ops->coarsen; 223 dm->ops->coarsen = DMCoarsen_SNESVI; 224 dmsnesvi->createglobalvector = dm->ops->createglobalvector; 225 dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI; 226 } else { 227 ierr = PetscContainerGetPointer(isnes,(void**)&dmsnesvi);CHKERRQ(ierr); 228 ierr = ISDestroy(&dmsnesvi->inactive);CHKERRQ(ierr); 229 } 230 ierr = DMClearGlobalVectors(dm);CHKERRQ(ierr); 231 ierr = ISGetLocalSize(inactive,&dmsnesvi->n);CHKERRQ(ierr); 232 dmsnesvi->inactive = inactive; 233 dmsnesvi->dm = dm; 234 PetscFunctionReturn(0); 235 } 236 237 #undef __FUNCT__ 238 #define __FUNCT__ "DMDestroyVI" 239 /* 240 DMDestroyVI - Frees the DM_SNESVI object contained in the DM 241 - also resets the function pointers in the DM for createinterpolation() etc to use the original DM 242 */ 243 PetscErrorCode DMDestroyVI(DM dm) 244 { 245 PetscErrorCode ierr; 246 247 PetscFunctionBegin; 248 if (!dm) PetscFunctionReturn(0); 249 ierr = PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)PETSC_NULL);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /* --------------------------------------------------------------------------------------------------------*/ 254 255 256 257 258 #undef __FUNCT__ 259 #define __FUNCT__ "SNESCreateIndexSets_VINEWTONRSLS" 260 PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS* ISact,IS* ISinact) 261 { 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = SNESVIGetActiveSetIS(snes,X,F,ISact);CHKERRQ(ierr); 266 ierr = ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */ 271 #undef __FUNCT__ 272 #define __FUNCT__ "SNESCreateSubVectors_VINEWTONRSLS" 273 PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec* newv) 274 { 275 PetscErrorCode ierr; 276 Vec v; 277 278 PetscFunctionBegin; 279 ierr = VecCreate(((PetscObject)snes)->comm,&v);CHKERRQ(ierr); 280 ierr = VecSetSizes(v,n,PETSC_DECIDE);CHKERRQ(ierr); 281 ierr = VecSetFromOptions(v);CHKERRQ(ierr); 282 *newv = v; 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 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 323 implemented in this algorithm. It basically identifies the active constraints and does 324 a linear solve on the other variables (those not associated with the active constraints). */ 325 #undef __FUNCT__ 326 #define __FUNCT__ "SNESSolve_VINEWTONRSLS" 327 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 328 { 329 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 330 PetscErrorCode ierr; 331 PetscInt maxits,i,lits; 332 PetscBool lssucceed; 333 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 334 PetscReal fnorm,gnorm,xnorm=0,ynorm; 335 Vec Y,X,F; 336 KSPConvergedReason kspreason; 337 338 PetscFunctionBegin; 339 snes->numFailures = 0; 340 snes->numLinearSolveFailures = 0; 341 snes->reason = SNES_CONVERGED_ITERATING; 342 343 maxits = snes->max_its; /* maximum number of iterations */ 344 X = snes->vec_sol; /* solution vector */ 345 F = snes->vec_func; /* residual vector */ 346 Y = snes->work[0]; /* work vectors */ 347 348 ierr = SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);CHKERRQ(ierr); 349 ierr = SNESLineSearchSetVecs(snes->linesearch, X, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 350 ierr = SNESLineSearchSetUp(snes->linesearch);CHKERRQ(ierr); 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(((PetscObject)X)->comm,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(((PetscObject)snes)->comm,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_(((PetscObject)(vi->xu))->comm));CHKERRQ(ierr); 456 ierr = VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));CHKERRQ(ierr); 457 ierr = VecView(X,PETSC_VIEWER_BINARY_(((PetscObject)X)->comm));CHKERRQ(ierr); 458 ierr = VecView(F,PETSC_VIEWER_BINARY_(((PetscObject)F)->comm));CHKERRQ(ierr); 459 ierr = ISView(IS_inact,PETSC_VIEWER_BINARY_(((PetscObject)IS_inact)->comm));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_VINEWTONRSLS(snes,nis_inact,&F_inact);CHKERRQ(ierr); 468 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);CHKERRQ(ierr); 469 ierr = SNESCreateSubVectors_VINEWTONRSLS(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 = PetscObjectTypeCompare((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 = PetscObjectTypeCompare((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->precheck) { 569 PetscBool changed_y = PETSC_FALSE; 570 ierr = (*snes->ops->precheck)(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 = SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);CHKERRQ(ierr); 584 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 585 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); 586 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 587 if (snes->domainerror) { 588 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 589 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr); 593 594 if (!lssucceed) { 595 if (++snes->numFailures >= snes->maxFailures) { 596 PetscBool ismin; 597 snes->reason = SNES_DIVERGED_LINE_SEARCH; 598 ierr = SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);CHKERRQ(ierr); 599 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 600 break; 601 } 602 } 603 /* Update function and solution vectors */ 604 fnorm = gnorm; 605 /* Monitor convergence */ 606 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 607 snes->iter = i+1; 608 snes->norm = fnorm; 609 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 610 SNESLogConvHistory(snes,snes->norm,lits); 611 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 612 /* Test for convergence, xnorm = || X || */ 613 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 614 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 615 if (snes->reason) break; 616 } 617 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 618 if (i == maxits) { 619 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 620 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 621 } 622 PetscFunctionReturn(0); 623 } 624 625 #undef __FUNCT__ 626 #define __FUNCT__ "SNESVISetRedundancyCheck" 627 PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx) 628 { 629 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 630 631 PetscFunctionBegin; 632 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 633 vi->checkredundancy = func; 634 vi->ctxP = ctx; 635 PetscFunctionReturn(0); 636 } 637 638 #if defined(PETSC_HAVE_MATLAB_ENGINE) 639 #include <engine.h> 640 #include <mex.h> 641 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "SNESVIRedundancyCheck_Matlab" 645 PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS* is_redact,void* ctx) 646 { 647 PetscErrorCode ierr; 648 SNESMatlabContext *sctx = (SNESMatlabContext*)ctx; 649 int nlhs = 1, nrhs = 5; 650 mxArray *plhs[1], *prhs[5]; 651 long long int l1 = 0, l2 = 0, ls = 0; 652 PetscInt *indices=PETSC_NULL; 653 654 PetscFunctionBegin; 655 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 656 PetscValidHeaderSpecific(is_act,IS_CLASSID,2); 657 PetscValidPointer(is_redact,3); 658 PetscCheckSameComm(snes,1,is_act,2); 659 660 /* Create IS for reduced active set of size 0, its size and indices will 661 bet set by the Matlab function */ 662 ierr = ISCreateGeneral(((PetscObject)snes)->comm,0,indices,PETSC_OWN_POINTER,is_redact);CHKERRQ(ierr); 663 /* call Matlab function in ctx */ 664 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 665 ierr = PetscMemcpy(&l1,&is_act,sizeof(is_act));CHKERRQ(ierr); 666 ierr = PetscMemcpy(&l2,is_redact,sizeof(is_act));CHKERRQ(ierr); 667 prhs[0] = mxCreateDoubleScalar((double)ls); 668 prhs[1] = mxCreateDoubleScalar((double)l1); 669 prhs[2] = mxCreateDoubleScalar((double)l2); 670 prhs[3] = mxCreateString(sctx->funcname); 671 prhs[4] = sctx->ctx; 672 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");CHKERRQ(ierr); 673 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 674 mxDestroyArray(prhs[0]); 675 mxDestroyArray(prhs[1]); 676 mxDestroyArray(prhs[2]); 677 mxDestroyArray(prhs[3]); 678 mxDestroyArray(plhs[0]); 679 PetscFunctionReturn(0); 680 } 681 682 #undef __FUNCT__ 683 #define __FUNCT__ "SNESVISetRedundancyCheckMatlab" 684 PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char* func,mxArray* ctx) 685 { 686 PetscErrorCode ierr; 687 SNESMatlabContext *sctx; 688 689 PetscFunctionBegin; 690 /* currently sctx is memory bleed */ 691 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 692 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 693 sctx->ctx = mxDuplicateArray(ctx); 694 ierr = SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 #endif 699 700 /* -------------------------------------------------------------------------- */ 701 /* 702 SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use 703 of the SNESVI nonlinear solver. 704 705 Input Parameter: 706 . snes - the SNES context 707 . x - the solution vector 708 709 Application Interface Routine: SNESSetUp() 710 711 Notes: 712 For basic use of the SNES solvers, the user need not explicitly call 713 SNESSetUp(), since these actions will automatically occur during 714 the call to SNESSolve(). 715 */ 716 #undef __FUNCT__ 717 #define __FUNCT__ "SNESSetUp_VINEWTONRSLS" 718 PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes) 719 { 720 PetscErrorCode ierr; 721 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 722 PetscInt *indices; 723 PetscInt i,n,rstart,rend; 724 SNESLineSearch linesearch; 725 726 PetscFunctionBegin; 727 ierr = SNESSetUp_VI(snes);CHKERRQ(ierr); 728 729 /* Set up previous active index set for the first snes solve 730 vi->IS_inact_prev = 0,1,2,....N */ 731 732 ierr = VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);CHKERRQ(ierr); 733 ierr = VecGetLocalSize(snes->vec_sol,&n);CHKERRQ(ierr); 734 ierr = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 735 for (i=0;i < n; i++) indices[i] = rstart + i; 736 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);CHKERRQ(ierr); 737 738 /* set the line search functions */ 739 if (!snes->linesearch) { 740 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 741 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 742 } 743 PetscFunctionReturn(0); 744 } 745 /* -------------------------------------------------------------------------- */ 746 #undef __FUNCT__ 747 #define __FUNCT__ "SNESReset_VINEWTONRSLS" 748 PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes) 749 { 750 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data; 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = SNESReset_VI(snes);CHKERRQ(ierr); 755 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 756 PetscFunctionReturn(0); 757 } 758 759 /* -------------------------------------------------------------------------- */ 760 /*MC 761 SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method 762 763 Options Database: 764 + -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 765 - -snes_vi_monitor - prints the number of active constraints at each iteration. 766 767 Level: beginner 768 769 References: 770 - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large-Scale 771 Applications, Optimization Methods and Software, 21 (2006). 772 773 .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONRSLS, SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSet(), 774 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 775 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 776 777 M*/ 778 EXTERN_C_BEGIN 779 #undef __FUNCT__ 780 #define __FUNCT__ "SNESCreate_VINEWTONRSLS" 781 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 = PETSC_NULL; 793 snes->ops->converged = SNESDefaultConverged_VI; 794 795 snes->usesksp = PETSC_TRUE; 796 snes->usespc = PETSC_FALSE; 797 798 ierr = PetscNewLog(snes,SNES_VINEWTONRSLS,&vi);CHKERRQ(ierr); 799 snes->data = (void*)vi; 800 vi->checkredundancy = PETSC_NULL; 801 802 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESVISetVariableBounds_C","SNESVISetVariableBounds_VI",SNESVISetVariableBounds_VI);CHKERRQ(ierr); 803 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESVISetComputeVariableBounds_C","SNESVISetComputeVariableBounds_VI",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr); 804 PetscFunctionReturn(0); 805 } 806 EXTERN_C_END 807 808