module new_solver_mod
  use chain_mod
  use llist_mod
  use arnoldi_mod
  use ddot_mod

  implicit none
  logical, parameter, private  :: T = .true., F = .false.
  real,               private  :: scale
  !----------------------------------------------------------------
contains
  subroutine solver_prec(oper,step, prec,nprec, x,dat, niter,eps &
  ,             p0,rm0,err,resd,resm,xmov,rmov,verb,weight,kmask)
    optional :: p0,rm0,err,resd,resm,xmov,rmov,verb,weight,kmask
    interface
       integer function   oper(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function oper
       integer function   prec(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function prec
       integer function weight(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function weight
       integer function  kmask(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function kmask
       integer function step(forget,x,g,rr,gg) result(stat)
         logical              :: forget
         real, dimension(:)   :: x,g, rr,gg
       end function step
    end interface
    !----------------------------------------------------------------
    real, dimension(:),    intent(in)          :: dat, p0,rm0
    real, dimension(:),    intent(out)         :: x, err, resd, resm
    real,                  intent(in)          :: eps
    integer,               intent(in)          :: niter, nprec
    real, dimension(:,:),  intent(out)         :: xmov, rmov   
    logical,               intent(in)          :: verb

    real, dimension(size(x))                   :: p ,  g
    real, dimension(size(dat) + nprec), target :: rr, gg, tt
    real, dimension(:), pointer                :: rd, gd, td
    real, dimension(:), pointer                :: rm, gm, tm
    integer                                    :: iter, stat, nd
    logical                                    :: forget, v
    integer                                    :: dprd, dprm, dpx, dpg

    v = F ; if(present(verb)) v = verb
    nd = size(dat)

    rd => rr(1 : nd); rm => rr(1 + nd :)
    gd => gg(1 : nd); gm => gg(1 + nd :)
    td => tt(1 : nd); tm => tt(1 + nd :)

    if(present(weight)) then
       stat = weight(F,F,-dat,rd)  !! Rd = - W D
    else
       rd = -dat                   !! Rd = -   D
    end if
    rm = 0.                        !! Rm =     0
    if(present(rm0)) rm=rm0

    if(present(p0)) then
       stat =   prec(F,F,p0,tm)    !! t  =          P p0
       stat =   oper(F,F,tm,td)    !! t  =        L P p0
       if(present(weight)) then
          stat = weight(F,T,td,rd) !! Rd = Rd + W L P p0
       else
          rd = rd + td             !! Rd = Rd +   L P p0
       end if
       rm   = rm + eps * p0        !! Rm = Rm +   e I p0
       p    = p0
    else
       p = 0.
    end if

    forget = T; if(present(p0)) forget=F
    do iter = 1,niter

       if(present(weight)) then
          call chain(weight,oper,prec,T,F,g,rd,tm,td)
       else
          call chain(       oper,prec,T,F,g,rd,tm   )
       end if                      !! g =   (WLP)' Rd
       g = g + eps*rm              !! g = g+ e  I' Rm

       if(present(kmask)) then
          tm=g
          stat=kmask(F,F,tm,g)
       end if

       gm   =  eps*g               !! Gm =   e I g
       if(present(weight)) then
          call chain(weight,oper,prec,F,F,g,gd,tm,td)
       else
          call chain(       oper,prec,F,F,g,gd,tm   )
       end if                      !! Gd = (WLP) g

       stat = step(forget, p,g, rr,gg)
       if(stat ==1) then
          exit !! can't descend anymore
       end if
       stat = prec(F,F,p,x)         !! x = P p

       if(present(xmov)) xmov(:,iter) =  x(:size(xmov,1))
       if(present(rmov)) rmov(:,iter) = rr(:size(rmov,1))
       if(present(err ))  err(  iter) = dot_product(rd,rd)
       if(v) call solver_report(iter,x,g,rd,rm)

       forget=F
    end do

    if(present(resd)) resd = rd
    if(present(resm)) resm = rm(:size(resm))
  end subroutine solver_prec
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine solver_reg(oper,step, reg,nreg, x,dat, niter,eps &
  ,             x0,rm0,err,resd,resm,xmov,rmov,verb,weight,kmask)
    optional :: x0,rm0,err,resd,resm,xmov,rmov,verb,weight,kmask
    interface
       integer function   oper(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function oper
       integer function    reg(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function reg
       integer function weight(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function weight
       integer function  kmask(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function kmask
       integer function step(forget,x,g,rr,gg) result(stat)
         logical              :: forget
         real, dimension(:)   :: x,g, rr,gg
       end function step

    end interface
    !----------------------------------------------------------------
    real, dimension(:),    intent(in)          :: dat, x0,rm0
    real, dimension(:),    intent(out)         :: x, err, resd, resm
    real,                  intent(in)          :: eps
    integer,               intent(in)          :: niter, nreg
    real, dimension(:,:),  intent(out)         :: xmov, rmov
    logical,               intent(in)          :: verb

    real, dimension(size(x))                   :: g
    real, dimension(size(dat) + nreg), target  :: rr, gg, tt
    real, dimension(:), pointer                :: rd, gd, td
    real, dimension(:), pointer                :: rm, gm, tm
    integer                                    :: iter, stat, nd
    logical                                    :: forget, v
    integer                                    :: dprd, dprm, dpx, dpg

    v = F ; if(present(verb)) v = verb
    nd = size(dat)

    rd => rr(1 : nd); rm => rr(1 + nd :)
    gd => gg(1 : nd); gm => gg(1 + nd :)
    td => tt(1 : nd); tm => tt(1 + nd :)

    if(present(weight)) then
       stat = weight(F,F,-dat,rd)      !! Rd = - W D
    else
       rd = -dat                       !! Rd = -   D
    end if
    rm = 0.                            !! Rm =     0
    if(present(rm0)) rm=rm0

    if(present(x0)) then
       stat =   oper(F,F,    x0,td)    !!  t =        L x0

       if(present(weight)) then
          stat = weight(F,T,    td,rd) !! Rd = Rd + W L x0
       else
          rd = rd + td                 !! Rd = Rd +   L x0
       end if
       stat =    reg(F,T,eps*x0,rm)    !! Rm = Rm + e A x0
       x = x0
    else
       x = 0.
    end if

    forget = T; if(present(x0)) forget=F
    do iter = 1,niter

       if(present(weight)) then
          call chain(weight,oper,T,F,g,rd,td) 
       else
          stat = oper(T,F,g,rd)
       end if                         !! g =     (WL)' Rd
       stat =  reg(T,T,g,eps*rm)      !! g = g + e  A' Rm

       if(present(kmask)) then
          tm=g
          stat=kmask(F,F,tm,g)
       end if

       stat =  reg(F,F,eps*g,gm)      !! Gm = e A g
       if(present(weight)) then
          call chain(weight,oper,F,F,g,gd,td) 
       else
          stat = oper(F,F,    g,gd)   
       end if                         !! Gd =(WL) g

       stat = step(forget, x,g, rr,gg)
       if(stat ==1) then
          exit !! can't descend anymore
       end if


       if(present(xmov)) xmov(:,iter) =  x(:size(xmov,1))
       if(present(rmov)) rmov(:,iter) = rr(:size(rmov,1))
       if(present(err ))  err(  iter) = dot_product(rd,rd)
       if(v) call solver_report(iter,x,g,rd,rm)

       forget=F
    end do

    if(present(resd)) resd = rd
    if(present(resm)) resm = rm(:size(resm))
  end subroutine solver_reg
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine solver(oper,step, x,dat, niter&
  ,             x0,err,resd,xmov,rmov,verb,weight,kmask)
    optional :: x0,err,resd,xmov,rmov,verb,weight,kmask
    interface
       integer function   oper(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function oper
       integer function weight(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function weight
       integer function  kmask(adj,add, x,dat) result(stat)
         logical, intent(in)  :: adj,add
         real, dimension(:)   :: x, dat
       end function kmask
       integer function step(forget,x,g,rr,gg) result(stat)
         logical              :: forget
         real, dimension(:)   :: x,g, rr,gg
       end function step

    end interface
    !----------------------------------------------------------------
    real, dimension(:),    intent(in)          :: dat, x0
    real, dimension(:),    intent(out)         :: x, err, resd
    integer,               intent(in)          :: niter
    real, dimension(:,:),  intent(out)         :: xmov, rmov   
    logical,               intent(in)          :: verb

    real, dimension(size(x))                   :: g
    real, dimension(size(dat)), target         :: rr, gg
    real, dimension(size(dat)+size(x)), target :: tt
    real, dimension(:), pointer                :: rd, gd, td
    real, dimension(:), pointer                :: rm, gm, tm
    integer                                    :: iter, stat, nd
    logical                                    :: forget, v
    integer                                    :: dprd, dprm, dpx, dpg

    v = F ; if(present(verb)) v = verb
    nd = size(dat)

    rd => rr(1 : nd);
    gd => gg(1 : nd);
    td => tt(1 : nd);  tm=>tt(nd+1:)

    if(present(weight)) then
       stat = weight(F,F,-dat,rd)      !! Rd = - W D
    else
       rd = -dat                       !! Rd = -   D
    end if

    if(present(x0)) then
       stat =   oper(F,F,    x0,td)    !!  t =        L x0

       if(present(weight)) then
          stat = weight(F,T,    td,rd) !! Rd = Rd + W L x0
       else
          rd = rd + td                 !! Rd = Rd +   L x0
       end if
       x = x0 
    else
       x = 0.
    end if

    forget = T; if(present(x0)) forget=F
    do iter = 1,niter

       if(present(weight)) then
          td=rd
          stat = weight(T,F,rd,td) !! Rd=        W' Rd
       end if
       stat = oper(T,F,g,    rd)   !! g =     L' W' Rd

       if(present(kmask)) then
          tm=g
          stat=kmask(F,F,tm,g)
       end if

       stat = oper(F,F,    g,gd)   !! Gd =   L g
       if(present(weight)) then
          td=gd
          stat = weight(F,F,td,gd) !! Gd = W L g
       end if

       stat = step(forget, x,g, rr,gg)
       if(stat ==1) then
          exit !! can't descend anymore
       end if

       if(present(xmov)) xmov(:,iter) =  x(:size(xmov,1))
       if(present(rmov)) rmov(:,iter) = rr(:size(rmov,1))
       if(present(err ))  err(  iter) = dot_product(rd,rd)
       if(v) call solver_report(iter,x,g,rd)

       forget=F
    end do

    if(present(resd)) resd = rd
  end subroutine solver
  !----------------------------------------------------------------
  !----------------------------------------------------------------
  subroutine solver_report(iter,x,g,rd,rm)
    integer,            intent(in) :: iter
    real, dimension(:), intent(in) :: x,g
    real, dimension(:), pointer    :: rd,rm
    optional                       ::    rm
    double precision               :: dd,dm,dx,dg

    dd=dot_product(rd,rd)
    dx=dot_product(x,x)
    dg=dot_product(g,g)

    dm=0.
    if(present(rm)) dm=dot_product(rm,rm)

    if(iter==1) then
       scale=max(dd,dm)/1000.
       write(0,*) "DOT PRODUCT VALUES SCALED BY ",scale
    endif

    write (0,*) "iteration ", iter &
    ,           "   res dat", int(dd/scale) &
    ,           "   res mod", int(dm/scale) &
    ,           "       mod", int(dx/scale) &
    ,           "      grad", int(dg/scale)

  end subroutine solver_report

end module new_solver_mod