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 PetscCall(SNESSetWorkVecs(snes,3)); 29 } 30 W = snes->work[0]; 31 G = snes->work[1]; 32 F = snes->work[2]; 33 PetscCall(VecGetOwnershipRange(X,&s,NULL)); 34 PetscCall(SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its)); 35 PetscCall(SNESGetDM(snes,&dm)); 36 PetscCall(SNESGetFunction(snes,NULL,&func,&fctx)); 37 if (!coloring) { 38 /* create the coloring */ 39 PetscCall(DMHasColoring(dm,&flg)); 40 if (flg && !mat) { 41 PetscCall(DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring)); 42 } else { 43 if (!snes->jacobian) PetscCall(SNESSetUpMatrices(snes)); 44 PetscCall(MatColoringCreate(snes->jacobian,&mc)); 45 PetscCall(MatColoringSetDistance(mc,1)); 46 PetscCall(MatColoringSetFromOptions(mc)); 47 PetscCall(MatColoringApply(mc,&coloring)); 48 PetscCall(MatColoringDestroy(&mc)); 49 } 50 gs->coloring = coloring; 51 } 52 PetscCall(ISColoringGetIS(coloring,PETSC_USE_POINTER,&ncolors,&coloris)); 53 PetscCall(VecEqual(X,snes->vec_sol,&equal)); 54 if (equal && snes->normschedule == SNES_NORM_ALWAYS) { 55 /* assume that the function is already computed */ 56 PetscCall(VecCopy(snes->vec_func,F)); 57 } else { 58 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 59 PetscCall((*func)(snes,X,F,fctx)); 60 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 61 if (B) PetscCall(VecAXPY(F,-1.0,B)); 62 } 63 for (i=0;i<ncolors;i++) { 64 PetscCall(ISGetIndices(coloris[i],&idx)); 65 PetscCall(ISGetLocalSize(coloris[i],&size)); 66 PetscCall(VecCopy(X,W)); 67 PetscCall(VecGetArray(W,&wa)); 68 for (j=0;j<size;j++) { 69 wa[idx[j]-s] += h; 70 } 71 PetscCall(VecRestoreArray(W,&wa)); 72 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 73 PetscCall((*func)(snes,W,G,fctx)); 74 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 75 if (B) PetscCall(VecAXPY(G,-1.0,B)); 76 for (k=0;k<its;k++) { 77 dxt = 0.; 78 xt = 0.; 79 ft = 0.; 80 PetscCall(VecGetArray(W,&wa)); 81 PetscCall(VecGetArray(X,&xa)); 82 PetscCall(VecGetArrayRead(F,&fa)); 83 PetscCall(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 PetscCall(VecRestoreArray(X,&xa)); 101 PetscCall(VecRestoreArrayRead(F,&fa)); 102 PetscCall(VecRestoreArrayRead(G,&ga)); 103 PetscCall(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 PetscCall(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 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0)); 116 PetscCall((*func)(snes,X,F,fctx)); 117 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0)); 118 if (B) PetscCall(VecAXPY(F,-1.0,B)); 119 } 120 if (k<its-1) { 121 PetscCall(VecSwap(X,W)); 122 PetscCall(VecSwap(F,G)); 123 } 124 } 125 } 126 PetscCall(ISColoringRestoreIS(coloring,PETSC_USE_POINTER,&coloris)); 127 PetscFunctionReturn(0); 128 } 129