module base_solver_mod
  use  base_chains_mod
  implicit none

  real, dimension(:), allocatable, private :: dcal
  real, dimension(:), allocatable, private :: dtmp
contains
  !!---------------------------------------------------------------------------
  !!.. Conjugate gradient solver for linear positive definite system
  !!.. W(L'Lm - L'd) = 0
  !!.. Fop = L'L, d = L'd
  subroutine base_solver_cglpd_smp(m, d, Fop, niter, search, Wop, beta_max, m0, mmov, rmov, fmov)
    optional             :: Wop, beta_max, m0, mmov, rmov, fmov
    real, dimension(:)   :: m
    real, dimension(:)   :: d
    integer              :: niter
    logical              :: search
    real                 :: beta_max
    real, dimension(:)   :: m0    
    real, dimension(:,:) :: mmov                       !!.. model movie
    real, dimension(:,:) :: rmov                       !!.. residual movie
    real, dimension(:)   :: fmov                       !!.. objective function value movie
    interface
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface

    real, dimension(:), allocatable :: s
    real, dimension(:), allocatable :: r
    real, dimension(:), allocatable :: t
    real, dimension(:), allocatable :: b
    double precision                :: rdro, rdrp, sdso
    double precision                :: alph, beta 
    double precision                :: ax, bx, cx
    double precision                :: fa, fb, fc
    double precision                :: tol, fmin, xmin
    integer                         :: iter, stat, i

    allocate(s(size(m)))
    allocate(r(size(d)))
    allocate(t(size(d))) 
    allocate(b(size(d)))

    tol = 1.d-8

    m =  0.
    r = -d

    !!.. weight the residual r = W r
    if (present(Wop)) stat = Wop(.false., .false., -d, r)

    s = -r !!.. steepest descent

    alph = 0.
    beta = 0.

    do iter=1, niter

       !!.. update the search direction
       s    = -r + beta * s
       !!.. b = F s, data space
       stat = Fop(.false., .false., s, b)
       !!.. b = W F s, apply weighting      
       if (present(Wop)) then
          t = b
          stat = Wop(.false., .false., t, b)
       end if
       !!.. dot-product: (Wr)'Wr
       rdro = sum(dprod(r,r))
       !!.. dot-product: s'WFs                
       sdso = sum(dprod(s,b))
       !!.. step length based on quadratic approximation
       alph = rdro / max(sdso, 1.d-12)

       if (search) then
          !!.. initial brackets for line search
          ax   = 0.
          bx   = alph * 1.618
          call base_mnbrak_linfnc(r, b, ax, bx, cx, fa, fb, fc)
          !!.. now we do a line search: Golden section search
          call base_golden_linfnc(r, b, ax, bx, cx, tol, fmin, xmin)
          !!.. get the robust step length
          alph = xmin
       end if

       !!.. update the model
       m    = m + alph * s
       !!.. update the residual
       r    = r + alph * b
       !!.. dot-product: r'r
       rdrp = sum(dprod(r, r))
       !!.. step length for search direction
       !!.. beta is to be determined by the requirement that s_{k-1} and s_k
       !!.. must be conjugate with respect to F
       beta = rdrp / max(rdro, 1.d-12)

       if (present(beta_max)) then
          if (beta > beta_max) beta = beta*beta_max
       end if
    
       write(0,*)"--------------------report------------------------------------"
       write(0,*)iter, rdrp, beta

       if (present(mmov)) mmov(:,iter) = m
       if (present(rmov)) rmov(:,iter) = r
       if (present(fmov)) fmov(iter)   = real(rdrp)
       if (abs(rdro-rdrp)<=tol) exit
    end do  

    deallocate(s)
    deallocate(r)
    deallocate(t) 
    deallocate(b)
    return 
  end subroutine

  !!-----------------------------------------------------------------------------------------------
  !!.. with regularization
  !!..
  !!..  L      d
  !!..     m =    => W((L'L+A'A)m - L'd) = 0
  !!..  A      0
  !!..  Fop = L'L , Aop = A, d = L'd 
  subroutine base_solver_cglpd_reg(m, d, Fop, Aop, niter, eps, search, Wop, beta_max, m0, mmov, rmov, fmov)
    optional             :: Wop, beta_max, m0, mmov, rmov, fmov
    real, dimension(:)   :: m
    real, dimension(:)   :: d
    integer              :: niter
    real, dimension(:)   :: eps
    logical              :: search
    real                 :: beta_max
    real, dimension(:)   :: m0 
    real, dimension(:,:) :: mmov                       !!.. model movie
    real, dimension(:,:) :: rmov                       !!.. data residual movie
    real, dimension(:)   :: fmov                       !!.. objective function value movie
    interface
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Aop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface

    real, dimension(:), allocatable :: s
    real, dimension(:), allocatable :: r
    real, dimension(:), allocatable :: t1, t2
    real, dimension(:), allocatable :: b
    double precision                :: rdro, rdrp, sdso
    double precision                :: alph, beta
    double precision                :: ax, bx, cx
    double precision                :: fa, fb, fc
    double precision                :: tol, fmin, xmin
    integer                         :: iter, stat, i

    allocate(s( size(m)))
    allocate(r( size(m)))    !!.. residual
    allocate(t1(size(m)))
    allocate(t2(size(m))) 
    allocate( b(size(m)))

    tol = 1.d-8

    m = 0.
    r = -d

    !!.. r = W r
    if (present(Wop)) then
       stat = Wop(.false., .false., -d, r)
    end if

    s = -r !!.. steepest descent

    alph = 0.
    beta = 0.

    do iter=1, niter

       !!.. update the search direction
       s    = -r + beta * s

       !!.. b = F s, data space
       stat = Fop(.false., .false.,  s,  b)

       !!.. b = W F s, apply weighting      
!       if (present(Wop)) then
!          t1 = b
!          stat = Wop(.false., .false., t1, b)
!       end if

       !!.. t1 = A s
       stat = Aop(.false., .false.,  s, t1)
       !!.. t2 = A'A s
       stat = Aop( .true.,  .true., t2, t1)

       !!.. b = F s + A'A s
       b = b + eps*t2

       !!.. b = W F s, apply weighting      
       if (present(Wop)) then
          t1 = b
          stat = Wop(.false., .false., t1, b)
       end if

       !!.. dot-product: (Wr)'Wr
       rdro = sum(dprod(r,r))
       !!.. dot-product: s'WFs                
       sdso = sum(dprod(s,b))
       !!.. step length based on quadratic approximation
       alph = rdro / max(sdso, 1.d-12)

       if (search) then
          !!.. initial brackets for line search
          ax   = 0.
          bx   = alph * 1.618
          call base_mnbrak_linfnc(r, b, ax, bx, cx, fa, fb, fc)
          !!.. now we do a line search: Golden section search
          call base_golden_linfnc(r, b, ax, bx, cx, tol, fmin, xmin)
          !!.. get the robust step length
          alph = xmin
       end if

       !!.. update the model
       m    = m + alph * s
       !!.. update the residual
       r    = r + alph * b
       !!.. dot-product: r'r
       rdrp = sum(dprod(r, r))
       !!.. step length for search direction
       !!.. beta is to be determined by the requirement that s_{k-1} and s_k
       !!.. must be conjugate with respect to F
       beta = rdrp / max(rdro, 1.d-12)

       if (present(beta_max)) then
          if (beta > beta_max) beta = beta*beta_max
       end if

       write(0,*)"--------------------report------------------------------------"
       write(0,*)iter, rdrp, beta

       if (present(mmov)) mmov(:,iter) = m
       if (present(rmov)) rmov(:,iter) = r
       if (present(fmov)) fmov(iter)   = real(rdrp)
       if (abs(rdro-rdrp)<=tol) exit
    end do

    deallocate( s)
    deallocate( r)
    deallocate(t1)
    deallocate(t2)
    deallocate( b)
    return
  end subroutine





  !!------------------------------------------------------------------------------------
  !!.. Preconditioned conjugate gradient solver for linear positive definite system
  !!.. W(Lm - d) = 0
  !!------------------------------------------------------------------------------------
  subroutine base_solver_cglpd_prc(m, d, Fop, Sop, niter, search, Wop, beta_max, mmov, rmov, fmov)
    optional             :: Wop, beta_max, mmov, rmov, fmov
    real, dimension(:)   :: m
    real, dimension(:)   :: d
    integer              :: niter
    logical              :: search
    real                 :: beta_max
    real, dimension(:,:) :: mmov
    real, dimension(:,:) :: rmov
    real, dimension(:)   :: fmov

    interface
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Sop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface

    real, dimension(:), allocatable :: s
    real, dimension(:), allocatable :: r
    real, dimension(:), allocatable :: t
    real, dimension(:), allocatable :: b
    real, dimension(:), allocatable :: p
    double precision                :: rdro, rdrp, sdso
    double precision                :: alph, beta
    double precision                :: ax, bx, cx
    double precision                :: fa, fb, fc
    double precision                :: tol, xmin, fmin 
    integer                         :: iter, stat, i

    allocate(s(size(m)))
    allocate(r(size(d)))
    allocate(t(size(d)))
    allocate(b(size(d)))
    allocate(p(size(m)))

    tol = 1.d-8

    b =  d
    m =  0.
    r = -d

    !!.. weight the residual r = W r
    if (present(Wop)) stat = Wop(.false., .false., -d, r)

    s = -r !!.. steepest descent


    alph = 0.d-15
    beta = 0.d-15

    do iter=1, niter
       stat = Sop(.true. , .false., p, r)            !!.. p = S' r
       s    = - p + beta * s
       stat = Fop(.false., .false., s, b)            !!.. b = F s, data space
       if (present(Wop)) then
          t = b
          b = 0.
          stat = Wop(.false., .false., t, b)
       end if
       rdro = sum(dprod(r,p))
       sdso = sum(dprod(s,b))
       alph = rdro / max(sdso, 1.d-12)

       if (search) then
          !!.. initial brackets for line search
          ax   =  0.
          bx   =  alph*1.618
          write(0,*)"------------------initial-------------------------------------"
          write(0,*)"ax,bx",ax,bx
          !call base_mnbrak_genfnc(base_objfnc, Fop, m, d, s, ax, bx, cx, fa, fb, fc, Wop)
          call base_mnbrak_linfnc(r, b, ax, bx, cx, fa, fb, fc)
          write(0,*)"------------------starting------------------------------------"
          write(0,*)"x:",ax,bx,cx
          write(0,*)"f:",fa,fb,fc

          !!.. now we do a line search: Golden section search
          !call base_golden_genfnc(base_objfnc, Fop, m, d, s, ax, bx, cx, tol, fmin, xmin, Wop)
          call base_golden_linfnc(r, b, ax, bx, cx, tol, fmin, xmin)
          write(0,*)"-------------------updates------------------------------------"
          write(0,*)"alph=", alph
          write(0,*)"xmin=", xmin
          write(0,*)"fmin=", fmin

          !!.. get the robust step length
          alph = xmin
       end if

       m    = m + alph * s
       r    = r + alph * b
       stat = Sop(.true. , .false., p, r)
       rdrp = sum(dprod(r, p))
       beta = rdrp / max(rdro, 1.d-12)

       if (present(beta_max)) then
          if (beta > beta_max) beta = beta*beta_max
       end if
       write(0,*)"--------------------report------------------------------------"
       write(0,*)iter, fmin

       if (present(mmov)) mmov(:,iter) = m
       if (present(rmov)) rmov(:,iter) = r
       if (present(fmov)) fmov(iter)   = fmin
    end do

    deallocate(s)
    deallocate(r)
    deallocate(t)
    deallocate(b) 
    deallocate(p)
    return
  end subroutine
  !----------------------------------------------------------------------------
  !.. objective function
  !----------------------------------------------------------------------------
  double precision function base_objfnc(Fop, m, d, Wop)
    optional           :: Wop
    real, dimension(:) :: m
    real, dimension(:) :: d
    interface
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface
    integer            :: stat
    if (.not.allocated(dcal)) allocate(dcal(size(d)))
    !!.. dcal = F m
    stat = Fop(.false., .false., m, dcal)
    !!.. dcal = F m - d
    dcal = dcal - d
    !!.. dcal = W (F m - d)
    if (present(Wop)) then
       if (.not.allocated(dtmp)) allocate(dtmp(size(d)))
       dtmp = dcal
       dcal = 0.
       stat = Wop(.false., .false., dtmp, dcal)
    end if
    base_objfnc = sum(dprod(dcal, dcal))
  end function
  !!-----------------------------------------------------------------------------
  !!.. Golden section search method
  !!-----------------------------------------------------------------------------
  subroutine base_golden_genfnc(func, Fop, m, d, s, ax, bx, cx, tol, fmin, xmin, Wop)
    optional :: Wop
    interface
      double precision function func(Fop, m, d, Wop)
        optional :: Wop
        interface 
          integer function Fop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
          integer function Wop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
        end interface
        real, dimension(:) :: m
        real, dimension(:) :: d
      end function
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface
    real, dimension(:) :: m, d, s
    double precision   :: ax, bx, cx
    double precision   :: tol, fmin, xmin

    double precision   :: x0, x1, x2, x3
    double precision   :: f1, f2, R, C
    parameter (R=.61803399, C=1.-R)

    x0 = ax
    x3 = cx
    if (abs(cx-bx)>abs(bx-ax)) then
       x1 = bx
       x2 = bx + C*(cx-bx)
    else
       x2 = bx
       x1 = bx - C*(bx-ax)
    end if

    if (present(Wop)) then
       f1 = func(Fop, real(m+x1*s), d, Wop)
       f2 = func(Fop, real(m+x2*s), d, Wop)
    else
       f1 = func(Fop, real(m+x1*s), d)
       f2 = func(Fop, real(m+x2*s), d)
    end if

    do while(abs(x3-x0)>tol*(abs(x1)+abs(x2)))
       if (f2<f1) then
          x0 = x1
          x1 = x2
          x2 = R*x1 + C*x3
          f1 = f2
          if (present(Wop)) then
             f2 = func(Fop, real(m+x2*s), d, Wop)
          else
             f2 = func(Fop, real(m+x2*s), d)
          end if
            
       else
          x3 = x2
          x2 = x1
          x1 = R*x2 + C*x0
          f2 = f1
          if (present(Wop)) then
             f1 = func(Fop, real(m+x1*s), d, Wop)
          else
             f1 = func(Fop, real(m+x1*s), d)
          end if 
       end if
    end do
    if (f1<f2) then
       fmin = f1
       xmin = x1
    else
       fmin = f2
       xmin = x2
    end if
    return
  end subroutine

  subroutine base_golden_linfnc(rr, dr, ax, bx, cx, tol, fmin, xmin)
    real, dimension(:) :: rr, dr
    double precision   :: ax, bx, cx
    double precision   :: tol, fmin, xmin

    double precision   :: x0, x1, x2, x3
    double precision   :: f1, f2, R, C
    parameter (R=.61803399, C=1.-R)

    x0 = ax
    x3 = cx
    if (abs(cx-bx)>abs(bx-ax)) then
       x1 = bx
       x2 = bx + C*(cx-bx)
    else
       x2 = bx
       x1 = bx - C*(bx-ax)
    end if

    f1 = sum(dprod(rr+x1*dr, rr+x1*dr))
    f2 = sum(dprod(rr+x2*dr, rr+x2*dr))

    do while(abs(x3-x0)>tol*(abs(x1)+abs(x2)))
       if (f2<f1) then
          x0 = x1
          x1 = x2
          x2 = R*x1 + C*x3
          f1 = f2
          f2 = sum(dprod(rr+x2*dr, rr+x2*dr))
       else
          x3 = x2
          x2 = x1
          x1 = R*x2 + C*x0
          f2 = f1
          f1 = sum(dprod(rr+x1*dr, rr+x1*dr))
       end if
    end do
    if (f1<f2) then
       fmin = f1
       xmin = x1
    else
       fmin = f2
       xmin = x2
    end if
    return
  end subroutine
  !------------------------------------------------------------------------------
  subroutine base_brents_genfnc(func, Fop, mm, dd, ss, ax, bx, cx, tol, fmin, xmin, ITMAX, Wop)
    optional :: Wop
    interface
      double precision function func(Fop, m, d, Wop)
        optional :: Wop
        interface
          integer function Fop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
          integer function Wop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
        end interface
        real, dimension(:) :: m
        real, dimension(:) :: d
      end function
      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface
    real, dimension(:) :: mm, dd, ss
    double precision   :: ax, bx, cx
    double precision   :: tol, fmin, xmin
    integer            :: ITMAX

    integer            :: iter
    double precision   :: a,b,c,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
    double precision   :: CGOLD, ZEPS
    parameter (CGOLD=0.3819660, ZEPS=1.0d-10)
    
    a = min(ax, cx)
    b = max(ax, cx)
    v = bx
    w = v
    x = v
    e = 0.
    if (present(Wop)) then
       fx = func(Fop, real(mm+x*ss), dd, Wop)
    else   
       fx = func(Fop, real(mm+x*ss), dd)
    end if
    fv = fx
    fw = fx

    do iter=1, ITMAX
       write(0,*)"in Brent line search, iter=", iter, ITMAX
       xm = 0.5*(a+b)
       tol1 = tol*abs(x) + ZEPS
       tol2 = 2*tol1
       if (abs(x-xm)<=(tol2-0.5*(b-a))) then
          xmin = x
          fmin = fx
          return
       end if
       if (abs(e)>tol1) then
          r = (x-w)*(fx-fv)
          q = (x-v)*(fx-fw)
          p = (x-v)*q-(x-w)*r
          q = 2.*(q-r)
          if (q>0.) p = -q
          q = abs(q)
          etemp = e
          e = d
          if (abs(p)>=abs(0.5*q*etemp) .or. p<=q*(a-x) .or. p>=q*(b-x)) then
             if (x>=xm) then
                e = a-x
             else
                e = b-x
             end if
             d = CGOLD*e
          else
             d = p/q
             u = x+d
          end if 
          if (u-a<tol2 .or. b-u<tol2) d=sign(tol1,xm-x)
       end if
       if (abs(d)>=tol1) then
          u = x+d
       else
          u = x+sign(tol1,d)
       end if 
       if (present(Wop)) then
          fu = func(Fop, real(mm+u*ss), dd, Wop)
       else
          fu = func(Fop, real(mm+u*ss), dd)
       end if

       if (fu<=fx) then
          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
    fmin = fx
    return
  end subroutine 
  !------------------------------------------------------------------------------
  subroutine base_mnbrak_genfnc(func, Fop, m, d, s, ax, bx, cx, fa, fb, fc, Wop)
    optional :: Wop

    interface

      double precision function func(Fop, m, d, Wop)
        optional :: Wop
        interface
          integer function Fop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
          integer function Wop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
        end interface
        real, dimension(:) :: m
        real, dimension(:) :: d
      end function

      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function

      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function 

    end interface
    real, dimension(:) :: m, d, s
    double precision   :: ax, bx, cx
    double precision   :: fa, fb, fc
    
    double precision   :: GOLD, GLIMIT, TINY
    double precision   :: dum, fu, q, r, u, ulim
    parameter (GOLD=1.618034, GLIMIT=10., TINY=1.e-20)

    if (present(Wop)) then
       fa = func(Fop, real(m+ax*s), d, Wop)
       fb = func(Fop, real(m+bx*s), d, Wop)
    else
       fa = func(Fop, real(m+ax*s), d)
       fb = func(Fop, real(m+bx*s), d)
    end if

    if (fb>fa) then
       dum = ax
       ax  = bx
       bx  = dum
       dum = fb
       fb  = fa
       fa  = dum
    end if

    cx = bx + GOLD*(bx-ax)

    if (present(Wop)) then
       fc = func(Fop, real(m*cx*s), d, Wop)
    else
       fc = func(Fop, real(m+cx*s), d)
    end if
 
    do while (fb>=fc)
       r = (bx-ax)*(fb-fc)
       q = (bx-cx)*(fb-fa)
       u = bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r))
       ulim = bx+GLIMIT*(cx-bx)
       if ((bx-u)*(u-cx)>0.) then
          if (present(Wop)) then
             fu = func(Fop, real(m+u*s), d, Wop)
          else
             fu = func(Fop, real(m+u*s), d)
          end if
          if (fu<fc) then
             ax = bx
             fa = fb
             bx = u
             fb = fu
             return
          else if(fu>fb) then
             cx = u
             fc = fu
             return
          end if
          u  = cx+GOLD*(cx-bx)
          if (present(Wop)) then
             fu = func(Fop, real(m+u*s), d, Wop)
          else
             fu = func(Fop, real(m+u*s), d)
          end if
       else if ((cx-u)*(u-ulim)>0) then
          if (present(Wop)) then
             fu = func(Fop, real(m+u*s), d, Wop)
          else
             fu = func(Fop, real(m+u*s), d)
          end if
          if (fu<fc) then
             bx = cx
             cx = u
             u  = cx+GOLD*(cx-bx)
             fb = fc
             fc = fu
             if (present(Wop)) then
                fu = func(Fop, real(m+u*s), d, Wop)
             else
                fu = func(Fop, real(m+u*s), d)
             end if
          end if
       else if ((u-ulim)*(ulim-cx)>=0.) then
          u  = ulim
          if (present(Wop)) then
             fu = func(Fop, real(m+u*s), d, Wop)
          else
             fu = func(Fop, real(m+u*s), d)
          end if 
       else
          u  = cx+GOLD*(cx-bx)
          if (present(Wop)) then
             fu = func(Fop, real(m+u*s), d, Wop)
          else
             fu = func(Fop, real(m+u*s), d)
          end if
       end if
       ax = bx
       bx = cx
       cx = u
       fa = fb
       fb = fc
       fc = fu
    end do
    return
  end subroutine

  subroutine base_mnbrak_linfnc(rr, dr, ax, bx, cx, fa, fb, fc)
    real, dimension(:) :: rr, dr
    double precision   :: ax, bx, cx
    double precision   :: fa, fb, fc

    double precision   :: GOLD, GLIMIT, TINY
    double precision   :: dum, fu, q, r, u, ulim
    parameter (GOLD=1.618034, GLIMIT=10., TINY=1.e-20)

    fa = sum(dprod(rr+ax*dr, rr+ax*dr))
    fb = sum(dprod(rr+bx*dr, rr+bx*dr))

    if (fb>fa) then
       dum = ax
       ax  = bx
       bx  = dum
       dum = fb
       fb  = fa
       fa  = dum
    end if

    cx = bx + GOLD*(bx-ax)

    fc = sum(dprod(rr+cx*dr, rr+cx*dr)) 

    do while (fb>=fc)
       r = (bx-ax)*(fb-fc)
       q = (bx-cx)*(fb-fa)
       u = bx-((bx-cx)*q-(bx-ax)*r)/(2.*sign(max(abs(q-r),TINY),q-r))
       ulim = bx+GLIMIT*(cx-bx)
       if ((bx-u)*(u-cx)>0.) then
          fu = sum(dprod(rr+u*dr, rr+u*dr))
          if (fu<fc) then
             ax = bx
             fa = fb
             bx = u
             fb = fu
             return
          else if(fu>fb) then
             cx = u
             fc = fu
             return
          end if
          u  = cx+GOLD*(cx-bx)
          fu = sum(dprod(rr+u*dr, rr+u*dr))

       else if ((cx-u)*(u-ulim)>0) then
          fu = sum(dprod(rr+u*dr, rr+u*dr))
          if (fu<fc) then
             bx = cx
             cx = u
             u  = cx+GOLD*(cx-bx)
             fb = fc
             fc = fu
             fu = sum(dprod(rr+u*dr, rr+u*dr))
          end if
       else if ((u-ulim)*(ulim-cx)>=0.) then
          u  = ulim
          fu = sum(dprod(rr+u*dr, rr+u*dr))
       else
          u  = cx+GOLD*(cx-bx)
          fu = sum(dprod(rr+u*dr, rr+u*dr))
       end if
       ax = bx
       bx = cx
       cx = u
       fa = fb
       fb = fc
       fc = fu
    end do
    return
  end subroutine

  subroutine base_solver_anneal(func, Fop, xmin, dobs, Tred, KTmax, KTacc_max, Smax, Tmin, Omin, wt, Wop)
    optional :: Wop
    interface
      double precision function func(Fop, m, d, Wop)
        optional :: Wop
        interface
          integer function Fop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
          integer function Wop(adj, add, m, d)
            logical            :: adj, add
            real, dimension(:) :: m, d
          end function
        end interface
        real, dimension(:) :: m
        real, dimension(:) :: d
      end function

      integer function Fop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function

      integer function Wop(adj, add, m, d)
        logical            :: adj, add
        real, dimension(:) :: m, d
      end function
    end interface

    real, dimension(:)              :: xmin
    real, dimension(:)              :: dobs
    real                            :: Tred
    integer                         :: KTmax
    integer                         :: KTacc_max
    integer                         :: Smax
    real                            :: Tmin
    real                            :: Omin 
    real                            :: wt
 
    real, dimension(:), allocatable :: xnow, xnew
    integer                         :: loop_inner, loop_outer, s, ktacc, kt, niter, flagT, i
    real                            :: Enow, Enew, T, dE, Prob, Probacc
    real                            :: pt, dx, scale

    allocate(xnow(size(xmin)))
    allocate(xnew(size(xmin)))

    loop_outer = 1;
    loop_inner = 1;

    xnow  = xmin
    xnew  = xnow
    scale = func(Fop, xnow, dobs, Wop)

    Enow  = 1
    T     = 1
    s     = 0
    niter = 0

