module FminFMM

!!*****************************************************************************
!! Procedury Fmin, BrentMin z knihy
!! "Computer Methods for Mathematical  Computations",
!! by George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler
!! Prentice-Hall, 1977.
!!
!! Prepis pro Fortran 90:
!! Last updated: 28 February 2005 by Ralph Carmichael ralph@pdas.com 
!! Public Domain Aeronautical Software 
!! P.O. Box 1438 Santa Cruz CA 95061 USA 
!!
!! Drobne upravy pro g95
!! Jan Cely, UFKL,Prirodovedecka fakulta MU, leden 2007,brezen 2008
!!*****************************************************************************
!!
!!	Vsechny procedury pracuji v jednoduche (SP) i dvojite (DP) presnosti
!!.............................................................................
!!
!!	SUBROUTINE BrentMin(ax,bx,F,tol,maxIter,neval,errCode,xMin,fMin)
!!
!!	Brentova procedura pro nalezeni minima funkce F(x) na intervalu <ax,bx>
!!	VSTUP
!!	ax, bx		levy a pravy okraj intervalu pro lokalizaci minima
!!	F			funkce F(x)
!!	tol			pozadovany interval pro lokalizaci minima
!!	maxIter		maximalni povoleny pocet iteraci
!!	VYSTUP
!!	neval		INTEGER, pocet provedenych vypoctu funkce F(x)
!!	errCOde		INTEGER, 0 = OK, 1 = dosazen max. pocet iteraci
!!	xMin		poloha minima
!!	fMin		F(xMin)
!!
!!*****************************************************************************
!!
!!	FUNCTION Fmin(ax,bx,f,tol)
!!
!!	Hleda minimum funkce F(x) na intervalu <ax,bx>
!!	Vola proceduru BrentMin s s hodnotou MaxIter = 100
!!	VSTUP
!!	ax, bx		levy a pravy okraj intervalu pro lokalizaci minima
!!	F			funkce F(x)
!!	tol			pozadovany interval pro lokalizaci minima
!!*****************************************************************************


implicit none

private

INTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(4)
INTEGER,PARAMETER:: DP = SELECTED_REAL_KIND(8)

character(len=20),parameter,public :: FminFMMVer = "1.0" ! 7.4.2008

PUBLIC:: Fmin,BrentMin
private :: FminSingle,BrentMinSingle,FminDouble,BrentMinDouble

INTERFACE Fmin
	MODULE PROCEDURE FminSingle,FminDouble
END INTERFACE
  
INTERFACE BrentMin
    MODULE PROCEDURE BrentMinSingle,BrentMinDouble
END INTERFACE


!!****************************************************************************
!!
contains

!!*****************************************************************************
!!							Single
!!*****************************************************************************
!!
FUNCTION FminSingle (ax,bx,f,tol) RESULT (xopt)
! ---------------------------------------------------------------------------
! PURPOSE -

  INTERFACE
    FUNCTION F(x) RESULT(ff)  ! integrand function subprogram
    IMPLICIT NONE
      	iNTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(4)
      REAL(SP),INTENT(IN)::x
      REAL(SP):: ff
    END FUNCTION F
  END INTERFACE

  REAL(SP),INTENT(IN):: ax   !  left endpoint of initial interval
  REAL(SP),INTENT(IN):: bx   !  right endpoint of initial interval
  REAL(SP),INTENT(IN):: tol  !  desired length of interval of uncertainity
  REAL(SP) :: xopt ! the result
  
  INTEGER:: errCode
  REAL(SP):: fMin
  INTEGER,PARAMETER:: MAX_ITER = 100
  INTEGER:: neval
!----------------------------------------------------------------------------
  CALL BrentMinSingle(ax,bx,F,tol,MAX_ITER,neval,errCode,xopt,fMin)
  
  RETURN
END Function FminSingle   
!------------------------------------------------------------------------------

