module opt_com_mod

use bmad
use bmadz_struct
use bmadz_interface
use bunchcross_mod
use zquad_lens_mod
use pretz_mod
use twiss_max_mod
use constraints_mod
use nonlin_mod
use precision_def
USE nrtype;USE nrutil, ONLY : nrerror, assert_eq, diagmult
USE nr
use opt_mod


type (indep_var_struct), save  ::  indep_var_o, sext_var_o
type (lat_struct), save, target :: lat_o
type (constraint_struct), save ::  sp_o
type (db_struct), save ::        cesr_o
type (pc_struct), save ::          pc_o
type (twiss_max_struct), save ::   twiss_max_o
type (pretz_struct), save ::       pretz_o
type (global_struct), save ::      global_o
type (zquad_struct), save ::       quad_o
type (nonlin_ele_struct), save ::  nonlin_o
type (moment_struct), save ::      moment_o
type (coord_struct), allocatable, save :: co_o(:)

contains

subroutine bmad_opt (lat, indep_var, cesr, con, pc, quad,  &
     co_a, global, twiss_max, pretz, nonlin, moment, sext_var)


  implicit none

  type (indep_var_struct)  indep_var, sext_var, comp_indep_var
  type (an_indep_var)  thevar
  type (lat_struct)  lat
  type (constraint_struct)  con
  type (db_struct) cesr
  type (pc_struct) pc
  type (twiss_max_struct) twiss_max
  type (pretz_struct)  pretz
  type (global_struct)  global
  type (zquad_struct)  quad
  type (coord_struct), allocatable :: co_a(:)
  type (nonlin_ele_struct)  nonlin
  type (moment_struct)  moment

  real(8) fig, opti
  real(8) v0(indep_var_maxx), vdel(indep_var_maxx), vec(indep_var_maxx)
  real(8) comp_vec(indep_var_maxx+10)
  real(8) fig_start, fig_end
  real(rp), pointer :: a_ptr

  integer n, gen, pop
  integer ix_ele, attrib
  integer num

  logical err

!
  if (indep_var%n_var < 1) then
    print *, 'Error in BMAD_OPT: No independent variables passed,'
    print *, '                   cannot continue.  Check $USE_LIST.'
    call err_exit
  endif


  call reallocate_coord( co_o, lat%n_ele_max )
  call allocate_nonlin_ele( nonlin_o, lat%n_ele_max )
  call allocate_moment( moment_o, lat%n_ele_max )

  gen = con%n_cycles
  pop = 5 * indep_var%n_var

  sp_o = con
  lat_o = lat
  cesr_o = cesr
  pc_o = pc
  twiss_max_o = twiss_max
  pretz_o = pretz
  global_o = global
  quad_o = quad
  nonlin_o = nonlin
  co_o = co_a
  indep_var_o = indep_var
  sext_var_o = sext_var

  do n = 1, indep_var%n_var 
    vec(n) = var_attribute(lat, indep_var%v(n))
    v0(n) = vec(n)
    vdel(n) = indep_var%v(n)%delta
  enddo

  call calcf( lat, con, pc, quad, co_a, global, twiss_max, pretz, nonlin,  &
       indep_var )
  call fom( lat, cesr, con, pc, quad, co_a, global, twiss_max, pretz,  &
       nonlin, moment, indep_var, sext_var )

  lat_o = lat

  fig_start = con%figure_of_merit
  call showme_topten(con, lat, 0)

  if (con%plot) call plotdo_bmadz('/xs', lat, co_a, pc)


  fig = opti (vec, gen,pop,indep_var%n_var,demerit,v0,vdel)
  print *, 'fig = ', fig

  do n = 1, indep_var%n_var
    thevar = indep_var_o%v(n)
    print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', vec(n)
    a_ptr => var_attribute(lat, indep_var_o%v(n))
    a_ptr = vec(n)
    
    ix_ele = indep_var_o%v(n)%rindex
    attrib = indep_var_o%v(n)%val_id
    call update_ele (vec(n), lat, ix_ele, attrib, num) 

  enddo

  if(con%sol_com >= 0 .and. con%sol_com <= 2) then
    call compensation_indep_var(con%sol_com, con%compensating_ele, lat_o, indep_var_o, comp_indep_var, comp_vec)
    
    do n = indep_var%n_var+1,2+con%sol_com
      thevar = comp_indep_var%v(n)
      print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', comp_vec(n)
      a_ptr => var_attribute(lat, comp_indep_var%v(n))
      a_ptr = comp_vec(n)
      
      ix_ele = comp_indep_var%v(n)%rindex
      attrib = comp_indep_var%v(n)%val_id
      call update_ele (comp_vec(n), lat, ix_ele, attrib, num) 
    enddo
  endif

  call lattice_bookkeeper(lat, err)
  call calcf(lat, con, pc, quad, co_a, global, twiss_max, pretz, nonlin, indep_var)

  call fom(lat, cesr, con, pc, quad, co_a, global, twiss_max, pretz,  nonlin, moment, indep_var, sext_var)

  fig_end = con%figure_of_merit

  print *, 'fig_start = ', fig_start
  print *, 'fig_end = ', fig_end


  return
