1 #include <../src/snes/impls/gs/gsimpl.h>
2
SNESComputeNGSDefaultSecant(SNES snes,Vec X,Vec B,void * ctx)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;
26
27 PetscFunctionBegin;
28 if (snes->nwork < 3) {
29 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
30 }
31 W = snes->work[0];
32 G = snes->work[1];
33 F = snes->work[2];
34 ierr = VecGetOwnershipRange(X,&s,NULL);CHKERRQ(ierr);
35 ierr = SNESNGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr);
36 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
37 ierr = SNESGetFunction(snes,NULL,&func,&fctx);CHKERRQ(ierr);
38 if (!coloring) {
39 /* create the coloring */
40 ierr = DMHasColoring(dm,&flg);CHKERRQ(ierr);
41 if (flg && !mat) {
42 ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);CHKERRQ(ierr);
43 } else {
44 if (!snes->jacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);}
45 ierr = MatColoringCreate(snes->jacobian,&mc);CHKERRQ(ierr);
46 ierr = MatColoringSetDistance(mc,1);CHKERRQ(ierr);
47 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
48 ierr = MatColoringApply(mc,&coloring);CHKERRQ(ierr);
49 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
50 }
51 gs->coloring = coloring;
52 }
53 ierr = ISColoringGetIS(coloring,PETSC_USE_POINTER,&ncolors,&coloris);CHKERRQ(ierr);
54 ierr = VecEqual(X,snes->vec_sol,&equal);CHKERRQ(ierr);
55 if (equal && snes->normschedule == SNES_NORM_ALWAYS) {
56 /* assume that the function is already computed */
57 ierr = VecCopy(snes->vec_func,F);CHKERRQ(ierr);
58 } else {
59 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
60 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr);
61 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
62 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);}
63 }
64 for (i=0;i<ncolors;i++) {
65 ierr = ISGetIndices(coloris[i],&idx);CHKERRQ(ierr);
66 ierr = ISGetLocalSize(coloris[i],&size);CHKERRQ(ierr);
67 ierr = VecCopy(X,W);CHKERRQ(ierr);
68 ierr = VecGetArray(W,&wa);CHKERRQ(ierr);
69 for (j=0;j<size;j++) {
70 wa[idx[j]-s] += h;
71 }
72 ierr = VecRestoreArray(W,&wa);CHKERRQ(ierr);
73 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
74 ierr = (*func)(snes,W,G,fctx);CHKERRQ(ierr);
75 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
76 if (B) {ierr = VecAXPY(G,-1.0,B);CHKERRQ(ierr);}
77 for (k=0;k<its;k++) {
78 dxt = 0.;
79 xt = 0.;
80 ft = 0.;
81 ierr = VecGetArray(W,&wa);CHKERRQ(ierr);
82 ierr = VecGetArray(X,&xa);CHKERRQ(ierr);
83 ierr = VecGetArrayRead(F,&fa);CHKERRQ(ierr);
84 ierr = VecGetArrayRead(G,&ga);CHKERRQ(ierr);
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 ierr = VecRestoreArray(X,&xa);CHKERRQ(ierr);
102 ierr = VecRestoreArrayRead(F,&fa);CHKERRQ(ierr);
103 ierr = VecRestoreArrayRead(G,&ga);CHKERRQ(ierr);
104 ierr = VecRestoreArray(W,&wa);CHKERRQ(ierr);
105
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 ierr = MPIU_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
113 if (alldone) break;
114 }
115 if (i < ncolors-1 || k < its-1) {
116 ierr = PetscLogEventBegin(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
117 ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr);
118 ierr = PetscLogEventEnd(SNES_NGSFuncEval,snes,X,B,0);CHKERRQ(ierr);
119 if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);}
120 }
121 if (k<its-1) {
122 ierr = VecSwap(X,W);CHKERRQ(ierr);
123 ierr = VecSwap(F,G);CHKERRQ(ierr);
124 }
125 }
126 }
127 ierr = ISColoringRestoreIS(coloring,PETSC_USE_POINTER,&coloris);CHKERRQ(ierr);
128 PetscFunctionReturn(0);
129 }
130