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