end subroutine bmad_opt

!-----------------------------------------------------------------
!-----------------------------------------------------------------
!-----------------------------------------------------------------

real(8) function demerit(vec, end_flag)

  implicit none

  type (indep_var_struct) comp_indep_var
  type (an_indep_var)  thevar
  type (ele_struct)  ele
  type (ele_struct), pointer :: ele2
  type (control_struct), pointer :: ctl

  integer i, iv, ix
  integer :: n, iter = 0
  integer lun, end_flag
  integer ix_ele, attrib
  integer num

  real(8) :: vec(*), dem_min  = 1.0e35
  real(8) comp_vec(indep_var_maxx+3)
  real(rp), pointer :: a_ptr

  character*60 new_bmad_file
  character*16 attrib_name
  logical :: write_opt = .true.
  logical :: type_list = .true.
  logical err_flag
  logical err

!

  do n = 1, indep_var_o%n_var
    a_ptr => var_attribute(lat_o, indep_var_o%v(n))
    a_ptr = vec(n)
    ix_ele = indep_var_o%v(n)%rindex
    attrib = indep_var_o%v(n)%val_id
    call update_ele (vec(n), lat_o, ix_ele, attrib, num)
  enddo
  call lattice_bookkeeper(lat_o, err)
  call calcf(lat_o, sp_o, pc_o, quad_o, co_o, global_o,  &
       twiss_max_o, pretz_o, nonlin_o, indep_var_o, err_flag)
  call fom(lat_o, cesr_o, sp_o, pc_o, quad_o, co_o, global_o,  &
       twiss_max_o, pretz_o, nonlin_o, moment_o, indep_var_o, sext_var_o)

  demerit = sp_o%figure_of_merit
  if(demerit <= 0.99999*dem_min .or. (mod(iter,250)==0)) then
    print *, '****************************************************'
    print *, ' So far the minimum is ', min(dem_min, demerit)
    print *, '****************************************************'
    if (demerit <= 0.99999*dem_min) then
      dem_min = demerit
      if (write_opt) then
        lun = lunget()
        open(lun, file='opt.out')
        write(lun, *) 'demerit = ', demerit
        write(lun, *) 'num_variables =', indep_var_o%n_var
        do n = 1, indep_var_o%n_var
          theVar = indep_var_o%v(n)
          write(lun, '(4a, g16.8)')  &
               trim(theVar%ele_name), '[', trim(theVar%val_name), '] = ', vec(n)
        enddo

        do n = 1, indep_var_o%n_var
          ele = lat_o%ele(indep_var_o%v(n)%rindex)
          if (ele%lord_status == group_lord$) then
            write (lun, *)
            write (lun, *) '! Group: ', ele%name
            do i = 1, ele%n_slave
              ele2 => pointer_to_slave(ele, i, ctl)
              iv = ctl%ix_attrib ! pointer to attribute
              attrib_name = attribute_name(ele2, iv)
              write (lun, '(5a, 1pe16.6)') ' ! ',  &
                   trim(ele2%name), '[',  &
                   trim(attrib_name), '] =', ele2%value(iv)
            enddo
          endif
        enddo

        close(lun)
      endif

      if(err_flag) then
        print *, '****************************************************'
        print *, ' Warning!!! Bmad Error flag is set!                 '
        print *, '****************************************************'
      endif
      call showme_topten(sp_o, lat_o, 0)

      if (sp_o%plot) call plotdo_bmadz('/xs', lat_o, co_o, pc_o)


      if (sp_o%new_bmad_file<=' ') then
        new_bmad_file = 'new_bmad_input.'
      else
        new_bmad_file = sp_o%new_bmad_file
      endif


      if (sp_o%sol_com >= 0 .and. sp_o%sol_com <= 2) then
        comp_vec(1:indep_var_o%n_var)= vec(1:indep_var_o%n_var)
        call compensation_indep_var(sp_o%sol_com, sp_o%compensating_ele, &
             lat_o, indep_var_o, comp_indep_var, comp_vec)
        call opt_out(lat_o%input_file_name,new_bmad_file,comp_indep_var, &
             comp_vec(1:comp_indep_var%n_var))
      else
        call opt_out(lat_o%input_file_name,new_bmad_file,indep_var_o, vec)  
      endif
  

      if (sp_o%nonlinearity) call write_moments(lat_o, co_o, nonlin_o, moment_o)
        
    endif
    iter = 0
  endif

  iter = iter + 1

