1 #ifndef lcfiplus_algo_h
2 #define lcfiplus_algo_h 1
15 const double R=0.61803399;
16 const double C=(1.0-
R);
17 const double GOLD=1.618034;
25 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
26 #define SHFT2(a,b,c) (a)=(b);(b)=(c);
27 #define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
28 #define SIGN(a,b) ((b)>=0.0 ? fabs(a) : -fabs(a))
29 #define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
32 double golden(
double ax,
double bx,
double cx, T* obj,
double tol,
double& xmin) {
33 double f1,f2,x0,x1,x2,x3;
36 if (fabs(cx-bx) > fabs(bx-ax)) {
45 while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) {
48 SHFT2(f1,f2,(*obj)(&x2));
51 SHFT2(f2,f1,(*obj)(&x1));
65 double brent(
double ax,
double bx,
double cx, T* obj,
double tol,
double& xmin) {
67 double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
70 a=(ax < cx ? ax : cx);
71 b=(ax > cx ? ax : cx);
74 for (iter=1; iter<=
ITMAX; iter++) {
76 tol2=2.0*(tol1=tol*fabs(x)+
ZEPS);
77 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
90 if (fabs(p) >= fabs(0.5*q*etemp) || p<= q*(a-x) || p >= q*(b-x))
91 d=
CGOLD*(e=(x >= xm ? a-x : b-x));
95 if (u-a < tol2 || b-u < tol2)
99 d=
CGOLD*(e=(x >= xm ? a-x : b-x));
101 u=(fabs(d) >= tol1 ? x+d : x+
SIGN(tol1,d));
111 if (fu <= fw || w == x) {
116 }
else if (fu <= fv || v == x || v == w) {
122 fprintf(stderr,
"too many iterations in brent()\n");
127 template<
class T,
class U>
128 double dbrent(
double ax,
double bx,
double cx, T* obj, U* dobj,
double tol,
double& xmin) {
130 double a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
131 double fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
133 a=(ax < cx ? ax : cx);
134 b=(ax > cx ? ax : cx);
137 dw=dv=dx=(*dobj)(&x);
138 for (iter=1; iter<=
ITMAX; iter++) {
140 tol1=tol*fabs(x)+
ZEPS;
142 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
146 if (fabs(e) > tol1) {
149 if (dw != dx) d1=(w-x)*dx/(dx-dw);
150 if (dv != dx) d2=(v-x)*dx/(dx-dv);
153 ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
154 ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
159 d=(fabs(d1) < fabs(d2) ? d1 : d2);
164 if (fabs(d) <= fabs(0.5*olde)) {
166 if (u-a < tol2 || b-u < tol2)
169 d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
172 d=0.5*(e=(dx >= 0.0 ? a-x : b-x ));
175 d=0.5*(e=(dx >= 0.0 ? a-x : b-x ));
177 if (fabs(d) >= tol1) {
192 MOV3(v,fv,dv, w,fw,dw);
193 MOV3(w,fw,dw, x,fx,dx);
194 MOV3(x,fx,dx, u,fu,du);
198 if (fu <= fw || w == x) {
199 MOV3(v,fv,dv, w,fw,dw);
200 MOV3(w,fw,dw, u,fu,du);
201 }
else if (fu < fv || v == x || v == w) {
202 MOV3(v,fv,dv, u,fu,du);
206 fprintf(stderr,
"too many iterations in dbrent\n");
211 void lnsrch(
int n,
double xold[],
double fold,
double g[],
double p[],
double x[],
212 double* f,
double stpmax,
int* check, T* obj) {
213 const double ALF=1.0e-4;
214 const double TOLX=1.0e-7;
216 static double maxarg1,maxarg2;
217 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
220 double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,
224 for (sum=0.0,i=0; i<n; i++) sum += p[i]*p[i];
227 for (i=0; i<n; i++) p[i] *= stpmax/sum;
228 for (slope=0.0,i=0; i<n; i++)
231 for (i=0; i<n; i++) {
232 temp=fabs(p[i])/
FMAX(fabs(xold[i]),1.0);
233 if (temp > test) test=temp;
238 for (i=0; i<n; i++) x[i]=xold[i]+alam*p[i];
241 for (i=0; i<n; i++) x[i]=xold[i];
244 }
else if (*f <= fold+ALF*alam*slope)
return;
247 tmplam = -slope/(2.0*(*f-fold-slope));
249 rhs1 = *f-fold-alam*slope;
250 rhs2=f2-fold2-alam2*slope;
251 a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
252 b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
253 if (a == 0.0) tmplam = -slope/(2.0*b);
255 disc=b*b-3.0*a*slope;
256 if (disc<0.0) fprintf(stderr,
"Roundoff problem in lnsrch.\n");
257 else tmplam=(-b+sqrt(disc))/(3.0*a);
266 alam=
FMAX(tmplam,0.1*alam);
270 template<
class T,
class U>
271 void dfpmin(
double p[],
int n,
double gtol,
int* iter,
double* fret, T* obj, U* dobj) {
272 bool verbose =
false;
274 const double EPS=3.0e-8;
275 const double TOLX=4*EPS;
276 const double STPMAX=100.0;
278 static double maxarg1,maxarg2;
279 #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))
282 double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
283 double* dg, *g, *hdg, **hessin, *pnew, *xi;
288 hessin =
new double*[n];
290 hessin[i] =
new double[n];
291 pnew =
new double[n];
298 for (i=0; i<n; i++) {
299 for (j=0; j<n; j++) hessin[i][j]=0.0;
304 stpmax=STPMAX*
FMAX(sqrt(sum),(
double)n);
306 for (its=1; its<=
ITMAX; its++) {
308 cout <<
"dfpmin: iteration " << its <<
", p = ";
309 for (i=0; i<n; i++)cout << p[i] <<
" ";
311 for (i=0; i<n; i++)cout << xi[i] <<
" ";
315 lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,obj);
317 for (i=0; i<n; i++) {
322 for (i=0; i<n; i++) {
323 temp=fabs(xi[i])/
FMAX(fabs(p[i]),1.0);
324 if (temp > test) test=temp;
326 if (verbose) cout <<
"testtol = " << test;
328 if (verbose) cout << endl;
333 for (i=0; i<n; ++i) {
341 for (i=0; i<n; i++) dg[i]=g[i];
345 for (i=0; i<n; i++) {
346 temp=fabs(g[i])*
FMAX(fabs(p[i]),1.0)/den;
347 if (temp > test) test=temp;
350 if (verbose) cout <<
" testgtol = " << test << endl;
356 for (i=0; i<n; ++i) {
364 for (i=0; i<n; i++) dg[i]=g[i]-dg[i];
365 for (i=0; i<n; i++) {
367 for (j=0; j<n; j++) hdg[i] += hessin[i][j]*dg[j];
369 fac=fae=sumdg=sumxi=0.0;
370 for (i=0; i<n; i++) {
373 sumdg += (dg[i]*dg[i]);
374 sumxi += (xi[i]*xi[i]);
376 if (fac > sqrt(EPS*sumdg*sumxi)) {
379 for (i=0; i<n; i++) dg[i]=fac*xi[i]-fad*hdg[i];
380 for (i=0; i<n; i++) {
381 for (j=i; j<n; j++) {
382 hessin[i][j] += fac*xi[i]*xi[j]
383 -fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
384 hessin[j][i]=hessin[i][j];
388 for (i=0; i<n; i++) {
390 for (j=0; j<n; j++) xi[i] -= hessin[i][j]*g[j];
393 fprintf(stderr,
"too many iterations in dfpmin");
397 for (i=0; i<n; ++i) {
#define SHFT3(a, b, c, d)
Definition: algo.h:27
void dfpmin(double p[], int n, double gtol, int *iter, double *fret, T *obj, U *dobj)
Definition: algo.h:271
const double CGOLD
Definition: algo.h:20
#define SIGN(a, b)
Definition: algo.h:28
#define MOV3(a, b, c, d, e, f)
Definition: algo.h:29
const double GOLD
Definition: algo.h:17
const double R
Definition: algo.h:15
const double C
Definition: algo.h:16
const double ZEPS
Definition: algo.h:22
double dbrent(double ax, double bx, double cx, T *obj, U *dobj, double tol, double &xmin)
Definition: algo.h:128
double brent(double ax, double bx, double cx, T *obj, double tol, double &xmin)
Definition: algo.h:65
#define SHFT2(a, b, c)
Definition: algo.h:26
const int ITMAX
Definition: algo.h:23
void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[], double *f, double stpmax, int *check, T *obj)
Definition: algo.h:211
double golden(double ax, double bx, double cx, T *obj, double tol, double &xmin)
Definition: algo.h:32
const double TINY
Definition: algo.h:19
const double GLIMIT
Definition: algo.h:18
#define SHFT(a, b, c, d)
Definition: algo.h:25