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:   PetscErrorCode    ierr;
      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,*wa;
     25:   const PetscScalar *fa,*ga;

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

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