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: }