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

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