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

  5: static PetscErrorCode SNESNGSDestroy_Private(ISColoring coloring)
  6: {

 10:   ISColoringDestroy(&coloring);
 11:   return(0);
 12: }

 16: PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx)
 17: {
 19:   SNES_NGS       *gs = (SNES_NGS*)snes->data;
 20:   PetscInt       i,j,k,ncolors;
 21:   DM             dm;
 22:   PetscBool      flg;
 23:   ISColoring     coloring;
 24:   MatColoring    mc;
 25:   Vec            W,G,F;
 26:   PetscScalar    h=gs->h;
 27:   IS             *coloris;
 28:   PetscScalar    f,g,x,w,d;
 29:   PetscReal      dxt,xt,ft,ft1=0;
 30:   const PetscInt *idx;
 31:   PetscInt       size,s;
 32:   PetscReal      atol,rtol,stol;
 33:   PetscInt       its;
 34:   PetscErrorCode (*func)(SNES,Vec,Vec,void*);
 35:   void           *fctx;
 36:   PetscContainer colorcontainer;
 37:   PetscBool      mat = gs->secant_mat,equal,isdone,alldone;
 38:   PetscScalar    *xa,*fa,*wa,*ga;

 41:   if (snes->nwork < 3) {
 42:     SNESSetWorkVecs(snes,3);
 43:   }
 44:   W = snes->work[0];
 45:   G = snes->work[1];
 46:   F = snes->work[2];
 47:   VecGetOwnershipRange(X,&s,NULL);
 48:   SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
 49:   SNESGetDM(snes,&dm);
 50:   SNESGetFunction(snes,NULL,&func,&fctx);
 51:   PetscObjectQuery((PetscObject)snes,"SNESNGSColoring",(PetscObject*)&colorcontainer);
 52:   if (!colorcontainer) {
 53:     /* create the coloring */
 54:     DMHasColoring(dm,&flg);
 55:     if (flg && !mat) {
 56:       DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);
 57:     } else {
 58:       if (!snes->jacobian) {SNESSetUpMatrices(snes);}
 59:       MatColoringCreate(snes->jacobian,&mc);
 60:       MatColoringSetDistance(mc,1);
 61:       MatColoringSetFromOptions(mc);
 62:       MatColoringApply(mc,&coloring);
 63:       MatColoringDestroy(&mc);
 64:     }
 65:     PetscContainerCreate(PetscObjectComm((PetscObject)snes),&colorcontainer);
 66:     PetscContainerSetPointer(colorcontainer,(void *)coloring);
 67:     PetscContainerSetUserDestroy(colorcontainer,(PetscErrorCode (*)(void *))SNESNGSDestroy_Private);
 68:     PetscObjectCompose((PetscObject)snes,"SNESNGSColoring",(PetscObject)colorcontainer);
 69:     PetscContainerDestroy(&colorcontainer);
 70:   } else {
 71:     PetscContainerGetPointer(colorcontainer,(void **)&coloring);
 72:   }
 73:   ISColoringGetIS(coloring,&ncolors,&coloris);
 74:   VecEqual(X,snes->vec_sol,&equal);
 75:   if (equal && snes->normschedule == SNES_NORM_ALWAYS) {
 76:     /* assume that the function is already computed */
 77:     VecCopy(snes->vec_func,F);
 78:   } else {
 79:     PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
 80:     (*func)(snes,X,F,fctx);
 81:     PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
 82:     if (B) {VecAXPY(F,-1.0,B);}
 83:   }
 84:   VecGetArray(X,&xa);
 85:   VecGetArray(F,&fa);
 86:   VecGetArray(G,&ga);
 87:   VecGetArray(W,&wa);
 88:   for (i=0;i<ncolors;i++) {
 89:     ISGetIndices(coloris[i],&idx);
 90:     ISGetLocalSize(coloris[i],&size);
 91:     VecCopy(X,W);
 92:     for (j=0;j<size;j++) {
 93:       wa[idx[j]-s] += h;
 94:     }
 95:     PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
 96:     (*func)(snes,W,G,fctx);
 97:     PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
 98:     if (B) {VecAXPY(G,-1.0,B);}
 99:     for (k=0;k<its;k++) {
100:       dxt = 0.;
101:       xt = 0.;
102:       ft = 0.;
103:       for (j=0;j<size;j++) {
104:         f = fa[idx[j]-s];
105:         x = xa[idx[j]-s];
106:         g = ga[idx[j]-s];
107:         w = wa[idx[j]-s];
108:         if (PetscAbsScalar(g-f) > atol) {
109:           /* This is equivalent to d = x - (h*f) / PetscRealPart(g-f) */
110:           d = (x*g-w*f) / PetscRealPart(g-f);
111:         } else {
112:           d = x;
113:         }
114:         dxt += PetscRealPart(PetscSqr(d-x));
115:         xt += PetscRealPart(PetscSqr(x));
116:         ft += PetscRealPart(PetscSqr(f));
117:         xa[idx[j]-s] = d;
118:       }

120:       if (k == 0) ft1 = PetscSqrtReal(ft);
121:       if (k<its-1) {
122:         isdone = PETSC_FALSE;
123:         if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE;
124:         if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE;
125:         if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE;
126:         MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));
127:         if (alldone) break;
128:       }
129:       if (i < ncolors-1 || k < its-1) {
130:         PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
131:         (*func)(snes,X,F,fctx);
132:         PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
133:         if (B) {VecAXPY(F,-1.0,B);}
134:       }
135:       if (k<its-1) {
136:         VecSwap(X,W);
137:         VecSwap(F,G);
138:       }
139:     }
140:   }
141:   VecRestoreArray(X,&xa);
142:   VecRestoreArray(F,&fa);
143:   VecRestoreArray(G,&ga);
144:   VecRestoreArray(W,&wa);
145:   ISColoringRestoreIS(coloring,&coloris);
146:   return(0);
147: }