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_wei_areal
use OW_mpi
use steepest_mod
use base_spline_mod


implicit none

real                                          :: start, finish, start1, finish1, start2, finish2
integer                                       :: ni,i,ishot,j,k,titer,stat1
real, allocatable                             :: velint(:),vel0(:),pick1(:)
double precision, allocatable                 :: x(:), g(:)
real, allocatable                             :: l(:), u(:)
integer                                       :: mpierr,i_r,ierr,status

real                                          :: lambda,lam1,lam2,lam3,lamc,lam4
double precision                              :: ft,fn,f2,gt,gn,g2,h20,h320,h32
double precision, allocatable                 :: xn(:), dx(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 call MPI_INIT_(mpierr)
 call sep_init()
 call owwei_mpi_init()
 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 set_no_putch()
 call owparam()
 call owparm_bcast1()
 allocate(x(img%z%n*img%xm%n*img%ym%n));x=dble(0.)
 allocate(g(img%z%n*img%ym%n*img%xm%n));g=dble(0.)
 if (hor1(1:1).ne." ") then
    allocate(pick1(img%xm%n*img%ym%n));pick1=9999.
 end if
 call init_files()
 if (wpar%i_mpi==0) then
    start=0.;finish=0.
    call owexec_time(start)
    allocate(xn(img%z%n*img%xm%n*img%ym%n));xn=dble(0.)
    allocate(velint(v%z%n*v%xm%n*v%ym%n),vel0(img%z%n*img%ym%n*img%xm%n))
    allocate(dx(img%z%n*img%xm%n*img%ym%n));dx=dble(0.)

    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*1/x
          u=p_up*1/x
       end if
    end if 

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Read velocities
    call sreed("vel",velint,v%z%n*v%xm%n*v%ym%n*4)
    call owgrab_velim(velint,vel0,ni)
    deallocate(velint)
   !----------------------------------------------
   ! Initial velocity model = background velocity
    x=dble(vel0)
    xn=dble(vel0)
    call to_history("n4",iters,"Velout")
    call srite("Velout",vel0,4*img%z%n*img%xm%n*img%ym%n)
    deallocate(vel0)

    !----------------------------------------------
    ! Horizon with no velocity update
    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)
    end if

    do i_r=1,wpar%n_mpi-1
       write(0,*) wpar%i_mpi,i_r
       if (allocated(pick1)) then
!          write(0,*) "before sending pick1",i_r
          call MPI_SEND_(pick1, img%xm%n*img%ym%n,MPI_REAL,i_r,105,MPI_COMM_WORLD,ierr)
!          write(0,*) "after sending pick1",i_r
       end if
!       write(0,*) "before sending x",i_r
       call MPI_SEND_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,i_r,110,MPI_COMM_WORLD,ierr)
!       write(0,*) "after sending x",i_r
    end do
 else
    if (allocated(pick1)) then
!       write(0,*) "before receiving pick1",wpar%i_mpi
       call MPI_RECV_(pick1, img%xm%n*img%ym%n,MPI_REAL,0,105,MPI_COMM_WORLD,status,ierr)
!       write(0,*) "after receiving pick1",wpar%i_mpi
       call steep_mpi_init(pick1)
!       write(0,*) "after steep_mpi_init",wpar%i_mpi
    end if
!    write(0,*) "before receiving x",wpar%i_mpi
    call MPI_RECV_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,0,110,MPI_COMM_WORLD,status,ierr)
!    write(0,*) "after receiving x",wpar%i_mpi
 end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Main code loop starts here  
!! changed or to and
 do while (err > tola.or.iters <= max_iters) 
    call fgeval_init(iters)
