program steepest_ls3D

! line_search = Option switch for line seaching (1=2-pt,0=simple) 
!      = 0 = none 
!      = 1 = simple 
!      = 2 = 2-point parabolic 
!      = 3 = 2-point parabolic 

use sep 
use OW_type
use OW_grabvel2_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use steepest_mod
use base_spline_mod
!use fx2d
use lloyd_vel

implicit none

character*80                                  :: hor1, listagem
real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: vmin,vmax,p_up,p_lo,tol(2)
integer                                       :: ni,i,ishot,j,k,iters,titer,stat1
real, allocatable                             :: velint(:),taper(:,:),pick1(:)
real, allocatable                             :: vel0(:)
logical                                       :: cond,bounds,cbounds
double precision, allocatable                 :: x(:), g(:)
real, allocatable                             :: l(:), u(:)

integer                                       :: rect1, rect2, rect3
integer                                       :: numtasks, rank, rc,mpierr

integer                                       :: max_iters, line_search, a
real                                          :: err, err1, epsteep, factor, lambda_min, tola, nrm, hh
real                                          :: lambda,lam1,lam2,lam3,lamc,lam4,gamma
double precision                              :: ft,fn,f2,gt,gn,g2,h20,h320,h32
double precision, allocatable                 :: xn(:), dx(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! include '/opt/MPI/mpich-1.2.7p1/include/mpif.h'
 call from_param("list",listagem,"listagem")
 open(unit=3,file=listagem)

 call sep_init()
 call sep_begin_prog()
 start=0.;finish=0.
 call owexec_time(start)
 call owparam()
 call init_files()
 call from_param("line_search",line_search,3)
 call from_param("max_iters",max_iters,10)
 call from_param("gamma",gamma,1.)

allocate(taper(kxm%n,kxm%n));taper=0.
allocate(g(img%z%n*img%ym%n*img%xm%n));g=dble(0.)
allocate(velint(v%z%n*v%xm%n*v%ym%n),vel0(img%z%n*img%ym%n*img%xm%n))
allocate(x(img%z%n*img%xm%n*img%ym%n));x=dble(0.)
allocate(xn(img%z%n*img%xm%n*img%ym%n));xn=dble(0.)
allocate(dx(img%z%n*img%xm%n*img%ym%n));dx=dble(0.)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(taper)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Read velocities
!
 call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
 call owgrab_velim(velint,vel0,ni)
 deallocate(velint)

!----------------------------------------------
! Set bounds
 call from_param("bounds",bounds,.false.)
 if (bounds) then
    allocate(l(img%z%n*img%ym%n*img%xm%n));l=0.
    allocate(u(img%z%n*img%ym%n*img%xm%n));u=0.
!    call from_param("cbounds",cbounds,.false.)
!    if (cbounds) then
!       call from_param("vmin",vmin)
!       call from_param("vmax",vmax)
!       l=1/vmax
!       u=1/vmin
!    else
       call from_param("p_up",p_up)
       call from_param("p_lo",p_lo)
       if (p_up >= 1.) p_up=1+p_up/100
       if (p_lo >= 1.) p_lo=1-p_lo/100
       l=p_lo*x
       u=p_up*x
!    end if
 end if 
!----------------------------------------------
! Horizon with no velocity update
if (exist_file("horv0")) then
   call from_param('horv0',hor1)
   open(51,file=hor1,status='old')
   allocate(pick1(img%xm%n*img%ym%n));pick1=9999.
   do i=1,img%xm%n*img%ym%n
      read(51,*,end=111) pick1(i)
   end do
111 close(51)
end if

! ********************************************** 
!  Initialization of code variables 
 iters=1
 err=1
 epsteep=1e-8 
 tol(1)=1e-8
 tol(2)=1e-5

! Line search variables 
 factor=0.5
 lambda_min=0.00001
! ********************************************** 
!----------------------------------------------
! Initial velocity model = background velocity
 x=dble(vel0)
 xn=dble(vel0)

 call to_history("n4",iters,"Velout")
 call srite("Velout",sngl(x),4*img%z%n*img%xm%n*img%ym%n)
!   Set tolerance for stopping criterium 
 err1=tiny(1.)
! ********************************************** 
 tola=0; 
! ********************************************** 
 do while (err > tola.and.iters <= max_iters) ! Main code loop starts here  
     if (exist_file("horv0")) then
       call fgeval_init(sngl(x),taper,iters,pick1)
     else
       call fgeval_init(sngl(x),taper,iters)
     end if
     call feval(fn,.true.)
     call geval(g)
     call fgeval_close()
!    Error 
     err=norm(g)
!    Set tolerance for stopping criterium 
     if (iters==1) tola=tol(1)*err+tol(2)
 
!    Update the Newton search direction 
     dx=-(g/maxval(abs(g)))*minval(x)
     write(0,*) minval(dx),maxval(dx)
!     call to_history("n4",iters,"grad")
!     call srite("grad",sngl(dx),4*img%z%n*img%xm%n)

! ****************************************************************  
! Line Searching below 
! ****************************************************************  
     lam3=1.
     gn=fn*fn
999    x=xn+lam3*dx
       if (any(mask=(x <= l.or.x >= u))) then
          lam3=factor*lam3
          write(0,*) "1) Bound reached - recomputing lam3", lam3
          goto 999
       end if
!      if (exist_file("horv0")) then
!         call ow_smoovel(x,vel0,pick1)
!         call oweditvel(pick1=pick1,vel=vel,vel0=vel0)
!      else
!         call ow_smoovel(x,vel0)
!         call oweditvel(vel=vel)
!     end if
     if (exist_file("horv0")) then
       call fgeval_init(sngl(x),taper,iters,pick1)
     else
       call fgeval_init(sngl(x),taper,iters)
     end if
     call feval(ft,.false.)
     call fgeval_close()
     gt=ft*ft
     do while (gt > gn.and.lam3 > lambda_min) 
1000      lam3=factor*lam3
          x=xn+lam3*dx
          if (any(mask=(x <= l.or.x >= u))) then
             write(0,*) "2) Bound reached - recomputing lam3", lam3
             goto 1000
          end if
!        if (exist_file("horv0")) then
!           call ow_smoovel(x,vel0,pick1)
!        else
!           call ow_smoovel(x,vel0)
!        end if
        if (exist_file("horv0")) then
          call fgeval_init(sngl(x),taper,iters,pick1)
        else
          call fgeval_init(sngl(x),taper,iters)
        end if
        call feval(ft,.false.)
        call fgeval_close()
        gt=ft*ft
        write(*,01) "01) fn = ",int(fn)," ft = ",int(ft)
     end do     
     
     lam2=factor*lam3; 
     x=xn+lam2*dx 
!     if (exist_file("horv0")) then
!        call ow_smoovel(x,vel0,pick1)
!     else
!        call ow_smoovel(x,vel0)
!     end if
     if (exist_file("horv0")) then
       call fgeval_init(sngl(x),taper,iters,pick1)
     else
       call fgeval_init(sngl(x),taper,iters)
     end if
     call feval(f2,.false.)
     call fgeval_close()
     g2=f2*f2
     h20=(g2-gn)/lam2
     h32=(gt-g2)/(lam3-lam2)
     h320=(h32-h20)/lam3
 
     f2=0.
     lamc=0.5*(lam2-h20/h320); 
     x=xn+lamc*dx
!     if (exist_file("horv0")) then
!        call ow_smoovel(x,vel0,pick1)
!     else
!        call ow_smoovel(x,vel0)
!     end if
     if (exist_file("horv0")) then
       call fgeval_init(sngl(x),taper,iters,pick1)
     else
       call fgeval_init(sngl(x),taper,iters)
     end if
     call feval(f2,.false.)
     call fgeval_close()
     g2=f2*f2
    
     write(*,02) "02) fn = ",int(fn)," ft = ",int(ft)," f2 = ",int(f2)
     if (g2 < gt) then
        lambda=lamc
        fn=f2
        gn=g2
     else 
        lambda=lam3
        fn=ft
        gn=gt
     end if   
 
! ****************************************************************** 
! Line search complete 
! ****************************************************************** 
     xn=xn+lambda*dx ! Update Newton iterate
     x=xn 
!     if (exist_file("horv0")) then
!        call ow_smoovel(x,vel0,pick1)
!     else
!        call ow_smoovel(x,vel0)
!     end if
     write(*,04) 'Iteration = ',iters,' error=',int(err)
     write(*,03) "03) fn = ",int(fn)
 
     iters=iters+1;            !Update iterations counter 
     call to_history("n4",iters,"Velout")
     call srite("Velout",sngl(x),4*img%z%n*img%xm%n*img%ym%n)
 
 end do

! Success statement based on error and iters ******** 
 if (err > err1.or.iters == max_iters) then
    write(0,*) 'Convergence may not have occurred; because either the error is greater than 20% grad_0 or max_iters was reached' 
 else if (err < err1) then
    write(0,*)  'The error is less than 20% grad_0; convergence has likely been achieved' 
 end if 

 if (allocated(taper)) deallocate(taper)
 if (allocated(g)) deallocate(g)
 if (allocated(x)) deallocate(x)
 if (allocated(xn)) deallocate(xn)
 if (allocated(dx)) deallocate(dx)
01 format(a9,i7,a6,i9)
02 format(a9,i7,2(a6,i9))
03 format(a9,i7)
04 format(a12,i4,a7,i11)
close (3)
 call sep_end_prog()

end program