end function demerit

!-----------------------------------------------------------------
!-----------------------------------------------------------------
!-----------------------------------------------------------------

subroutine mrqmin_funcs(x,a,yfit,dyda)

  implicit none
  
  type (indep_var_struct) comp_indep_var
  type (an_indep_var)  thevar
  type (ele_struct)  ele

  real(rp), dimension(:), intent(in) :: x,a
  real(rp), dimension(:,:), intent(out) :: dyda
  real(rp), dimension(:), intent(out) :: yfit
  real(rp), pointer :: a_ptr

  integer n
  integer ma
  integer ix_ele, attrib, num

  logical err  
!

  ma = indep_var_o%n_var
  do n = 1, indep_var_o%n_var
    a_ptr => var_attribute(lat_o, indep_var_o%v(n))
    a_ptr = a(n)
    ix_ele = indep_var_o%v(n)%rindex
    attrib = indep_var_o%v(n)%val_id
    call update_ele(a(n), lat_o, ix_ele, attrib, num)
  enddo
  call lattice_bookkeeper(lat_o, err)
  call calcf(lat_o, sp_o, pc_o, quad_o, co_o, global_o,  &
       twiss_max_o, pretz_o, nonlin_o, indep_var_o)
  call fom(lat_o, cesr_o, sp_o, pc_o, quad_o, co_o, global_o,  &
       twiss_max_o, pretz_o, nonlin_o, moment_o, indep_var_o, sext_var_o)

  yfit(1:sp_o%n_constraint) = sp_o%c(1:sp_o%n_constraint)%contribution

  call bmad_calc_dyda(lat_o, sp_o, pc_o, quad_o, co_o, global_o, &
       twiss_max_o, pretz_o, cesr_o, nonlin_o, indep_var_o, &
       dyda, moment_o, sext_var_o)

  return
end subroutine mrqmin_funcs

!-----------------------------------------------------------------
!-----------------------------------------------------------------
!-----------------------------------------------------------------

  subroutine frprmn_opt (lat, indep_var, cesr, con, pc, quad,  &
       co, global, twiss_max, pretz, nonlin, moment, sext_var, ftol, iter, fret)

    implicit none

    type (indep_var_struct)  indep_var, sext_var
    type (indep_var_struct)  comp_indep_var
    type (an_indep_var)  thevar
    type (lat_struct)  lat
    type (constraint_struct)  con
    type (db_struct) cesr
    type (pc_struct) pc
    type (twiss_max_struct) twiss_max
    type (pretz_struct)  pretz
    type (global_struct)  global
    type (zquad_struct)  quad
    type (coord_struct), allocatable :: co(:)
    type (nonlin_ele_struct)  nonlin
    type (moment_struct)  moment

