module new_big_solver_mod
  use chain_mod
  use llist_mod
  use arnoldi_mod
	use ddot_mod
  implicit none
  logical, parameter, private  :: T = .true., F = .false.
  !----------------------------------------------------------------
contains

  subroutine solver_prec(oper,step, prec,nprec, x,dat, niter,eps &
  ,             p0,rm0,err,resd,resm,xmov,rmov,verb,weight,kmask)
    optional :: p0,err,resd,resm,xmov,rmov,verb,weight,kmask,rm0
    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
		real                                        :: mys
		double precision                            :: dd,dm,dx,dg

    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 = 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
			 if(v) then
					dd=dot_product(rd,rd)
					dm=dot_product(rm,rm)
					dx=dot_product(x,x)
					dg=dot_product(g,g)
					if(iter==1) then
						mys=max(dd,dm)/1000.
            write(0,*) "DOT PRODUCT VALUES SCALED BY ",mys
					endif
					  write (0,*) "iteration ", iter,  &
          " res dat", int(dd/mys), " res mod", int(dm/mys),  &
          " mod ", int(dx/mys), " grad", int(dg/mys)
			 end if

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

!			call srite("gmov.H",g,size(g)*4)
!			call srite("pmov.H",p,size(p)*4)
       stat = prec(F,F,p,x)         !! x = P p
!			call srite("rmov.H",rr,size(rr)*4)

       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)
       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,err,resd,resm,xmov,rmov,verb,weight,kmask,rm0
    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
		real                                        :: mys
		double precision                            :: dd,dm,dx,dg

    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)   !! g =     L' W' Rd
			 end if
			 
       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)   !! Gd =   L g
       end if

			 if(v) then
					dd=dot_product(rd,rd)
					dm=dot_product(rm,rm)
					dx=dot_product(x,x)
					dg=dot_product(g,g)
					if(iter==1) then
						mys=max(dd,dm)/1000.
            write(0,*) "DOT PRODUCT VALUES SCALED BY ",mys
					endif
					  write (0,*) "iteration ", iter,  &
          " res dat", int(dd/mys), " res mod", int(dm/mys),  &
          " mod ", int(dx/mys), " grad", int(dg/mys)
			 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)
       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
		real                                        :: mys
		double precision                            :: dd,dm,dx,dg

    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

			 if(v) then
					dd=dot_product(rd,rd)
					dx=dot_product(x,x)
					dg=dot_product(g,g)
					if(iter==1) then
						mys=max(dd,dm)/1000.
            write(0,*) "DOT PRODUCT VALUES SCALED BY ",mys
					endif
					  write (0,*) "iteration ", iter,  &
          " res dat", int(dd/mys),   &
          " mod ", int(dx/mys), " grad", int(dg/mys)
			 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)
       forget=F
    end do

    if(present(resd)) resd = rd
  end subroutine solver

end module new_big_solver_mod
