Actual source code: gssecant.c

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

  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;

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

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