1:  #include <../src/snes/impls/gs/gsimpl.h>

  3: PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx)
  4: {
  6:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
  7:   PetscInt       i,j,k,ncolors;
  8:   DM             dm;
  9:   PetscBool      flg;
 10:   ISColoring     coloring = gs->coloring;
 11:   MatColoring    mc;
 12:   Vec            W,G,F;
 13:   PetscScalar    h=gs->h;
 14:   IS             *coloris;
 15:   PetscScalar    f,g,x,w,d;
 16:   PetscReal      dxt,xt,ft,ft1=0;
 17:   const PetscInt *idx;
 18:   PetscInt       size,s;
 19:   PetscReal      atol,rtol,stol;
 20:   PetscInt       its;
 21:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
 22:   void           *fctx;
 23:   PetscBool      mat = gs->secant_mat,equal,isdone,alldone;
 24:   PetscScalar    *xa,*fa,*wa,*ga;

 27:   if (snes->nwork < 3) {
 28:     SNESSetWorkVecs(snes,3);
 29:   }
 30:   W = snes->work[0];
 31:   G = snes->work[1];
 32:   F = snes->work[2];
 33:   VecGetOwnershipRange(X,&s,NULL);
 34:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
 35:   SNESGetDM(snes,&dm);
 36:   SNESGetFunction(snes,NULL,&func,&fctx);
 37:   if (!coloring) {
 38:     /* create the coloring */
 39:     DMHasColoring(dm,&flg);
 40:     if (flg && !mat) {
 41:       DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);
 42:     } else {
 43:       if (!snes->jacobian) {SNESSetUpMatrices(snes);}
 44:       MatColoringCreate(snes->jacobian,&mc);
 45:       MatColoringSetDistance(mc,1);
 46:       MatColoringSetFromOptions(mc);
 47:       MatColoringApply(mc,&coloring);
 48:       MatColoringDestroy(&mc);
 49:     }
 50:     gs->coloring = coloring;
 51:   }
 52:   ISColoringGetIS(coloring,&ncolors,&coloris);
 53:   VecEqual(X,snes->vec_sol,&equal);
 54:   if (equal && snes->normschedule == SNES_NORM_ALWAYS) {
 55:     /* assume that the function is already computed */
 56:     VecCopy(snes->vec_func,F);
 57:   } else {
 58:     PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
 59:     (*func)(snes,X,F,fctx);
 60:     PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
 61:     if (B) {VecAXPY(F,-1.0,B);}
 62:   }
 63:   VecGetArray(X,&xa);
 64:   VecGetArray(F,&fa);
 65:   VecGetArray(G,&ga);
 66:   VecGetArray(W,&wa);
 67:   for (i=0;i<ncolors;i++) {
 68:     ISGetIndices(coloris[i],&idx);
 69:     ISGetLocalSize(coloris[i],&size);
 70:     VecCopy(X,W);
 71:     for (j=0;j<size;j++) {
 72:       wa[idx[j]-s] += h;
 73:     }
 74:     PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
 75:     (*func)(snes,W,G,fctx);
 76:     PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
 77:     if (B) {VecAXPY(G,-1.0,B);}
 78:     for (k=0;k<its;k++) {
 79:       dxt = 0.;
 80:       xt = 0.;
 81:       ft = 0.;
 82:       for (j=0;j<size;j++) {
 83:         f = fa[idx[j]-s];
 84:         x = xa[idx[j]-s];
 85:         g = ga[idx[j]-s];
 86:         w = wa[idx[j]-s];
 87:         if (PetscAbsScalar(g-f) > atol) {
 88:           /* This is equivalent to d = x - (h*f) / PetscRealPart(g-f) */
 89:           d = (x*g-w*f) / PetscRealPart(g-f);
 90:         } else {
 91:           d = x;
 92:         }
 93:         dxt += PetscRealPart(PetscSqr(d-x));
 94:         xt += PetscRealPart(PetscSqr(x));
 95:         ft += PetscRealPart(PetscSqr(f));
 96:         xa[idx[j]-s] = d;
 97:       }

 99:       if (k == 0) ft1 = PetscSqrtReal(ft);
100:       if (k<its-1) {
101:         isdone = PETSC_FALSE;
102:         if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE;
103:         if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE;
104:         if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE;
105:         MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));
106:         if (alldone) break;
107:       }
108:       if (i < ncolors-1 || k < its-1) {
109:         PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
110:         (*func)(snes,X,F,fctx);
111:         PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
112:         if (B) {VecAXPY(F,-1.0,B);}
113:       }
114:       if (k<its-1) {
115:         VecSwap(X,W);
116:         VecSwap(F,G);
117:       }
118:     }
119:   }
120:   VecRestoreArray(X,&xa);
121:   VecRestoreArray(F,&fa);
122:   VecRestoreArray(G,&ga);
123:   VecRestoreArray(W,&wa);
124:   ISColoringRestoreIS(coloring,&coloris);
125:   return(0);
126: }