!********************************************************
!*        Find minimum of a real function Y=F(X)        *
!* ---------------------------------------------------- *
!* REFERENCE: "Numerical Recipes, The Art of Scientific *
!*             Computing By W.H. Press, B.P. Flannery,  *
!*             S.A. Teukolsky and W.T. Vetterling,      *
!*             Cambridge University Press, 1986"        *
!*             [BIBLI 08].                              *
!* ---------------------------------------------------- *
!* Sample run:                                          *
!*                                                      *
!* Find a minimum of function X*SIN(X)-2*COS(X)         *
!* between X=4 and X=6.                                 *
!*                                                      *
!* X1=4.d0  X2=5.d0  X3=6.d0  TOL=1.d-6                 *
!*                                                      *
!* Function minimum is:  -5.53452877399827              *
!*                                                      *
!* for X=   5.23293671913279                            *
!*                                                      *
!********************************************************
PROGRAM SEEK_MINIMUM

REAL*8 X1,X2,X3,XMINI,YMINI,TOL

X1=4.d0; X2=5.d0; X3=6.d0
TOL=1.d-6

YMINI=BRENT(X1,X2,X3,TOL,XMINI)

print *,''
print *,' Function minimum is ',YMINI
print *,''
print *,' for X=',XMINI
print *,''

stop
END


!Function to be analyzed
REAL*8 FUNCTION F(X)
  REAL*8 X  
  F=X*SIN(X)-2.*COS(X)
RETURN
END  


REAL*8 FUNCTION BRENT(AX,BX,CX,TOL,XMIN)
!Given a function F, and a bracketing triplet of abscissas
!AX,BX,CX (such that BX is between AX and CX, and F(BX) is less 
!than both F(AX) and F(CX)), this routine isolates the minimum 
!to a fractional precision of about TOL using Brent's method.
!The abscissa of the minimum is returned in XMIN, and the minimum
!function value is returned as BRENT, the returned function value.
PARAMETER(ITMAX=100,CGOLD=.3819660,ZEPS=1.0E-10)
!Maximum allowed number of iterations; golden ration; and a small
!number which protects against trying to achieve fractional accuracy
!for a minimum that happens to be exactly zero.
IMPLICIT REAL*8 A-H,O-Z

A=MIN(AX,CX)
B=MAX(AX,CX)
V=BX
W=V
X=V
D=0.
E=0.
FX=F(X)
FV=FX
FW=FX
DO 11 ITER=1,ITMAX	                                !main loop
  XM=0.5*(A+B)
  TOL1=TOL*ABS(X)+ZEPS
  TOL2=2.*TOL1
  IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3  !Test for done here
  IF(ABS(E).GT.TOL1) THEN     !Construct a trial parabolic fit
    R=(X-W)*(FX-FV)
	Q=(X-V)*(FX-FW)
    P=(X-V)*Q-(X-W)*R
	Q=2.D0*(Q-R)         !bug corrected 07/24/2006 (0.2 instead of 2.D0)
	IF(Q.GT.0)  P=-P
	Q=ABS(Q)
	ETEMP=E
    E=D
	IF (ABS(P).GE.ABS(.5*Q*ETEMP).OR.P.LE.Q*(A-X).OR.  &
	P.GE.Q*(B-X))  GOTO 1
!   The above conditions determine the acceptability of the 
!   parabolic fit. Here it is o.k.:
    D=P/Q
	U=X+D
	IF(U-A.LT.TOL2.OR.B-U.LT.TOL2)  D=SIGN(TOL1,XM-X)
	GOTO 2
  ENDIF
1 IF(X.GE.XM) THEN
    E=A-X
  ELSE
    E=B-X
  ENDIF
  D=CGOLD*E
2 IF(ABS(D).GE.TOL1) THEN
    U=X+D
  ELSE
    U=X+SIGN(TOL1,D)
  ENDIF
  FU=F(U)  !This the one function evaluation per iteration
  IF(FU.LE.FX) THEN
    IF(U.GE.X) THEN
	  A=X
    ELSE
	  B=X
	ENDIF
	V=W
	FV=FW
	W=X
	FW=FX
	X=U
	FX=FU
  ELSE
    IF(U.LT.X) THEN
	  A=U
	ELSE
	  B=U
	ENDIF
	IF(FU.LE.FW.OR.W.EQ.X) THEN
	  V=W
	  FV=FW
	  W=U
	  FW=FU
    ELSE IF(FU.LE.FV.OR.V.EQ.X.OR.V.EQ.W) THEN
	  V=U
	  FV=FU
    ENDIF
  ENDIF
11 CONTINUE
Pause ' Brent exceed maximum iterations.'
3 XMIN=X   !exit section
  BRENT=FX
  RETURN
  END
  
!end of file Brent.f90