Actual source code: gssecant.c
petsc-3.14.6 2021-03-30
1: #include <../src/snes/impls/gs/gsimpl.h>
3: PETSC_EXTERN PetscErrorCode SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx)
4: {
5: PetscErrorCode ierr;
6: SNES_NGS *gs = (SNES_NGS*)snes->data;
7: PetscInt i,j,k,ncolors;
8: DM dm;
9: PetscBool flg;
10: ISColoring coloring = gs->coloring;
11: MatColoring mc;
12: Vec W,G,F;
13: PetscScalar h=gs->h;
14: IS *coloris;
15: PetscScalar f,g,x,w,d;
16: PetscReal dxt,xt,ft,ft1=0;
17: const PetscInt *idx;
18: PetscInt size,s;
19: PetscReal atol,rtol,stol;
20: PetscInt its;
21: PetscErrorCode (*func)(SNES,Vec,Vec,void*);
22: void *fctx;
23: PetscBool mat = gs->secant_mat,equal,isdone,alldone;
24: PetscScalar *xa,*wa;
25: const PetscScalar *fa,*ga;
28: if (snes->nwork < 3) {
29: SNESSetWorkVecs(snes,3);
30: }
31: W = snes->work[0];
32: G = snes->work[1];
33: F = snes->work[2];
34: VecGetOwnershipRange(X,&s,NULL);
35: SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);
36: SNESGetDM(snes,&dm);
37: SNESGetFunction(snes,NULL,&func,&fctx);
38: if (!coloring) {
39: /* create the coloring */
40: DMHasColoring(dm,&flg);
41: if (flg && !mat) {
42: DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);
43: } else {
44: if (!snes->jacobian) {SNESSetUpMatrices(snes);}
45: MatColoringCreate(snes->jacobian,&mc);
46: MatColoringSetDistance(mc,1);
47: MatColoringSetFromOptions(mc);
48: MatColoringApply(mc,&coloring);
49: MatColoringDestroy(&mc);
50: }
51: gs->coloring = coloring;
52: }
53: ISColoringGetIS(coloring,PETSC_USE_POINTER,&ncolors,&coloris);
54: VecEqual(X,snes->vec_sol,&equal);
55: if (equal && snes->normschedule == SNES_NORM_ALWAYS) {
56: /* assume that the function is already computed */
57: VecCopy(snes->vec_func,F);
58: } else {
59: PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
60: (*func)(snes,X,F,fctx);
61: PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
62: if (B) {VecAXPY(F,-1.0,B);}
63: }
64: for (i=0;i<ncolors;i++) {
65: ISGetIndices(coloris[i],&idx);
66: ISGetLocalSize(coloris[i],&size);
67: VecCopy(X,W);
68: VecGetArray(W,&wa);
69: for (j=0;j<size;j++) {
70: wa[idx[j]-s] += h;
71: }
72: VecRestoreArray(W,&wa);
73: PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
74: (*func)(snes,W,G,fctx);
75: PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
76: if (B) {VecAXPY(G,-1.0,B);}
77: for (k=0;k<its;k++) {
78: dxt = 0.;
79: xt = 0.;
80: ft = 0.;
81: VecGetArray(W,&wa);
82: VecGetArray(X,&xa);
83: VecGetArrayRead(F,&fa);
84: VecGetArrayRead(G,&ga);
85: for (j=0;j<size;j++) {
86: f = fa[idx[j]-s];
87: x = xa[idx[j]-s];
88: g = ga[idx[j]-s];
89: w = wa[idx[j]-s];
90: if (PetscAbsScalar(g-f) > atol) {
91: /* This is equivalent to d = x - (h*f) / PetscRealPart(g-f) */
92: d = (x*g-w*f) / PetscRealPart(g-f);
93: } else {
94: d = x;
95: }
96: dxt += PetscRealPart(PetscSqr(d-x));
97: xt += PetscRealPart(PetscSqr(x));
98: ft += PetscRealPart(PetscSqr(f));
99: xa[idx[j]-s] = d;
100: }
101: VecRestoreArray(X,&xa);
102: VecRestoreArrayRead(F,&fa);
103: VecRestoreArrayRead(G,&ga);
104: VecRestoreArray(W,&wa);
106: if (k == 0) ft1 = PetscSqrtReal(ft);
107: if (k<its-1) {
108: isdone = PETSC_FALSE;
109: if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE;
110: if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE;
111: if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE;
112: MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));
113: if (alldone) break;
114: }
115: if (i < ncolors-1 || k < its-1) {
116: PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);
117: (*func)(snes,X,F,fctx);
118: PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);
119: if (B) {VecAXPY(F,-1.0,B);}
120: }
121: if (k<its-1) {
122: VecSwap(X,W);
123: VecSwap(F,G);
124: }
125: }
126: }
127: ISColoringRestoreIS(coloring,PETSC_USE_POINTER,&coloris);
128: return(0);
129: }