!+
SUBROUTINE BrentMinSingle(ax,bx,F,tol,maxIter,neval,errCode,xMin,fMin)
! ---------------------------------------------------------------------------
! PURPOSE -

  INTERFACE
    FUNCTION F(x) RESULT(ff)  ! integrand function subprogram
    IMPLICIT NONE
      INTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(6,20)
      REAL(SP),INTENT(IN)::x
      REAL(SP):: ff
    END FUNCTION F
  END INTERFACE

  REAL(SP),INTENT(IN):: ax   	!  left endpoint of initial interval
  REAL(SP),INTENT(IN):: bx  	!  right endpoint of initial interval
  REAL(SP),INTENT(IN):: tol  	!  desired length of interval of uncertainity
  INTEGER,INTENT(IN):: maxIter  ! max number of iterations allowed. 25 is good
  INTEGER,INTENT(OUT):: neval
  INTEGER,INTENT(OUT):: errCode 	! =0 is OK; =1 too many iterations
  REAL(SP),INTENT(OUT):: xMin,fMin 	! the last and best value of the minimum                                  
      
  REAL(SP), PARAMETER :: C= (3.0-SQRT(5.0))/2 ! zlaty rez

  REAL(SP) :: a,b,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
  REAL(SP) :: fu,fv,fw,fx,x
  INTEGER :: iter
  REAL(SP),PARAMETER:: ZERO=0.0, ONE=1.0, TWO=2.0, THREE=3.0, HALF=0.5
!----------------------------------------------------------------------------
  errCode=0
  eps=SQRT(EPSILON(ONE))

  a=ax ! initialization
  b=bx
  v=a+c*(b-a)
  w=v
  x=v
  e=ZERO
  fx=F(x)
  fv=F(v)
  neval=2
  fw=fx

  DO iter=1,maxIter   							! iterate until minimum is found
    xm=HALF*(a+b)
    tol1=eps*ABS(x)+tol/THREE
    tol2=tol1+tol1
    IF (ABS(x-xm)<=(tol2-HALF*(b-a)) ) EXIT    	! check stopping criterion.
												! -------------------------
    IF (ABS(e) > tol1) THEN   ! is golden section necessary???
      r=(x-w)*(fx-fv)   ! trial parabolic fit
      q=(x-v)*(fx-fw)
      p=(x-v)*q-(x-w)*r
      q=TWO*(q-r)
      IF (q > ZERO) p=-p
      q=ABS(q)
      r=e
      e=d
      IF (ABS(p)>=ABS(HALF*q*r) .OR. p<=q*(a-x) .OR. p>=q*(b-x)) GO TO 40

      d=p/q   ! a parabolic interpolation step
      u=x+d
      IF ((u-a)<tol2 .OR. (b-u)<tol2) D=SIGN(tol1,xm-x)
      GO TO 50
    END IF

!40       CONTINUE
 40 Continue
 	IF (x >= xm) THEN   ! a golden section step
      e=a-x
    ELSE
      e=b-x
    END IF
    d=c*e

50  Continue
	IF (ABS(d) >= tol1) THEN   ! f must not be evaluated too close to x
      u=x+d
    ELSE
      u=x+SIGN(tol1,d)
    END IF
    fu=F(u)
    neval=neval+1

    IF (fu <= fx) THEN   ! update a,b,v,w, and x
      IF (u >= x) THEN
        a=x
      ELSE
        b=x
      END IF
      v=w
      fv=fw
      w=x
      fw=fx
      x=u
      fx=fu
    ELSE
      IF (u<x) THEN
        a=u
      ELSE
        b=u
      END IF
      IF (fu<=fw .OR. w==x) THEN
        v=w
        fv=fw
        w=u
        fw=fu
      ELSE IF (fu<=fv .OR. v==x .OR. v==w) THEN
        v=u
        fv=fu
      END IF
    END IF
  END DO

   xMin=x   ! end of main loop
   fMin=fx
  RETURN
END Subroutine BrentMinSingle   

!!============================================================================
!!							Double
!!============================================================================

FUNCTION FminDouble (ax,bx,f,tol) RESULT (xopt)
! ---------------------------------------------------------------------------
! PURPOSE -

  INTERFACE
    FUNCTION F(x) RESULT(ff)  ! integrand function subprogram
    IMPLICIT NONE
      INTEGER,PARAMETER:: DP = SELECTED_REAL_KIND(10,50)
      REAL(DP),INTENT(IN)::x
      REAL(DP):: ff
    END FUNCTION F
  END INTERFACE

  REAL(DP),INTENT(IN):: ax   !  left endpoint of initial interval
  REAL(DP),INTENT(IN):: bx   !  right endpoint of initial interval
  REAL(DP),INTENT(IN):: tol  !  desired length of interval of uncertainity
  REAL(DP) :: xopt ! the result
  
  INTEGER:: errCode
  REAL(DP):: fMin
  INTEGER,PARAMETER:: MAX_ITER = 100
  INTEGER:: neval
