1 #include <../src/snes/impls/gs/gsimpl.h> 2 3 PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx) 4 { 5 SNES_NGS *gs = (SNES_NGS*)snes->data; 6 PetscInt i,j,k,ncolors; 7 DM dm; 8 PetscBool flg; 9 ISColoring coloring = gs->coloring; 10 MatColoring mc; 11 Vec W,G,F; 12 PetscScalar h=gs->h; 13 IS *coloris; 14 PetscScalar f,g,x,w,d; 15 PetscReal dxt,xt,ft,ft1=0; 16 const PetscInt *idx; 17 PetscInt size,s; 18 PetscReal atol,rtol,stol; 19 PetscInt its; 20 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 21 void *fctx; 22 PetscBool mat = gs->secant_mat,equal,isdone,alldone; 23 PetscScalar *xa,*wa; 24 const PetscScalar *fa,*ga; 25 26 PetscFunctionBegin; 27 if (snes->nwork < 3) { 28 CHKERRQ(SNESSetWorkVecs(snes,3)); 29 } 30 W = snes->work[0]; 31 G = snes->work[1]; 32 F = snes->work[2]; 33 CHKERRQ(VecGetOwnershipRange(X,&s,NULL)); 34 CHKERRQ(SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its)); 35 CHKERRQ(SNESGetDM(snes,&dm)); 36 CHKERRQ(SNESGetFunction(snes,NULL,&func,&fctx)); 37 if (!coloring) { 38 /* create the coloring */ 39 CHKERRQ(DMHasColoring(dm,&flg)); 40 if (flg && !mat) { 41 CHKERRQ(DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring)); 42 } else { 43 if (!snes->jacobian) CHKERRQ(SNESSetUpMatrices(snes)); 44 CHKERRQ(MatColoringCreate(snes->jacobian,&mc)); 45 CHKERRQ(MatColoringSetDistance(mc,1)); 46 CHKERRQ(MatColoringSetFromOptions(mc)); 47 CHKERRQ(MatColoringApply(mc,&coloring)); 48 CHKERRQ(MatColoringDestroy(&mc)); 49 } 50 gs->coloring = coloring; 51 } 52 CHKERRQ(ISColoringGetIS(coloring,PETSC_USE_POINTER,&ncolors,&coloris)); 53 CHKERRQ(VecEqual(X,snes->vec_sol,&equal)); 54 if (equal && snes->normschedule == SNES_NORM_ALWAYS) { 55 /* assume that the function is already computed */ 56 CHKERRQ(VecCopy(snes->vec_func,F)); 57 } else { 58 CHKERRQ(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 59 CHKERRQ((*func)(snes,X,F,fctx)); 60 CHKERRQ(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 61 if (B) CHKERRQ(VecAXPY(F,-1.0,B)); 62 } 63 for (i=0;i<ncolors;i++) { 64 CHKERRQ(ISGetIndices(coloris[i],&idx)); 65 CHKERRQ(ISGetLocalSize(coloris[i],&size)); 66 CHKERRQ(VecCopy(X,W)); 67 CHKERRQ(VecGetArray(W,&wa)); 68 for (j=0;j<size;j++) { 69 wa[idx[j]-s] += h; 70 } 71 CHKERRQ(VecRestoreArray(W,&wa)); 72 CHKERRQ(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 73 CHKERRQ((*func)(snes,W,G,fctx)); 74 CHKERRQ(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 75 if (B) CHKERRQ(VecAXPY(G,-1.0,B)); 76 for (k=0;k<its;k++) { 77 dxt = 0.; 78 xt = 0.; 79 ft = 0.; 80 CHKERRQ(VecGetArray(W,&wa)); 81 CHKERRQ(VecGetArray(X,&xa)); 82 CHKERRQ(VecGetArrayRead(F,&fa)); 83 CHKERRQ(VecGetArrayRead(G,&ga)); 84 for (j=0;j<size;j++) { 85 f = fa[idx[j]-s]; 86 x = xa[idx[j]-s]; 87 g = ga[idx[j]-s]; 88 w = wa[idx[j]-s]; 89 if (PetscAbsScalar(g-f) > atol) { 90 /* This is equivalent to d = x - (h*f) / PetscRealPart(g-f) */ 91 d = (x*g-w*f) / PetscRealPart(g-f); 92 } else { 93 d = x; 94 } 95 dxt += PetscRealPart(PetscSqr(d-x)); 96 xt += PetscRealPart(PetscSqr(x)); 97 ft += PetscRealPart(PetscSqr(f)); 98 xa[idx[j]-s] = d; 99 } 100 CHKERRQ(VecRestoreArray(X,&xa)); 101 CHKERRQ(VecRestoreArrayRead(F,&fa)); 102 CHKERRQ(VecRestoreArrayRead(G,&ga)); 103 CHKERRQ(VecRestoreArray(W,&wa)); 104 105 if (k == 0) ft1 = PetscSqrtReal(ft); 106 if (k<its-1) { 107 isdone = PETSC_FALSE; 108 if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE; 109 if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE; 110 if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE; 111 CHKERRMPI(MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes))); 112 if (alldone) break; 113 } 114 if (i < ncolors-1 || k < its-1) { 115 CHKERRQ(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 116 CHKERRQ((*func)(snes,X,F,fctx)); 117 CHKERRQ(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 118 if (B) CHKERRQ(VecAXPY(F,-1.0,B)); 119 } 120 if (k<its-1) { 121 CHKERRQ(VecSwap(X,W)); 122 CHKERRQ(VecSwap(F,G)); 123 } 124 } 125 } 126 CHKERRQ(ISColoringRestoreIS(coloring,PETSC_USE_POINTER,&coloris)); 127 PetscFunctionReturn(0); 128 } 129