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