!----------------------------------------------------------------------------
  CALL BrentMinDouble(ax,bx,F,tol,MAX_ITER,neval,errCode,xopt,fMin)
  
  RETURN
END Function FminDouble   
!******************************************************************************

!+
SUBROUTINE BrentMinDouble(ax,bx,F,tol,maxIter,neval,errCode,xMin,fMin)
! ---------------------------------------------------------------------------
! PURPOSE -

  INTERFACE
    FUNCTION F(x) RESULT(ff)  ! integrand function subprogram
    IMPLICIT NONE
      INTEGER,PARAMETER:: DP = SELECTED_REAL_KIND(8)
      REAL(DP),INTENT(IN)::x
      REAL(DP):: ff
    END FUNCTION F
  END INTERFACE

  REAL(DP),INTENT(IN):: ax   !  left endpoint of initial interval
  REAL(DP),INTENT(IN):: bx   !  right endpoint of initial interval
  REAL(DP),INTENT(IN):: tol  !  desired length of interval of uncertainity
  INTEGER,INTENT(IN):: maxIter   ! max number of iterations allowed. 25 is good
  INTEGER,INTENT(OUT):: neval
  INTEGER,INTENT(OUT):: errCode   ! =0 is OK; =1 too many iterations
  REAL(DP),INTENT(OUT):: xMin,fMin ! the last and best value of the minimum                                  
      
  REAL(DP), PARAMETER :: C= (3.0-SQRT(5.0))/2	! zlaty rez

  REAL(DP) :: a,b,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
  REAL(DP) :: fu,fv,fw,fx,x
  INTEGER :: iter
  REAL(DP),PARAMETER:: ZERO=0.0, ONE=1.0, TWO=2.0, THREE=3.0, HALF=0.5
!----------------------------------------------------------------------------
  errCode=0
  eps=SQRT(EPSILON(ONE))

  a=ax ! initialization
  b=bx
  v=a+c*(b-a)
  w=v
  x=v
  e=ZERO
  fx=F(x)
  fv=F(v)
  neval=2
  fw=fx

DO iter=1,maxIter   				! iterate until minimum is found
    xm=HALF*(a+b)
    tol1=eps*ABS(x)+tol/THREE
    tol2=tol1+tol1
    IF (ABS(x-xm)<=(tol2-HALF*(b-a)) ) EXIT    	! check stopping criterion.
												! -------------------------
    IF (ABS(e) > tol1) THEN   		! is golden section necessary???
      r=(x-w)*(fx-fv)   			! trial parabolic fit
      q=(x-v)*(fx-fw)
      p=(x-v)*q-(x-w)*r
      q=TWO*(q-r)
      IF (q > ZERO) p=-p
      q=ABS(q)
      r=e
      e=d
      IF (ABS(p)>=ABS(HALF*q*r) .OR. p<=q*(a-x) .OR. p>=q*(b-x)) GO TO 40

      d=p/q   ! a parabolic interpolation step
      u=x+d
      IF ((u-a)<tol2 .OR. (b-u)<tol2) D=SIGN(tol1,xm-x)
      GO TO 50
    END IF

!40       CONTINUE
 40 Continue
 	IF (x >= xm) THEN   ! a golden section step
      e=a-x
    ELSE
      e=b-x
    END IF
    d=c*e

50  Continue
	IF (ABS(d) >= tol1) THEN   ! f must not be evaluated too close to x
      u=x+d
    ELSE
      u=x+SIGN(tol1,d)
    END IF
    fu=F(u)
    neval=neval+1

    IF (fu <= fx) THEN   ! update a,b,v,w, and x
      IF (u >= x) THEN
        a=x
      ELSE
        b=x
      END IF
      v=w
      fv=fw
      w=x
      fw=fx
      x=u
      fx=fu
    ELSE
	  IF (u<x) THEN
        a=u
      ELSE
        b=u
      END IF
      IF (fu<=fw .OR. w==x) THEN
        v=w
        fv=fw
        w=u
        fw=fu
      ELSE IF (fu<=fv .OR. v==x .OR. v==w) THEN
        v=u
        fv=fu
      END IF
    END IF
END DO

   xMin=x   ! end of main loop
   fMin=fx
  RETURN
END Subroutine BrentMinDouble 
!!*****************************************************************************


end module FminFMM
