module oc_mpi_solverreg_mod
  !!=======================
  !! Paul Sava    (C) 2000
  !! paul@sep.stanford.edu
  !! oclib v2.0
  !!=======================
  use oc_global_mod
  use oc_tags_mod
  use oc_utils_mod
  use oc_file_mod
  use oc_scale_mod
  use oc_filealgebra_mod
  implicit none

  integer, private :: ocniter,ocmaxmem
  real,    private :: oceps
  logical, private :: ocverb, ocmmovie, ocdmovie, ocwght
  integer, private :: ocimpi
contains
  !! L = linear         operator
  !! A = regularization operator
  !! W = weighting      operator
  !! S = step

  !----------------------------------------------------------------
  subroutine oc_mpi_solverreg_init(impi,niter,eps,maxmem,verb,mmovie,dmovie,wght)
    integer           :: impi
    integer, optional :: niter, maxmem
    real,    optional :: eps
    logical, optional :: verb, mmovie, dmovie, wght

    ocimpi = impi

    ocniter=10;     if(present(niter))     ocniter=niter
    ocmaxmem=0;     if(present(maxmem))   ocmaxmem=maxmem
    ocverb=F;       if(present(verb))       ocverb=verb
    ocmmovie=F;     if(present(mmovie))   ocmmovie=mmovie
    ocdmovie=F;     if(present(dmovie))   ocdmovie=dmovie
    oceps=0.01;     if(present(eps))         oceps=eps
    ocwght=F;       if(present(wght))       ocwght=wght

    if(ocimpi==0) then
       call ocline()
       write(0,*) "niter=",ocniter
       write(0,*) "  eps=",oceps
       write(0,*) " wght=",ocwght
       call ocline()
    end if
  end subroutine oc_mpi_solverreg_init
  !----------------------------------------------------------------
  subroutine oc_mpi_solverreg( L,A,S, m_,d_, nreg, W ,k_,m0_,res_,op1,op2,op3,op4,op5,op6,op7,op8,op9)
    optional                                 :: W, k_,m0_,res_
    character(len=128),    intent(in)        ::    k_,m0_,res_
    character(len=128),    intent(in)        :: m_, d_
    integer, external,     optional          :: op1,op2,op3,op4,op5,op6,op7,op8,op9
    integer, dimension(:), intent(in)        :: nreg

    type(fileinfo)                           :: file
    character(len=128)                       ::    rd_,rm_ !! residual
    character(len=128)                       :: g_,gd_,gm_ !! gradient and conjugate gradient
    character(len=128)                       :: s_,sd_,sm_ !! previous step
    character(len=128)                       ::mm_,md_     !! model/data movie
    character(len=128)                       :: name
    character(len=128)                       ::    tm_,td_

    integer                                  :: i, stat
    logical                                  :: forget
    real, dimension(1)                       :: res

    interface
       integer function L( adj,add, m_,d_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(st)
         integer, external, optional     :: op1,op2,op3,op4,op5,op6,op7,op8,op9
         logical,           intent(in)   :: adj, add
         character(len=128),intent(in)   :: m_, d_
       end function L
       integer function A( adj,add, m_,d_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(st)
         integer, external, optional     :: op1,op2,op3,op4,op5,op6,op7,op8,op9
         logical,           intent(in)   :: adj, add
         character(len=128),intent(in)   :: m_, d_
       end function A
       integer function W( adj,add, m_,d_, op1,op2,op3,op4,op5,op6,op7,op8,op9) result(st)
         integer, external, optional     :: op1,op2,op3,op4,op5,op6,op7,op8,op9
         logical,           intent(in)   :: adj, add
         character(len=128),intent(in)   :: m_, d_
       end function W
       integer function S( forget,m_,g_,rr_,gg_,s_,ss_,maxmem) result(st)
         integer, optional,               intent(in) :: maxmem
         logical,                         intent(in) :: forget
         character(len=128),              intent(in) ::  m_,  g_,  s_
         character(len=128), dimension(:),intent(in) :: rr_, gg_, ss_
       end function S
    end interface

    if(ocimpi==0) then
       call oc_allocatefile(file,m_,ocmaxmem)

       name=".RS_Rm";    rm_ = oc_allocate(F, name, file%esize, nreg) !! model residual
       name=".RS_Rd";    rd_ = oc_clone(F,  d_,name, ocmaxmem) !! data  residual

       name=".RS_g";      g_ = oc_clone(F,  m_,name, ocmaxmem)
       name=".RS_Gm";    gm_ = oc_allocate(F, name, file%esize, nreg) !! model conj grad
       name=".RS_Gd";    gd_ = oc_clone(F,  d_,name, ocmaxmem) !! data  conj grad

       name=".RS_s";      s_ = oc_clone(F,  m_,name, ocmaxmem)
       name=".RS_Sm";    sm_ = oc_allocate(F, name, file%esize, nreg) !! model previous step
       name=".RS_Sd";    sd_ = oc_clone(F,  d_,name, ocmaxmem) !! data  previous step

       name=".RS_Td";     td_ = oc_clone(F, d_,name, ocmaxmem)
       name=".RS_Tm";     tm_ = oc_clone(F, m_,name, ocmaxmem)

       if(ocmmovie) then
          name=".RS_Mmovie"; mm_ = oc_clone(F, m_,name, ocmaxmem)
          call oc_adddim(mm_, ocniter)
       end if
       if(ocdmovie) then
          name=".RS_Dmovie"; md_ = oc_clone(F, d_,name, ocmaxmem)
          call oc_adddim(md_, ocniter)
       end if

       if( present(res_)) then
          call oc_shapeheader(res_,4,(/ocniter/),(/0./),(/1./))
       end if

       call oc_assign( rm_, real(rZERO), ocmaxmem)
       call oc_assign( sm_, real(rZERO), ocmaxmem)
       call oc_assign( gm_, real(rZERO), ocmaxmem)

       if(ocwght) then
          call message('Td = -  d')
          call oc_linear( td_, (/d_/), (/-rONE/), ocmaxmem)
          call message('Rd = W Td')
          stat = W( F,F, td_,rd_,op1,op2,op3,op4,op5,op6,op7,op8,op9)
       else
          call message('Rd = -  d')
          call oc_linear( rd_, (/d_/), (/-rONE/), ocmaxmem) 
       end if

       call oc_scale_init(oceps,ocmaxmem)
    end if

    if(present(m0_)) then                                           !! Rd = WL m0 - Wd
       if(ocwght) then                                              !! Rm = eA m0
          if(ocimpi==0) then
             if(ocverb) call message('Td = L  m0')
          end if
          stat =           L( F, F, m0_,td_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Rd = W  Td')
             stat =        W( F, F, td_,rd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Tm = A  m0')
             stat =        A( F, F, m0_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Rm =eI  Tm')
             stat = oc_scale( F, F, tm_,rm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       else
          if(ocimpi==0) then
             if(ocverb) call message('Rd = L  m0')
          end if
          stat =           L( F, F, m0_,rd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Tm = A  m0')
             stat =        A( F, F, m0_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Rm =eI  Tm')
             stat = oc_scale( F, F, tm_,rm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       end if
       if(ocimpi==0) call oc_linear( m_, (/m0_/),(/rONE/), ocmaxmem)
    else
       if(ocimpi==0) call oc_assign( m_ ,     real(rZERO), ocmaxmem)
    end if

    forget=T; if(present(m0_)) forget=F
    do i = 1, ocniter
       if(ocimpi==0 .and. ocverb) call message('iteration',i)
       !----------------------------------------------------------------       
       if(ocimpi==0 .and. ocverb) call message('g= L*W* Rd + eA* Rm')
       if(ocwght) then
          if(ocimpi==0) then
             if(ocverb) call message('Td = W* Rd')
             stat =        W( T, F, td_,rd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message(' g = L* Td')
          end if
          stat =           L( T, F,  g_,td_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Tm = A* Rm')
             stat =        A( T, F, tm_,rm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message(' g+=eA* Rm')
             stat = oc_scale( T, T,  g_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       else
          if(ocimpi==0) then
             if(ocverb) call message(' g = L* Rd')
          end if
          stat =           L( T, F,  g_,rd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Tm = A* Rm')
             stat =        A( T, F, tm_,rm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message(' g+=eA* Rm')
             stat = oc_scale( T, T,  g_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       end if
       !----------------------------------------------------------------
       if(ocimpi==0) then
          if(present(k_)) call oc_mask(k_,g_,ocmaxmem) !! don't change known
       end if
       !----------------------------------------------------------------       
       if(ocimpi==0 .and. ocverb) call message('Gd= W L g; Gm= eA g')
       if(ocwght) then
          if(ocimpi==0) then
             if(ocverb) call message('Td = L   g')
          end if
          stat =           L( F, F,  g_,td_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Gd = W  Td')
             stat =        W( F, F, td_,gd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Tm =eI   g')
             stat = oc_scale( F, F,  g_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Gm = A  Tm')
             stat =        A( F, F, tm_,gm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       else
          if(ocimpi==0) then
             if(ocverb) call message('Gd = L   g')
          end if
          stat =           L( F, F,  g_,gd_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          if(ocimpi==0) then
             if(ocverb) call message('Tm =eI   g')
             stat = oc_scale( F, F,  g_,tm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
             if(ocverb) call message('Gm = A  Tm')
             stat =        A( F, F, tm_,gm_, op1,op2,op3,op4,op5,op6,op7,op8,op9)
          end if
       end if
       !----------------------------------------------------------------
       if(ocimpi==0 .and. ocverb) call message('descend')
       if(ocimpi==0) stat = S(forget,m_,g_,(/rd_,rm_/),(/gd_,gm_/),s_,(/sd_,sm_/),ocmaxmem)
       !----------------------------------------------------------------
       if(ocimpi==0) then
          if(ocmmovie) then
             if(ocverb) call message('mmovie append')
             call oc_append(mm_,  m_)
          end if
          if(ocdmovie) then
             if(ocverb) call message('dmovie append')
             call oc_append(md_, rd_)
          end if
          !----------------------------------------------------------------
          if(present(res_)) then
             res = oc_norm((/rd_/), ocmaxmem)
             if( present(res_)) call sep_write(res, res_)
          end if
          !----------------------------------------------------------------
          if(stat ==1) then
             if( present(res_)) call oc_shapeheader(res_,4,(/i/),(/0./),(/1./))
             call message('DESCENT STOPPED')
             exit !! can't descend anymore
          end if
          !----------------------------------------------------------------
       end if
       forget = F !! forget=false except for i=1
    end do !! end iterations

  end subroutine oc_mpi_solverreg
  !----------------------------------------------------------------
end module oc_mpi_solverreg_mod