!  external demerit
    real(rp) v0(indep_var_maxx), vdel(indep_var_maxx)
    real(rp) fig_start, fig_end
    real(rp) ftol, fret, p(indep_var_maxx)
    real(rp) fig
    real(rp) comp_vec(indep_var_maxx+3)
    real(rp), pointer :: a_ptr

    integer n, iter, ix_ele, attrib, num

    logical err 

    character*60 new_bmad_file
!

!  gen = con%n_cycles

    call reallocate_coord( co_o, lat%n_ele_max )
    call allocate_nonlin_ele( nonlin_o, lat%n_ele_max )
    call allocate_moment( moment_o, lat%n_ele_max )

    sp_o = con
    lat_o = lat
    cesr_o = cesr
    pc_o = pc
    twiss_max_o = twiss_max
    pretz_o = pretz
    global_o = global
    quad_o = quad
    nonlin_o = nonlin
    co_o = co
    indep_var_o = indep_var
    sext_var_o = sext_var
    
    do n = 1, indep_var%n_var
      p(n) = var_attribute(lat, indep_var%v(n))
      v0(n) = p(n)
      vdel(n) = indep_var%v(n)%delta
    enddo

    call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin,  &
         indep_var)
    call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz,  &
         nonlin, moment, indep_var, sext_var)

    lat_o = lat

    fig_start = con%figure_of_merit
    call showme_topten(con, lat, 0)
    
    if (con%plot) call plotdo_bmadz('/xs', lat, co, pc)

    call frprmn_z(p(1:indep_var%n_var), ftol, iter, fig)

    fret = fig
    print *, 'fig = ', fig

    do n = 1, indep_var%n_var
      thevar = indep_var_o%v(n)
      print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', p(n)
      a_ptr => var_attribute(lat, indep_var_o%v(n))
      a_ptr = p(n)
      ix_ele = indep_var_o%v(n)%rindex
      attrib = indep_var_o%v(n)%val_id
      call update_ele (p(n), lat, ix_ele, attrib, num) 
    enddo

    if (con%sol_com >= 0 .and. con%sol_com <= 2) then
      comp_vec(1:indep_var%n_var)= p(1:indep_var%n_var)
      call compensation_indep_var(con%sol_com,con%compensating_ele, lat, indep_var, comp_indep_var, comp_vec)
      do n = 1, indep_var%n_var+1,2+con%sol_com
        thevar = comp_indep_var%v(n)
        print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', p(n)
        a_ptr => var_attribute(lat, indep_var_o%v(n))
        a_ptr = p(n)
        ix_ele = comp_indep_var%v(n)%rindex
        attrib = comp_indep_var%v(n)%val_id
        call update_ele (comp_vec(n), lat, ix_ele, attrib, num) 
      enddo
    endif
  call lattice_bookkeeper(lat, err)
    call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var)
    
    call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz,  nonlin, moment, indep_var, sext_var)


    if (con%new_bmad_file<=' ') then
      new_bmad_file = 'new_bmad_input.'
    else
      new_bmad_file = con%new_bmad_file
    endif

    if (con%sol_com >= 0 .and. con%sol_com <= 2) then
      comp_vec(1:indep_var%n_var)= p(1:indep_var%n_var)
      call compensation_indep_var(con%sol_com,con%compensating_ele, &
           lat, indep_var, comp_indep_var, comp_vec)
      call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, &
           comp_vec(1:comp_indep_var%n_var))
    else
      call opt_out(lat%input_file_name,new_bmad_file,indep_var, p)  
    endif



    fig_end = con%figure_of_merit

    print *, 'fig_start = ', fig_start
    print *, 'fig_end = ', fig_end
    
    return
  end subroutine frprmn_opt

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------


  subroutine marquardt_opt (lat, indep_var, cesr, con, pc, quad,  &
       co, global, twiss_max, pretz, nonlin, moment, sext_var)

    implicit none
    
    type (indep_var_struct)  indep_var, sext_var
    type (indep_var_struct)  comp_indep_var
    type (an_indep_var)  thevar
    type (lat_struct)  lat
    type (constraint_struct)  con
    type (db_struct) cesr
    type (pc_struct) pc
    type (twiss_max_struct) twiss_max
    type (pretz_struct)  pretz
    type (global_struct)  global
    type (zquad_struct)  quad
    type (coord_struct), allocatable :: co(:)
    type (nonlin_ele_struct)  nonlin
    type (moment_struct)  moment

    real(rp) demerit
    real(rp) v0(indep_var_maxx), vdel(indep_var_maxx)
    real(rp) fig_start, fig_end
    integer n
    real(rp) p(indep_var_maxx)
    real(rp) fig
    real(rp) comp_vec(indep_var_maxx+3)
    
    real(rp), allocatable, save :: x(:), y(:), sig(:), a(:)
    real(rp), dimension(:,:), allocatable, save :: covar, alpha
    real(rp), save :: chisq, alambda
    real(rp), pointer :: a_ptr

    integer i
    integer,save :: ma
    integer,save ::loops
    integer ix_ele, attrib, num
    logical, allocatable, save :: maska(:)
    logical,save :: init_needed = .true., finished
    logical err

    character*60 new_bmad_file
