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(:)
integer, external                             :: omp_get_thread_num,omp_get_num_threads
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

 write(0,*) "before MPI_INIT"
 call MPI_INIT_(mpierr)
    if (mpierr .ne. MPI_SUCCESS) then
    write(0,*) 'Error starting MPI program. Terminating.'
    call MPI_ABORT_(MPI_COMM_WORLD, rc, mpierr)
 end if
 call sep_init()
 call sep_begin_prog()
 call set_no_putch()

 call owwei_mpi_init()
 call owparam()
 if(wpar%i_mpi==0) then
    call MPI_SEP_SEND_ARGS(wpar%n_mpi,1,0)
 else
    call MPI_SEP_RECEIVE_ARGS()
 end if

 allocate(g(img%z%n*img%ym%n*img%xm%n));g=dble(0.)
 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.)
 if (wpar%i_mpi==0) then
    start=0.;finish=0.
    call owexec_time(start)
    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.)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! 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",sngl(x),4*img%z%n*img%xm%n*img%ym%n)

    do i_r=1,wpar%n_mpi-1
       call MPI_Send_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,i_r,110,MPI_COMM_WORLD,ierr)
    end do
 else
    call MPI_Recv_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,0,110,MPI_COMM_WORLD,status,ierr)
 end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Main code loop starts here  
! changed or to and
 do while (err > tola.or.iters <= max_iters) 
    if (wpar%i_mpi==0) then
       call sseek("rec",0,0)
       call sseek("sou",0,0)
       fn=0;gn=0;g=0;g_r=0;fn_r=0
    else
       fn=0;g=0;
    end if
    do idip=1,shtblk ! Reading in/transfering data
       if (shtblk - idip +1 <= shtres) then
          shtidx = wpar%n_mpi
       else
          shtidx = wpar%n_mpi-1
       end if
       if (wpar%i_mpi==0) then
          do i_r=1,shtidx
             wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
             ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
             call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,120,MPI_COMM_WORLD,ierr)
             call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,130,MPI_COMM_WORLD,ierr)
          end do
          shtblm=shtidx
       else
          call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,120,MPI_COMM_WORLD,status,ierr)
          call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,130,MPI_COMM_WORLD,status,ierr)
       end if
       if (wpar%i_mpi > 0) then
          call fgeval_init(iters)
          call feval_exe(ws,wr,sngl(x),taper,fn,.true.)
          call geval(ws,wr,sngl(x),taper,g)
          
