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 284 PetscFunctionReturn(0); 285 } 286 287 /* Resets the snes PC and KSP when the active set sizes change */ 288 #undef __FUNCT__ 289 #define __FUNCT__ "SNESVIResetPCandKSP" 290 PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat) 291 { 292 PetscErrorCode ierr; 293 KSP snesksp; 294 295 PetscFunctionBegin; 296 ierr = SNESGetKSP(snes,&snesksp);CHKERRQ(ierr); 297 ierr = KSPReset(snesksp);CHKERRQ(ierr); 298 299 /* 300 KSP kspnew; 301 PC pcnew; 302 const MatSolverPackage stype; 303 304 305 ierr = KSPCreate(((PetscObject)snes)->comm,&kspnew);CHKERRQ(ierr); 306 kspnew->pc_side = snesksp->pc_side; 307 kspnew->rtol = snesksp->rtol; 308 kspnew->abstol = snesksp->abstol; 309 kspnew->max_it = snesksp->max_it; 310 ierr = KSPSetType(kspnew,((PetscObject)snesksp)->type_name);CHKERRQ(ierr); 311 ierr = KSPGetPC(kspnew,&pcnew);CHKERRQ(ierr); 312 ierr = PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);CHKERRQ(ierr); 313 ierr = PCSetOperators(kspnew->pc,Amat,Pmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 314 ierr = PCFactorGetMatSolverPackage(snesksp->pc,&stype);CHKERRQ(ierr); 315 ierr = PCFactorSetMatSolverPackage(kspnew->pc,stype);CHKERRQ(ierr); 316 ierr = KSPDestroy(&snesksp);CHKERRQ(ierr); 317 snes->ksp = kspnew; 318 ierr = PetscLogObjectParent(snes,kspnew);CHKERRQ(ierr); 319 ierr = KSPSetFromOptions(kspnew);CHKERRQ(ierr);*/ 320 PetscFunctionReturn(0); 321 } 322 323 /* Variational Inequality solver using reduce space method. No semismooth algorithm is 324 implemented in this algorithm. It basically identifies the active constraints and does 325 a linear solve on the other variables (those not associated with the active constraints). */ 326 #undef __FUNCT__ 327 #define __FUNCT__ "SNESSolve_VINEWTONRSLS" 328 PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes) 329 { 330 SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data; 331 PetscErrorCode ierr; 332 PetscInt maxits,i,lits; 333 PetscBool lssucceed; 334 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 335 PetscReal fnorm,gnorm,xnorm=0,ynorm; 336 Vec Y,X,F; 337 KSPConvergedReason kspreason; 338 339 PetscFunctionBegin; 340 snes->numFailures = 0; 341 snes->numLinearSolveFailures = 0; 342 snes->reason = SNES_CONVERGED_ITERATING; 343 344 maxits = snes->max_its; /* maximum number of iterations */ 345 X = snes->vec_sol; /* solution vector */ 346 F = snes->vec_func; /* residual vector */ 347 Y = snes->work[0]; /* work vectors */ 348 349 ierr = SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);CHKERRQ(ierr); 350 ierr = SNESLineSearchSetVecs(snes->linesearch, X, PETSC_NULL, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 351 ierr = SNESLineSearchSetUp(snes->linesearch);CHKERRQ(ierr); 352 353 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 354 snes->iter = 0; 355 snes->norm = 0.0; 356 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 357 358 ierr = SNESVIProjectOntoBounds(snes,X);CHKERRQ(ierr); 359 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 360 if (snes->domainerror) { 361 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 362 PetscFunctionReturn(0); 363 } 364 ierr = SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);CHKERRQ(ierr); 365 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 366 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 367 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(((PetscObject)X)->comm,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 368 369 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 370 snes->norm = fnorm; 371 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 372 SNESLogConvHistory(snes,fnorm,0); 373 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 374 375 /* set parameter for default relative tolerance convergence test */ 376 snes->ttol = fnorm*snes->rtol; 377 /* test convergence */ 378 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 379 if (snes->reason) PetscFunctionReturn(0); 380 381 382 for (i=0; i<maxits; i++) { 383 384 IS IS_act,IS_inact; /* _act -> active set _inact -> inactive set */ 385 IS IS_redact; /* redundant active set */ 386 VecScatter scat_act,scat_inact; 387 PetscInt nis_act,nis_inact; 388 Vec Y_act,Y_inact,F_inact; 389 Mat jac_inact_inact,prejac_inact_inact; 390 PetscBool isequal; 391 392 /* Call general purpose update function */ 393 if (snes->ops->update) { 394 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 395 } 396 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 397 398 399 /* Create active and inactive index sets */ 400 401 /*original 402 ierr = SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&IS_inact);CHKERRQ(ierr); 403 */ 404 ierr = SNESVIGetActiveSetIS(snes,X,F,&IS_act);CHKERRQ(ierr); 405 406 if (vi->checkredundancy) { 407 (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);CHKERRQ(ierr); 408 if (IS_redact) { 409 ierr = ISSort(IS_redact);CHKERRQ(ierr); 410 ierr = ISComplement(IS_redact,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 411 ierr = ISDestroy(&IS_redact);CHKERRQ(ierr); 412 } 413 else { 414 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 415 } 416 } else { 417 ierr = ISComplement(IS_act,X->map->rstart,X->map->rend,&IS_inact);CHKERRQ(ierr); 418 } 419 420 421 /* Create inactive set submatrix */ 422 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 423 424 if (0) { /* Dead code (temporary developer hack) */ 425 IS keptrows; 426 ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); 427 if (keptrows) { 428 PetscInt cnt,*nrows,k; 429 const PetscInt *krows,*inact; 430 PetscInt rstart=jac_inact_inact->rmap->rstart; 431 432 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 433 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 434 435 ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); 436 ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); 437 ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); 438 ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); 439 for (k=0; k<cnt; k++) { 440 nrows[k] = inact[krows[k]-rstart]; 441 } 442 ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); 443 ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); 444 ierr = ISDestroy(&keptrows);CHKERRQ(ierr); 445 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 446 447 ierr = ISCreateGeneral(((PetscObject)snes)->comm,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); 448 ierr = ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); 449 ierr = MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); 450 } 451 } 452 ierr = DMSetVI(snes->dm,IS_inact);CHKERRQ(ierr); 453 /* remove later */ 454 455 /* 456 ierr = VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));CHKERRQ(ierr); 457 ierr = VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));CHKERRQ(ierr); 458 ierr = VecView(X,PETSC_VIEWER_BINARY_(((PetscObject)X)->comm));CHKERRQ(ierr); 459 ierr = VecView(F,PETSC_VIEWER_BINARY_(((PetscObject)F)->comm));CHKERRQ(ierr); 460 ierr = ISView(IS_inact,PETSC_VIEWER_BINARY_(((PetscObject)IS_inact)->comm));CHKERRQ(ierr); 461 */ 462 463 /* Get sizes of active and inactive sets */ 464 ierr = ISGetLocalSize(IS_act,&nis_act);CHKERRQ(ierr); 465 ierr = ISGetLocalSize(IS_inact,&nis_inact);CHKERRQ(ierr); 466 467 /* Create active and inactive set vectors */ 468 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);CHKERRQ(ierr); 469 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);CHKERRQ(ierr); 470 ierr = SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);CHKERRQ(ierr); 471 472 /* Create scatter contexts */ 473 ierr = VecScatterCreate(Y,IS_act,Y_act,PETSC_NULL,&scat_act);CHKERRQ(ierr); 474 ierr = VecScatterCreate(Y,IS_inact,Y_inact,PETSC_NULL,&scat_inact);CHKERRQ(ierr); 475 476 /* Do a vec scatter to active and inactive set vectors */ 477 ierr = VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 478 ierr = VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 479 480 ierr = VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 481 ierr = VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 482 483 ierr = VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 484 ierr = VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 485 486 /* Active set direction = 0 */ 487 ierr = VecSet(Y_act,0);CHKERRQ(ierr); 488 if (snes->jacobian != snes->jacobian_pre) { 489 ierr = MatGetSubMatrix(snes->jacobian_pre,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);CHKERRQ(ierr); 490 } else prejac_inact_inact = jac_inact_inact; 491 492 ierr = ISEqual(vi->IS_inact_prev,IS_inact,&isequal);CHKERRQ(ierr); 493 if (!isequal) { 494 ierr = SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);CHKERRQ(ierr); 495 flg = DIFFERENT_NONZERO_PATTERN; 496 } 497 498 /* ierr = ISView(IS_inact,0);CHKERRQ(ierr); */ 499 /* ierr = ISView(IS_act,0);CHKERRQ(ierr);*/ 500 /* ierr = MatView(snes->jacobian_pre,0); */ 501 502 503 504 ierr = KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact,flg);CHKERRQ(ierr); 505 ierr = KSPSetUp(snes->ksp);CHKERRQ(ierr); 506 { 507 PC pc; 508 PetscBool flg; 509 ierr = KSPGetPC(snes->ksp,&pc);CHKERRQ(ierr); 510 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);CHKERRQ(ierr); 511 if (flg) { 512 KSP *subksps; 513 ierr = PCFieldSplitGetSubKSP(pc,PETSC_NULL,&subksps);CHKERRQ(ierr); 514 ierr = KSPGetPC(subksps[0],&pc);CHKERRQ(ierr); 515 ierr = PetscFree(subksps);CHKERRQ(ierr); 516 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);CHKERRQ(ierr); 517 if (flg) { 518 PetscInt n,N = 101*101,j,cnts[3] = {0,0,0}; 519 const PetscInt *ii; 520 521 ierr = ISGetSize(IS_inact,&n);CHKERRQ(ierr); 522 ierr = ISGetIndices(IS_inact,&ii);CHKERRQ(ierr); 523 for (j=0; j<n; j++) { 524 if (ii[j] < N) cnts[0]++; 525 else if (ii[j] < 2*N) cnts[1]++; 526 else if (ii[j] < 3*N) cnts[2]++; 527 } 528 ierr = ISRestoreIndices(IS_inact,&ii);CHKERRQ(ierr); 529 530 ierr = PCBJacobiSetTotalBlocks(pc,3,cnts);CHKERRQ(ierr); 531 } 532 } 533 } 534 535 ierr = SNES_KSPSolve(snes,snes->ksp,F_inact,Y_inact);CHKERRQ(ierr); 536 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 537 if (kspreason < 0) { 538 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 539 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 540 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 541 break; 542 } 543 } 544 545 ierr = VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 546 ierr = VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 547 ierr = VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 548 ierr = VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 549 550 ierr = VecDestroy(&F_inact);CHKERRQ(ierr); 551 ierr = VecDestroy(&Y_act);CHKERRQ(ierr); 552 ierr = VecDestroy(&Y_inact);CHKERRQ(ierr); 553 ierr = VecScatterDestroy(&scat_act);CHKERRQ(ierr); 554 ierr = VecScatterDestroy(&scat_inact);CHKERRQ(ierr); 555 ierr = ISDestroy(&IS_act);CHKERRQ(ierr); 556 if (!isequal) { 557 ierr = ISDestroy(&vi->IS_inact_prev);CHKERRQ(ierr); 558 ierr = ISDuplicate(IS_inact,&vi->IS_inact_prev);CHKERRQ(ierr); 559 } 560 ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); 561 ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); 562 if (snes->jacobian != snes->jacobian_pre) { 563 ierr = MatDestroy(&prejac_inact_inact);CHKERRQ(ierr); 564 } 565 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 566 snes->linear_its += lits; 567 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 568 /* 569 if (snes->ops->precheck) { 570 PetscBool changed_y = PETSC_FALSE; 571 ierr = (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr); 572 } 573 574 if (PetscLogPrintInfo) { 575 ierr = SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 576 } 577 */ 578 /* Compute a (scaled) negative update in the line search routine: 579 Y <- X - lambda*Y 580 and evaluate G = function(Y) (depends on the line search). 581 */ 582 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 583 ynorm = 1; gnorm = fnorm; 584 ierr = SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);CHKERRQ(ierr); 585 ierr = SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);CHKERRQ(ierr); 586 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); 587 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 588 if (snes->domainerror) { 589 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 590 ierr = DMDestroyVI(snes->dm);CHKERRQ(ierr); 591 PetscFunctionReturn(0); 592 } 593 ierr = SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);CHKERRQ(ierr); 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 = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 608 snes->iter = i+1; 609 snes->norm = fnorm; 610 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 611 SNESLogConvHistory(snes,snes->norm,lits); 612 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 613 /* Test for convergence, xnorm = || X || */ 614 if (snes->ops->converged != SNESSkipConverged) { 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=PETSC_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(((PetscObject)snes)->comm,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 = PetscMalloc(sizeof(SNESMatlabContext),&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 = PetscMalloc(n*sizeof(PetscInt),&indices);CHKERRQ(ierr); 736 for (i=0;i < n; i++) indices[i] = rstart + i; 737 ierr = ISCreateGeneral(((PetscObject)snes)->comm,n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);CHKERRQ(ierr); 738 739 /* set the line search functions */ 740 if (!snes->linesearch) { 741 ierr = SNESGetSNESLineSearch(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 EXTERN_C_BEGIN 780 #undef __FUNCT__ 781 #define __FUNCT__ "SNESCreate_VINEWTONRSLS" 782 PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes) 783 { 784 PetscErrorCode ierr; 785 SNES_VINEWTONRSLS *vi; 786 787 PetscFunctionBegin; 788 snes->ops->reset = SNESReset_VINEWTONRSLS; 789 snes->ops->setup = SNESSetUp_VINEWTONRSLS; 790 snes->ops->solve = SNESSolve_VINEWTONRSLS; 791 snes->ops->destroy = SNESDestroy_VI; 792 snes->ops->setfromoptions = SNESSetFromOptions_VI; 793 snes->ops->view = PETSC_NULL; 794 snes->ops->converged = SNESDefaultConverged_VI; 795 796 snes->usesksp = PETSC_TRUE; 797 snes->usespc = PETSC_FALSE; 798 799 ierr = PetscNewLog(snes,SNES_VINEWTONRSLS,&vi);CHKERRQ(ierr); 800 snes->data = (void*)vi; 801 vi->checkredundancy = PETSC_NULL; 802 803 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESVISetVariableBounds_C","SNESVISetVariableBounds_VI",SNESVISetVariableBounds_VI);CHKERRQ(ierr); 804 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESVISetComputeVariableBounds_C","SNESVISetComputeVariableBounds_VI",SNESVISetComputeVariableBounds_VI);CHKERRQ(ierr); 805 PetscFunctionReturn(0); 806 } 807 EXTERN_C_END 808 809