1 #include <../src/snes/impls/gs/gsimpl.h> 2 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 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