program lbfgs3D
! 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_wei_areal
use OW_mpi
use steepest_mod
use base_spline_mod
use solver_lbfgsb_mod


implicit none

real                                          :: start, finish, start1, finish1, start2, finish2
integer                                       :: ni,i,ishot,j,k,titer,stat1
integer                                       :: ir,is,iww
complex,allocatable                           :: ws(:,:,:),wr(:,:,:)
real, allocatable                             :: velint(:),vel0(:),pick1(:),velout(:)
integer                                       :: mpierr,i_r,ierr,status,shtida
logical                                       :: scale_cond

real                                          :: minxx
double precision                              :: fn,gn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Variables of L-BFGSC
!
 integer                       :: nmax ! nmax - dimension of the problem and mmax - max.# of ltd.memory corrections.
 
!     Declare the variables needed by the code.
!     A description of all these variables is given at the end of 
!     the driver.
 
 character*70                  :: task, csave
 logical                       :: lsave(4),cond
 integer                       :: iprint,subiter
 integer,allocatable           :: nbd(:), wwa(:)
 integer                       :: isave(44)
 double precision              :: factr, pgtol
 double precision              :: dsave(29)
 double precision, allocatable :: x(:), l(:), u(:), g(:),wa(:)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 call MPI_INIT(mpierr)
 if (mpierr .ne. MPI_SUCCESS) then
    write(0,*) 'Error starting MPI program. Terminating.'
    call MPI_ABORT(MPI_COMM_WORLD, wpar%i_mpi, mpierr)
 end if
 call sep_init()
 call owwei_mpi_init() 
 call set_no_putch()
 if (wpar%i_mpi==0) call owparam()
 call owparm_bcast1()
write(0,*) "Main: Parameter broadcast"

 allocate(x(img%z%n*img%xm%n*img%ym%n));x=dble(0.)
 allocate(wa(img%z%n*img%xm%n*img%ym%n));wa=dble(0.)
 allocate(wwa(img%z%n*img%xm%n*img%ym%n));wwa=0
 allocate(nbd(img%z%n*img%xm%n*img%ym%n));nbd=0
 allocate(g(img%z%n*img%ym%n*img%xm%n));g=dble(0.)
 allocate(pick1(img%xm%n*img%ym%n))
 allocate(vel0(nzold*img%ym%n*img%xm%n));vel0=0
 allocate(velout(nzold*img%ym%n*img%xm%n));velout=0
 nmax=img%z%n*img%ym%n*img%xm%n
write(0,*) "Main: vectors allocated"
 call init_files()
write(0,*) "Main: files initialized"

 if (wpar%i_mpi==0) then
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Read velocities
    allocate(velint(v%z%n*v%xm%n*v%ym%n));velint=0
    call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
    call owgrab_velim(velint,vel0,ni)
    deallocate(velint)
    !----------------------------------------------
    ! Horizon with no velocity update
write(0,*) "Before reading horizon"

    if (exist_file("horv0")) then
       open(51,file=hor1,status='old')
       do i=1,img%xm%n*img%ym%n
          read(51,*,end=111) pick1(i)
       end do
111    close(51)
    else
       pick1=0
    end if
write(0,*) "Before sending picks"
    do i_r=1,wpar%n_mpi-1
       write(0,*) "SENDING picks and velocity to node ",i_r
       call MPI_SEND(pick1, img%xm%n*img%ym%n,MPI_REAL,i_r,105,MPI_COMM_WORLD,ierr)
       call MPI_SEND(vel0, nzold*img%ym%n*img%xm%n,MPI_REAL,i_r,110,MPI_COMM_WORLD,ierr)
    end do
 else
    call MPI_RECV(pick1, img%xm%n*img%ym%n,MPI_REAL,0,105,MPI_COMM_WORLD,status,ierr)
    call MPI_RECV(vel0, nzold*img%ym%n*img%xm%n,MPI_REAL,0,110,MPI_COMM_WORLD,status,ierr)
    write(0,*) "Receiving picks and velocity at node ",wpar%i_mpi
 end if

 allocate(ws(sou%hx%n,sou%hy%n,wblock%ns),wr(sou%hx%n,sou%hy%n,wblock%ns))
