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

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