1 #include <../src/snes/impls/gs/gsimpl.h> 2 3 PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes, Vec X, Vec B, void *ctx) { 4 SNES_NGS *gs = (SNES_NGS *)snes->data; 5 PetscInt i, j, k, ncolors; 6 DM dm; 7 PetscBool flg; 8 ISColoring coloring = gs->coloring; 9 MatColoring mc; 10 Vec W, G, F; 11 PetscScalar h = gs->h; 12 IS *coloris; 13 PetscScalar f, g, x, w, d; 14 PetscReal dxt, xt, ft, ft1 = 0; 15 const PetscInt *idx; 16 PetscInt size, s; 17 PetscReal atol, rtol, stol; 18 PetscInt its; 19 PetscErrorCode (*func)(SNES, Vec, Vec, void *); 20 void *fctx; 21 PetscBool mat = gs->secant_mat, equal, isdone, alldone; 22 PetscScalar *xa, *wa; 23 const PetscScalar *fa, *ga; 24 25 PetscFunctionBegin; 26 if (snes->nwork < 3) { PetscCall(SNESSetWorkVecs(snes, 3)); } 27 W = snes->work[0]; 28 G = snes->work[1]; 29 F = snes->work[2]; 30 PetscCall(VecGetOwnershipRange(X, &s, NULL)); 31 PetscCall(SNESNGSGetTolerances(snes, &atol, &rtol, &stol, &its)); 32 PetscCall(SNESGetDM(snes, &dm)); 33 PetscCall(SNESGetFunction(snes, NULL, &func, &fctx)); 34 if (!coloring) { 35 /* create the coloring */ 36 PetscCall(DMHasColoring(dm, &flg)); 37 if (flg && !mat) { 38 PetscCall(DMCreateColoring(dm, IS_COLORING_GLOBAL, &coloring)); 39 } else { 40 if (!snes->jacobian) PetscCall(SNESSetUpMatrices(snes)); 41 PetscCall(MatColoringCreate(snes->jacobian, &mc)); 42 PetscCall(MatColoringSetDistance(mc, 1)); 43 PetscCall(MatColoringSetFromOptions(mc)); 44 PetscCall(MatColoringApply(mc, &coloring)); 45 PetscCall(MatColoringDestroy(&mc)); 46 } 47 gs->coloring = coloring; 48 } 49 PetscCall(ISColoringGetIS(coloring, PETSC_USE_POINTER, &ncolors, &coloris)); 50 PetscCall(VecEqual(X, snes->vec_sol, &equal)); 51 if (equal && snes->normschedule == SNES_NORM_ALWAYS) { 52 /* assume that the function is already computed */ 53 PetscCall(VecCopy(snes->vec_func, F)); 54 } else { 55 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval, snes, X, B, 0)); 56 PetscCall((*func)(snes, X, F, fctx)); 57 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval, snes, X, B, 0)); 58 if (B) PetscCall(VecAXPY(F, -1.0, B)); 59 } 60 for (i = 0; i < ncolors; i++) { 61 PetscCall(ISGetIndices(coloris[i], &idx)); 62 PetscCall(ISGetLocalSize(coloris[i], &size)); 63 PetscCall(VecCopy(X, W)); 64 PetscCall(VecGetArray(W, &wa)); 65 for (j = 0; j < size; j++) { wa[idx[j] - s] += h; } 66 PetscCall(VecRestoreArray(W, &wa)); 67 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval, snes, X, B, 0)); 68 PetscCall((*func)(snes, W, G, fctx)); 69 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval, snes, X, B, 0)); 70 if (B) PetscCall(VecAXPY(G, -1.0, B)); 71 for (k = 0; k < its; k++) { 72 dxt = 0.; 73 xt = 0.; 74 ft = 0.; 75 PetscCall(VecGetArray(W, &wa)); 76 PetscCall(VecGetArray(X, &xa)); 77 PetscCall(VecGetArrayRead(F, &fa)); 78 PetscCall(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 PetscCall(VecRestoreArray(X, &xa)); 96 PetscCall(VecRestoreArrayRead(F, &fa)); 97 PetscCall(VecRestoreArrayRead(G, &ga)); 98 PetscCall(VecRestoreArray(W, &wa)); 99 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 PetscCall(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 PetscCall(PetscLogEventBegin(SNES_NGSFuncEval, snes, X, B, 0)); 111 PetscCall((*func)(snes, X, F, fctx)); 112 PetscCall(PetscLogEventEnd(SNES_NGSFuncEval, snes, X, B, 0)); 113 if (B) PetscCall(VecAXPY(F, -1.0, B)); 114 } 115 if (k < its - 1) { 116 PetscCall(VecSwap(X, W)); 117 PetscCall(VecSwap(F, G)); 118 } 119 } 120 } 121 PetscCall(ISColoringRestoreIS(coloring, PETSC_USE_POINTER, &coloris)); 122 PetscFunctionReturn(0); 123 } 124