! SEND/RECV SHOTS

    do is=1,shtblk
       call owwei_init(is,1,1)
       if (is == shtblk) then
          if (shtres.eq.0) then
             shtida = shtblm
          else
             shtida = shtres
          end if
       else
          shtida = shtblm
       end if
       do iww=1, wblock%nb
          if (wpar%i_mpi==0) then
             do i_r=1,shtida
                ws=0;wr=0
                write(0,*) "Sent shot ",(is-1)*shtblm+i_r, " to node ", i_r," freq.block = ",iww
                call gathers_io(iww,(is-1)*shtblm+i_r,wr(:,:,1:wblock%blen(iww)),"rec","reed")
                call gathers_io(iww,(is-1)*shtblm+i_r,ws(:,:,1:wblock%blen(iww)),"sou","reed")
                call MPI_SEND(wr,sou%hx%n*sou%hy%n*wblock%blen(iww),MPI_COMPLEX,i_r,120,MPI_COMM_WORLD,ierr)
                call MPI_SEND(ws,sou%hx%n*sou%hy%n*wblock%blen(iww),MPI_COMPLEX,i_r,130,MPI_COMM_WORLD,ierr)
                call MPI_SEND(i_r,1,MPI_INTEGER,i_r,131,MPI_COMM_WORLD,ierr)
             end do
          else
             ws=0;wr=0
             call MPI_RECV(wr,sou%hx%n*sou%hy%n*wblock%blen(iww),MPI_COMPLEX,0,120,MPI_COMM_WORLD,status,ierr)
             call MPI_RECV(ws,sou%hx%n*sou%hy%n*wblock%blen(iww),MPI_COMPLEX,0,130,MPI_COMM_WORLD,status,ierr)
             call MPI_RECV(i_r,1,MPI_INTEGER,0,131,MPI_COMM_WORLD,status,ierr)
             call gathers_io(iww,is,wr(:,:,1:wblock%blen(iww)),"rec","rite")
             call gathers_io(iww,is,ws(:,:,1:wblock%blen(iww)),"sou","rite")
             write(0,*) "RECV shot ",(is-1)*shtblm+i_r," in node ", wpar%i_mpi," freq.block = ",iww
          end if
       end do
    end do
 deallocate(wr,ws)
write(0,*) "==================================="
 call MPI_BARRIER(MPI_COMM_WORLD,ierr)
write(0,*) "==================================="
write(0,*) "Beginning of the optimization loops"
write(0,*) "==================================="
!--------------------------------------------
 ! Initial velocity model = background velocity
 do k=nzmin,nzmax;do j=1,img%ym%n;do i=1,img%xm%n
    x((k-nzmin)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i)=dble(vel0((k-1)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i))
 end do;end do;end do

 minxx=sngl(minval(x))
 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.
    if (cbounds) then
       l=vmin
       u=vmax
    else
       l=p_lo*x
       u=p_up*x
    end if
    nbd=2
 else
    nbd=0
 end if 

 call steep_mpi_init(pick1,minxx)
 if (scaling) call scale_init()


 iters=1
 titer=1
 task = 'START'
 cond=.true.
 iprint = 1
 factr=1.0d+7
 pgtol=1.0d-5

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Main code loop starts here  
!! 
 scale_cond=.true.
 do while(cond .and. iters <= max_iters)
     cond=.false.
     call MPI_BARRIER(MPI_COMM_WORLD,ierr)
     call setulb(nmax,mmax,x,l,u,nbd,fn,g,factr,pgtol,wa,wwa,task,iprint,csave,lsave,isave,dsave)
     if (task(1:2) .eq. 'FG') then
!        the minimization routine has returned to request the
!        function fun and gradient g values at the current x.
!        Compute function value fun for the sample problem.
!
!     x is the updated slowness
!
        cond=.true.

        do k=nzmin,nzmax;do j=1,img%ym%n;do i=1,img%xm%n
           vel0((k-1)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i)=sngl(x((k-nzmin)*img%ym%n*img%xm%n+(j-1)*img%xm%n+i))
        end do;end do;end do
        call to_history("n4",titer,"Velout")
        call srite("Velout",vel0,4*nzold*img%xm%n*img%ym%n)
        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
        fn=0;gn=0.;g=0.
        call fgeval_init(iters,sngl(x))
        call init_files()
        call scale_init1(scale_cond)
        call fgeval(vel0=sngl(x),ff=fn,gg=gn,grad=g)
        call to_history("n4",titer,"grad_")
        call grad_io("grad_","rite",titer,sngl(g))
        call MPI_BARRIER(MPI_COMM_WORLD,ierr)
        scale_cond=.false.
        titer=titer+1
     end if 
     if (task(1:5) .eq. 'NEW_X')  then
        write(0,113) " iter = ",iter," task = ",task(1:5)," J(x) = ", sngl(fn)
        iters=iters+1
        cond=.true.
     end if
 end do

01  format(a9,i7,a6,i9)
02  format(a9,i7,2(a6,i9))
03  format(a9,i7)
04  format(a12,i4,a7,i11)
113 format(a8,i4,a8,a5,a8,E)
 if (scaling) call scale_close()
 if (allocated(x)) deallocate(x)
 if (allocated(g)) deallocate(g)
 if (allocated(l)) deallocate(l)
 if (allocated(u)) deallocate(u)
 if (allocated(wa)) deallocate(wa)
 if (allocated(wwa)) deallocate(wwa)
 if (allocated(nbd)) deallocate(nbd)
 if (allocated(g)) deallocate(g)
 if (allocated(vel0)) deallocate(vel0)
 call steep_mpi_close()
 call sep_close()
 call MPI_FINALIZE(mpierr)
 stop
end program