!

!  gen = con%n_cycles

    call reallocate_coord( co_o, lat%n_ele_max )
    call allocate_nonlin_ele( nonlin_o, lat%n_ele_max )
    call allocate_moment( moment_o, lat%n_ele_max )

 
    sp_o = con
    lat_o = lat
    cesr_o = cesr
    pc_o = pc
    twiss_max_o = twiss_max
    pretz_o = pretz
    global_o = global
    quad_o = quad
    nonlin_o = nonlin
    co_o = co
    indep_var_o = indep_var
    sext_var_o = sext_var
    
    do n = 1, indep_var%n_var
      p(n) = var_attribute(lat, indep_var%v(n))
      v0(n) = p(n)
      vdel(n) = indep_var%v(n)%delta
    enddo


    call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin,  &
                                                                   indep_var)
    call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz,  &
                                     nonlin, moment, indep_var, sext_var)

    lat_o = lat

    fig_start = con%figure_of_merit
    call showme_topten(con, lat, 0)

    if (con%plot) call plotdo_bmadz('/xs', lat, co, pc)

    finished = .false.

    if (init_needed) then
      loops = 0
      n = con%n_constraint
      ma = indep_var%n_var
      if (any(con%c(1:n)%weight == 0.) .or. ma == 0) then
        do i=1,n
          if(con%c(i)%weight == 0) &
               print '(a12,1x,i3,1x,a16)',' Constraint ', i, ' has zero weight'
        end do
        if(any(con%c(1:n)%weight == 0))print *,' Marquardt requires all nonzero weights'
        if(ma == 0)print *,' Number of variables is zero'
        stop
      endif
      if (allocated(x)) deallocate(x,y,sig,a,covar,alpha,maska)
      allocate(x(1:n), y(1:n), sig(1:n))
      allocate(a(1:ma), maska(1:ma))
      allocate(covar(1:ma,1:ma), alpha(1:ma,1:ma))


      forall (i = 1:con%n_constraint) x(i) = i 
      y(1:con%n_constraint) = 0.
      sig(1:con%n_constraint) = 1./con%c(1:con%n_constraint)%weight
      a(1:indep_var%n_var) = p(1:indep_var%n_var)
      alambda = -1
      forall (i = 1:ma) maska(i) = .true.
      init_needed = .false.

    endif

    if (alambda > 1e10) then
      print *, 'MARQUARDT: AT MINIMUM'
      finished = .true.
      init_needed = .true.
      alambda = 0
      call mrqmin(x, y, sig, a, maska, covar, alpha, chisq, mrqmin_funcs, alambda)
    else

      loops = loops +1
      if (loops == con%n_loops) alambda = 0
      call mrqmin(x, y, sig, a, maska, covar, alpha, chisq, mrqmin_funcs, alambda)
    endif
  
    print *, 'chisq = ', chisq
    
    do n = 1, indep_var%n_var
      thevar = indep_var_o%v(n)
      print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', a(n)
      a_ptr => var_attribute(lat, indep_var_o%v(n))
      a_ptr = a(n)
      if(con%sol_com >=0)then
        ix_ele = comp_indep_var%v(n)%rindex
        attrib = comp_indep_var%v(n)%val_id
        call update_ele(comp_vec(n), lat, ix_ele, attrib, num) 
      endif
    enddo

    if (con%sol_com >= 0 .and. con%sol_com <= 2) then
      comp_vec(1:indep_var%n_var)= a(1:indep_var%n_var)
      call compensation_indep_var(con%sol_com, con%compensating_ele, &
           lat, indep_var, comp_indep_var, comp_vec)
      do n = indep_var%n_var+1,2+con%sol_com
        thevar = comp_indep_var%v(n)
        print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', comp_vec(n)
        a_ptr => var_attribute(lat, comp_indep_var%v(n))
        a_ptr = comp_vec(n)
        
        ix_ele = comp_indep_var%v(n)%rindex
        attrib = comp_indep_var%v(n)%val_id
        call update_ele (comp_vec(n), lat, ix_ele, attrib, num) 
      enddo
      
    endif 

     call lattice_bookkeeper(lat,err)


    call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, &
         nonlin, indep_var)

    call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz,  &
         nonlin, moment, indep_var, sext_var)


    if (con%new_bmad_file<=' ') then
      new_bmad_file = 'new_bmad_input.'
    else
      new_bmad_file = con%new_bmad_file
    endif

    if (con%sol_com >= 0 .and. con%sol_com <= 2) then
      comp_vec(1:indep_var%n_var)= a(1:indep_var%n_var)
      call compensation_indep_var(con%sol_com, con%compensating_ele, &
           lat, indep_var, comp_indep_var, comp_vec)
      call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, &
           comp_vec(1:comp_indep_var%n_var))
    else
      call opt_out(lat%input_file_name,new_bmad_file,indep_var, a)  
    endif



    fig_end = con%figure_of_merit

    print '(1x,a12,e12.4)' , 'fig_start = ', fig_start
    print '(1x,a10,e12.4,6x,a11,e12.4,6x,a9,i3)', &
         'fig_end = ', fig_end, ' alambda = ',alambda,'loops = ',loops

    return
  end subroutine marquardt_opt

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------

  SUBROUTINE frprmn_z(p,ftol,iter,fret)
    IMPLICIT NONE

    INTEGER(I4B), INTENT(OUT) :: iter
    REAL(rp), INTENT(IN) :: ftol
    REAL(rp), INTENT(OUT) :: fret
    REAL(rp), DIMENSION(:), INTENT(INOUT) :: p
    INTERFACE
      FUNCTION opt_func(p)
        use precision_def
        USE nrtype
        IMPLICIT NONE
        REAL(rp), DIMENSION(:), INTENT(IN) :: p
        REAL(rp) :: opt_func
      END FUNCTION opt_func