!          call MPI_Send_(fn,                   1, MPI_DOUBLE_PRECISION,0,140,MPI_COMM_WORLD,ierr)
!          call MPI_Send_(g, img%xm%n*img%ym%n*nz, MPI_DOUBLE_PRECISION,0,150,MPI_COMM_WORLD,ierr)
          call MPI_Reduce_(fn, fn_r,                    1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
          call MPI_Reduce_( g,  g_r, img%xm%n*img%ym%n*nz, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       else
          g = g + g_r
          fn= fn+ fn_r
          gn = gn + fn*fn
          call MPI_BCAST_(gn,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
       end if
    end do

    if (wpar%i_mpi==0) then
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! Smoothing the gradient
       write(0,*) "main => ",maxval(g),minval(g)
       if (axkn.ne.99999 .or. axkn.ne.99999 .or. azkn.ne.99999) then
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
          call owexec_time(start2)
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          g=dble(0.)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
       end if
       write(0,*) "main => ",maxval(g),minval(g)
       write(0,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60
       if (exist_file("horv0")) call oweditvel(pick1,g)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call to_history("n4",iter,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n*img%ym%n)






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

! ****************************************************************  
! 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_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,i_r,160,MPI_COMM_WORLD,ierr)
       end do
    else
       call MPI_Recv_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,0,160,MPI_COMM_WORLD,status,ierr)
    end if
    if (wpar%i_mpi==0) then
       call sseek ("rec",0,0)
       call sseek ("sou",0,0)
       ft=0;gt=0;fn_r=0
    else
       fn=0
    end if
    do idip=1,shtblk ! Reading in/transfering data
       if (shtblk - idip +1 <= shtres) then
          shtidx = wpar%n_mpi
       else
          shtidx = wpar%n_mpi-1
       end if
       if (wpar%i_mpi==0) then
          do i_r=1,shtidx
             wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
             ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
             call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,170,MPI_COMM_WORLD,ierr)
             call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,180,MPI_COMM_WORLD,ierr)
          end do
          shtblm=shtidx
       else
          call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,170,MPI_COMM_WORLD,status,ierr)
          call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,180,MPI_COMM_WORLD,status,ierr)
       end if
       if (wpar%i_mpi > 0) then
          call fgeval_init(iters)
          call feval_exe(ws,wr,sngl(x),taper,fn,.false.)
          
          call MPI_Reduce_(fn, fn_r,                    1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
       else
          ft = ft+ fn_r
          gt = gt + ft*ft
          call MPI_BCAST_(gt,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
       end if
    end do

    do while (gt > gn.and.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_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,i_r,190,MPI_COMM_WORLD,ierr)
          end do
       else
          call MPI_Recv_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,0,190,MPI_COMM_WORLD,status,ierr)
       end if
       if (wpar%i_mpi==0) then
          call sseek ("rec",0,0)
          call sseek ("sou",0,0)
          ft=0;gt=0;fn_r=0
       else
          fn=0
       end if
       do idip=1,shtblk ! Reading in/transfering data
          if (shtblk - idip +1 <= shtres) then
             shtidx = wpar%n_mpi
          else
             shtidx = wpar%n_mpi-1
          end if
          if (wpar%i_mpi==0) then
             do i_r=1,shtidx
                wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
                ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
                call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,200,MPI_COMM_WORLD,ierr)
                call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,210,MPI_COMM_WORLD,ierr)
             end do
             shtblm=shtidx
          else
             call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,200,MPI_COMM_WORLD,status,ierr)
             call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,210,MPI_COMM_WORLD,status,ierr)
          end if
          if (wpar%i_mpi > 0) then
             call fgeval_init(iters)
             call feval_exe(ws,wr,sngl(x),taper,fn,.false.)
             
             call MPI_Reduce_(fn, fn_r,                    1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
          else
             ft = ft+ fn_r
             gt = gt + ft*ft
             call MPI_BCAST_(gt,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
          end if
       end do
       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_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,i_r,220,MPI_COMM_WORLD,ierr)
        end do
     else
        call MPI_Recv_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,0,220,MPI_COMM_WORLD,status,ierr)
     end if
     if (wpar%i_mpi==0) then
        call sseek ("rec",0,0)
        call sseek ("sou",0,0)
        f2=0;g2=0;fn_r=0
     else
        fn=0
     end if
     do idip=1,shtblk ! Reading in/transfering data
        if (shtblk - idip +1 <= shtres) then
           shtidx = wpar%n_mpi
        else
           shtidx = wpar%n_mpi-1
        end if
        if (wpar%i_mpi==0) then
           do i_r=1,shtidx
              wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
              ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
              call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,230,MPI_COMM_WORLD,ierr)
              call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,240,MPI_COMM_WORLD,ierr)
           end do
           shtblm=shtidx
        else
           call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,230,MPI_COMM_WORLD,status,ierr)
           call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,240,MPI_COMM_WORLD,status,ierr)
        end if
        if (wpar%i_mpi > 0) then
           call fgeval_init(iters)
           call feval_exe(ws,wr,sngl(x),taper,fn,.false.)
           
           call MPI_Reduce_(fn, fn_r,                    1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
        else
           f2 = f2+ fn_r
           g2 = g2 + f2*f2
           call MPI_BCAST_(g2,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
        end if
     end do
 
!! stopped from here
     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_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,i_r,220,MPI_COMM_WORLD,ierr)
        end do
     else
        call MPI_Recv_(sngl(x), img%xm%n*img%ym%n*nz,MPI_REAL,0,220,MPI_COMM_WORLD,status,ierr)
     end if

     if (wpar%i_mpi==0) then
        call sseek ("rec",0,0)
        call sseek ("sou",0,0)
        f2=0;g2=0;fn_r=0
     else
        fn=0
     end if
     do idip=1,shtblk ! Reading in/transfering data
        if (shtblk - idip +1 <= shtres) then
           shtidx = wpar%n_mpi
        else
           shtidx = wpar%n_mpi-1
        end if
        if (wpar%i_mpi==0) then
           do i_r=1,shtidx
              wr=0;   call gathers_io((idip-1)*shtblm+i_r,wr(ninx:nfnx,niny:nfny,:),"rec","reed")
              ws=0;   call gathers_io((idip-1)*shtblm+i_r,ws(ninx:nfnx,niny:nfny,:),"sou","reed")
              call MPI_Send_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,230,MPI_COMM_WORLD,ierr)
              call MPI_Send_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,i_r,240,MPI_COMM_WORLD,ierr)
           end do
           shtblm=shtidx
        else
           call MPI_Recv_(wr,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,230,MPI_COMM_WORLD,status,ierr)
           call MPI_Recv_(ws,kxm%n*kym%n*sou%w%n, MPI_COMPLEX,0,240,MPI_COMM_WORLD,status,ierr)
        end if
        if (wpar%i_mpi > 0) then
           call fgeval_init(iters)
           call feval_exe(ws,wr,sngl(x),taper,fn,.false.)
           
           call MPI_Reduce_(fn, fn_r,                    1, MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
        else
           f2 = f2+ fn_r
           g2 = g2 + f2*f2
           call MPI_BCAST_(g2,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,mpierr) 
        end if
     end do

     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

! 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