write(0,*)"scale", scale

    do while (loop_outer == 1)
       ktacc = 0
       kt    = 0
       do while (loop_inner == 1)
          do i=1, size(xnow)
             call random_number(pt)
             dx      = wt*(pt-0.5)*xnow(i)
             xnew(i) = xnow(i) + dx
             niter   = niter + 1
             kt      = kt + 1
             Enew    = real(func(Fop, xnew, dobs, Wop)) / scale
             dE      = Enew - Enow
             Prob    = exp(- dE / T)
             write(0,*)"Enew=", Enew, i, niter, dx

             if (dE <= 0.) then
                xnow = xnew
                Enow = Enew
                ktacc = ktacc + 1
                if (ktacc >= KTacc_max) then
                   loop_inner = 0
                end if

             else
                call random_number(Probacc)
                if (Prob > Probacc) then
                   xnow = xnew
                   Enow = Enew
                   ktacc = ktacc + 1
                   if (ktacc >= KTacc_max) then
                      loop_inner = 0
                   end if
                end if
             end if
             if (kt >= KTmax) then
                flagT = 0
                s = s + 1
                if (s >= Smax) then
                   loop_outer = 0
                end if
             end if
             if (Enow < Omin) then
                !Minimum tolerance reached
                loop_inner = 0
                loop_outer = 0
             end if
             if (T < Tmin) then
                !Minimum temperature reached
                loop_inner = 0
                loop_outer = 0
             end if
          end do
       end do
       if (loop_inner == 0) write(0,*)"loop_inner=0"
       if (loop_outer == 1) then
          T = T*Tred
          loop_inner = 1
       end if
    end do

    xmin = xnow

    write(0,*)"niter=", niter
    write(0,*)"Enow=", Enow

  end subroutine




end module

















