1 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdm.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "GSDestroy_Private" 6 PetscErrorCode GSDestroy_Private(ISColoring coloring) 7 { 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 12 PetscFunctionReturn(0); 13 14 } 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "SNESComputeGSDefaultSecant" 18 PETSC_EXTERN PetscErrorCode SNESComputeGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx) 19 { 20 21 PetscErrorCode ierr; 22 PetscInt i,j,k,ncolors; 23 DM dm; 24 PetscBool flg; 25 ISColoring coloring; 26 MatColoring mc; 27 Vec W,G,F; 28 PetscScalar h=1e-6; 29 IS *coloris; 30 PetscScalar f,g,x,w,d; 31 const PetscInt *idx; 32 PetscInt size; 33 PetscReal atol,rtol,stol; 34 PetscInt its; 35 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 36 void *fctx; 37 PetscContainer colorcontainer; 38 39 PetscFunctionBegin; 40 if (snes->nwork < 3) { 41 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr); 42 } 43 W = snes->work[0]; 44 G = snes->work[1]; 45 F = snes->work[2]; 46 ierr = SNESGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr); 47 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 48 ierr = SNESGetFunction(snes,NULL,&func,&fctx);CHKERRQ(ierr); 49 ierr = PetscObjectQuery((PetscObject)snes,"SNESGSColoring",(PetscObject*)&colorcontainer);CHKERRQ(ierr); 50 if (!colorcontainer) { 51 /* create the coloring */ 52 ierr = DMHasColoring(dm,&flg);CHKERRQ(ierr); 53 if (flg) { 54 ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);CHKERRQ(ierr); 55 } else { 56 if (!snes->jacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);} 57 ierr = MatColoringCreate(snes->jacobian,&mc);CHKERRQ(ierr); 58 ierr = MatColoringSetDistance(mc,1);CHKERRQ(ierr); 59 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 60 ierr = MatColoringApply(mc,&coloring);CHKERRQ(ierr); 61 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 62 } 63 ierr = PetscContainerCreate(PetscObjectComm((PetscObject)snes),&colorcontainer);CHKERRQ(ierr); 64 ierr = PetscContainerSetPointer(colorcontainer,(void *)coloring);CHKERRQ(ierr); 65 ierr = PetscContainerSetUserDestroy(colorcontainer,(PetscErrorCode (*)(void *))GSDestroy_Private);CHKERRQ(ierr); 66 ierr = PetscObjectCompose((PetscObject)snes,"SNESGSColoring",(PetscObject)colorcontainer);CHKERRQ(ierr); 67 ierr = PetscContainerDestroy(&colorcontainer);CHKERRQ(ierr); 68 } else { 69 ierr = PetscContainerGetPointer(colorcontainer,(void **)&coloring);CHKERRQ(ierr); 70 } 71 ierr = ISColoringGetIS(coloring,&ncolors,&coloris);CHKERRQ(ierr); 72 for (i=0;i<ncolors;i++) { 73 for (k=0;k<its;k++) { 74 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr); 75 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);} 76 ierr = ISGetIndices(coloris[i],&idx);CHKERRQ(ierr); 77 ierr = ISGetLocalSize(coloris[i],&size);CHKERRQ(ierr); 78 ierr = VecCopy(X,W);CHKERRQ(ierr); 79 for (j=0;j<size;j++) { 80 ierr = VecSetValue(W,idx[j],h,ADD_VALUES);CHKERRQ(ierr); 81 } 82 ierr = (*func)(snes,W,G,fctx);CHKERRQ(ierr); 83 if (B) {ierr = VecAXPY(G,-1.0,B);CHKERRQ(ierr);} 84 for (j=0;j<size;j++) { 85 ierr = VecGetValues(F,1,&idx[j],&f);CHKERRQ(ierr); 86 ierr = VecGetValues(X,1,&idx[j],&x);CHKERRQ(ierr); 87 ierr = VecGetValues(G,1,&idx[j],&g);CHKERRQ(ierr); 88 ierr = VecGetValues(W,1,&idx[j],&w);CHKERRQ(ierr); 89 d = (x*g-w*f) / PetscRealPart(g-f); 90 ierr = VecSetValue(X,idx[j],d,INSERT_VALUES);CHKERRQ(ierr); 91 } 92 } 93 } 94 ierr = ISColoringRestoreIS(coloring,&coloris);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97