1 #include <../src/snes/impls/gs/gsimpl.h> 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "GSDestroy_Private" 5 PetscErrorCode GSDestroy_Private(ISColoring coloring) 6 { 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 11 PetscFunctionReturn(0); 12 13 } 14 15 #undef __FUNCT__ 16 #define __FUNCT__ "SNESComputeGSDefaultSecant" 17 PETSC_EXTERN PetscErrorCode SNESComputeGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx) 18 { 19 PetscErrorCode ierr; 20 SNES_GS *gs = (SNES_GS*)snes->data; 21 PetscInt i,j,k,ncolors; 22 DM dm; 23 PetscBool flg; 24 ISColoring coloring; 25 MatColoring mc; 26 Vec W,G,F; 27 PetscScalar h=gs->h; 28 IS *coloris; 29 PetscScalar f,g,x,w,d; 30 PetscReal dxt,xt,ft,ft1; 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 PetscBool mat = gs->secant_mat,equal; 39 40 PetscFunctionBegin; 41 if (snes->nwork < 3) { 42 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr); 43 } 44 W = snes->work[0]; 45 G = snes->work[1]; 46 F = snes->work[2]; 47 ierr = SNESGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr); 48 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 49 ierr = SNESGetFunction(snes,NULL,&func,&fctx);CHKERRQ(ierr); 50 ierr = PetscObjectQuery((PetscObject)snes,"SNESGSColoring",(PetscObject*)&colorcontainer);CHKERRQ(ierr); 51 if (!colorcontainer) { 52 /* create the coloring */ 53 ierr = DMHasColoring(dm,&flg);CHKERRQ(ierr); 54 if (flg && !mat) { 55 ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);CHKERRQ(ierr); 56 } else { 57 if (!snes->jacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);} 58 ierr = MatColoringCreate(snes->jacobian,&mc);CHKERRQ(ierr); 59 ierr = MatColoringSetDistance(mc,1);CHKERRQ(ierr); 60 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 61 ierr = MatColoringApply(mc,&coloring);CHKERRQ(ierr); 62 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 63 } 64 ierr = PetscContainerCreate(PetscObjectComm((PetscObject)snes),&colorcontainer);CHKERRQ(ierr); 65 ierr = PetscContainerSetPointer(colorcontainer,(void *)coloring);CHKERRQ(ierr); 66 ierr = PetscContainerSetUserDestroy(colorcontainer,(PetscErrorCode (*)(void *))GSDestroy_Private);CHKERRQ(ierr); 67 ierr = PetscObjectCompose((PetscObject)snes,"SNESGSColoring",(PetscObject)colorcontainer);CHKERRQ(ierr); 68 ierr = PetscContainerDestroy(&colorcontainer);CHKERRQ(ierr); 69 } else { 70 ierr = PetscContainerGetPointer(colorcontainer,(void **)&coloring);CHKERRQ(ierr); 71 } 72 ierr = ISColoringGetIS(coloring,&ncolors,&coloris);CHKERRQ(ierr); 73 ierr = VecEqual(X,snes->vec_sol,&equal);CHKERRQ(ierr); 74 if (equal) { 75 /* assume that the function is already computed */ 76 ierr = VecCopy(snes->vec_func,F);CHKERRQ(ierr); 77 } else { 78 ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 79 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr); 80 ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 81 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);} 82 } 83 for (i=0;i<ncolors;i++) { 84 ierr = ISGetIndices(coloris[i],&idx);CHKERRQ(ierr); 85 ierr = ISGetLocalSize(coloris[i],&size);CHKERRQ(ierr); 86 ierr = VecCopy(X,W);CHKERRQ(ierr); 87 for (j=0;j<size;j++) { 88 ierr = VecSetValue(W,idx[j],h,ADD_VALUES);CHKERRQ(ierr); 89 } 90 ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 91 ierr = (*func)(snes,W,G,fctx);CHKERRQ(ierr); 92 ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 93 if (B) {ierr = VecAXPY(G,-1.0,B);CHKERRQ(ierr);} 94 for (k=0;k<its;k++) { 95 dxt = 0.; 96 xt = 0.; 97 for (j=0;j<size;j++) { 98 ierr = VecGetValues(F,1,&idx[j],&f);CHKERRQ(ierr); 99 ierr = VecGetValues(X,1,&idx[j],&x);CHKERRQ(ierr); 100 ierr = VecGetValues(G,1,&idx[j],&g);CHKERRQ(ierr); 101 ierr = VecGetValues(W,1,&idx[j],&w);CHKERRQ(ierr); 102 if (PetscAbsScalar(g-f) > atol) { 103 d = (x*g-w*f) / PetscRealPart(g-f); 104 } else { 105 d = x; 106 } 107 dxt += (d-x)*(d-x); 108 xt += x*x; 109 ft += f*f; 110 ierr = VecSetValue(X,idx[j],d,INSERT_VALUES);CHKERRQ(ierr); 111 } 112 if (k == 0) ft1 = PetscSqrtScalar(ft); 113 if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) break; 114 if (PetscSqrtReal(ft) < atol) break; 115 if (rtol*ft1 > PetscSqrtReal(ft)) break; 116 if (i < ncolors-1 || k < its-1) { 117 ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 118 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr); 119 ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr); 120 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);} 121 } 122 if (k<its-1) { 123 ierr = VecSwap(X,W);CHKERRQ(ierr); 124 ierr = VecSwap(F,G);CHKERRQ(ierr); 125 } 126 } 127 } 128 ierr = ISColoringRestoreIS(coloring,&coloris);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131