!    call fgeval(vel0=sngl(x),ff=fn,gg=gn,grad=g,riteD=.true.)
    call fgeval(vel0=sngl(x),ff=fn,gg=gn,riteD=.true.);goto 1212
    call MPI_BCAST_(gn,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
    if (wpar%i_mpi==0) then
       dx=-(g/maxval(abs(g)))*minval(x) 	! Update the Newton search direction 
       err=norm(dx)				! Error 
       if (iters==1) tola=tol(1)*err+tol(2) 	! Set tolerance for stopping criterium 
       write(0,*) minval(dx),maxval(dx)
       ! ******************** 
       ! Line Searching below 
       ! ******************** 
       lam3=1.
       gn=fn*fn
900    x=xn+lam3*dx
       if (any(mask=(1/x <= l.or.1/x >= u))) then
          lam3=factor*lam3
          call MPI_BCAST_(lam3,1,MPI_REAL,0,MPI_COMM_WORLD,mpierr) 
          write(0,*) "1) Bound reached - recomputing lam3", lam3
          goto 900
       end if
       do i_r=1,wpar%n_mpi-1
          call MPI_SEND_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,i_r,160,MPI_COMM_WORLD,ierr)
       end do
    else
       call MPI_RECV_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,0,160,MPI_COMM_WORLD,status,ierr)
    end if
    call fgeval(vel0=sngl(x),ff=ft,gg=gt,riteD=.false.)
    call MPI_BCAST_(gt,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 

! changed and to or (16.01.10)
    do while (gt > gn.or.lam3 > lambda_min) 
       if (wpar%i_mpi==0) then
901       lam3=factor*lam3
          call MPI_BCAST_(lam3,1,MPI_REAL,0,MPI_COMM_WORLD,mpierr) 
          x=xn+lam3*dx
          if (any(mask=(1/x <= l.or.1/x >= u))) then
             write(0,*) "2) Bound reached - recomputing lam3", lam3
             goto 901
          end if
          do i_r=1,wpar%n_mpi-1
             call MPI_SEND_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,i_r,190,MPI_COMM_WORLD,ierr)
          end do
       else
          call MPI_RECV_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,0,190,MPI_COMM_WORLD,status,ierr)
       end if
       call fgeval(vel0=sngl(x),ff=ft,gg=gt,riteD=.false.)
       call MPI_BCAST_(gt,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
       if (wpar%i_mpi==0) write(*,01) "01) fn = ",int(fn)," ft = ",int(ft)
    end do     
     
     if (wpar%i_mpi==0) then
        lam2=factor*lam3; 
        x=xn+lam2*dx 
        do i_r=1,wpar%n_mpi-1
           call MPI_SEND_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,i_r,220,MPI_COMM_WORLD,ierr)
        end do
     else
        call MPI_RECV_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,0,220,MPI_COMM_WORLD,status,ierr)
     end if
     call fgeval(vel0=sngl(x),ff=f2,gg=g2,riteD=.false.)
     call MPI_BCAST_(g2,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
 
     if (wpar%i_mpi==0) then
        h20=(g2-gn)/lam2
        h32=(gt-g2)/(lam3-lam2)
        h320=(h32-h20)/lam3
        lamc=0.5*(lam2-h20/h320); 
        x=xn+lamc*dx
        do i_r=1,wpar%n_mpi-1
           call MPI_SEND_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,i_r,220,MPI_COMM_WORLD,ierr)
        end do
     else
        call MPI_RECV_(x, img%xm%n*img%ym%n*nz,MPI_DOUBLE_PRECISION,0,220,MPI_COMM_WORLD,status,ierr)
     end if
     call fgeval(vel0=sngl(x),ff=f2,gg=g2,riteD=.false.)
     call MPI_BCAST_(g2,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 

     if (wpar%i_mpi==0) then
        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 
        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 if  
 end do

 if (allocated(x)) deallocate(x)
 if (allocated(g)) deallocate(g)
 if (wpar%i_mpi==0) then
    ! Success statement based on error and iters ******** 
    if (err > err1.or.iters == max_iters) then
       write(0,*) 'Convergence may not have occurred; 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(xn)) deallocate(xn)
    if (allocated(dx)) deallocate(dx)
 end if
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)
1212  call steep_mpi_close()
 call sep_close()
 call MPI_FINALIZE_(mpierr)
end program
