1 #include <../src/snes/impls/gs/gsimpl.h>
2
SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,PetscCtx ctx)3 PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes, Vec X, Vec B, PetscCtx 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 SNESFunctionFn *func;
21 void *fctx;
22 PetscBool mat = gs->secant_mat, equal, isdone, alldone;
23 PetscScalar *xa, *wa;
24 const PetscScalar *fa, *ga;
25
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));
100
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 PetscCallMPI(MPIU_Allreduce(&isdone, &alldone, 1, MPI_C_BOOL, MPI_LAND, 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 }
125