!BL
      FUNCTION opt_dfunc(p)
        use precision_def
        USE nrtype
        IMPLICIT NONE
        REAL(rp), DIMENSION(:), INTENT(IN) :: p
        REAL(rp), DIMENSION(size(p)) :: opt_dfunc
      END FUNCTION opt_dfunc
    END INTERFACE
    INTEGER(I4B), PARAMETER :: ITMAX=200
    REAL(rp), PARAMETER :: EPS=1.0e-10
    INTEGER(I4B) :: its
    REAL(rp) :: dgg,fp,gam,gg
    REAL(rp), DIMENSION(size(p)) :: g,h,xi
    fp=opt_func(p)
    xi=opt_dfunc(p)
    g=-xi
    h=g
    xi=h
    do its=1,ITMAX
      iter=its
      call linmin(p,xi,fret, opt_func)
!      call linmin(p,xi,fret)
      if (2.0*abs(fret-fp) <= ftol*(abs(fret)+abs(fp)+EPS)) RETURN
      fp=fret
      xi=opt_dfunc(p)
      gg=dot_product(g,g)
      !    dgg=dot_product(xi,xi)
      dgg=dot_product(xi+g,xi)
      if (gg == 0.0) RETURN
      gam=dgg/gg
      g=-xi
      h=g+gam*h
      xi=h
    end do
    call nrerror('frprmn_z: maximum iterations exceeded')
  END SUBROUTINE frprmn_z

end module opt_com_mod