module row_mod

  use nchconest1
  use nchconest2
  use ncnhconest1
  use ncnhconest2
 
  implicit none
  logical, parameter, private :: T = .true., F = .false.  
  integer, private            :: ncoeffs1,ncoeffs
  real, private               :: mu1

contains

	subroutine row2_init (ncoeffs1_in, ncoeffs_in, mu1_in)

        !! subroutine arguments
        integer                         :: ncoeffs1_in, ncoeffs_in
        real                            :: mu1_in

        ncoeffs1=ncoeffs1_in         
        ncoeffs=ncoeffs_in         
        mu1=mu1_in

    end subroutine row2_init

    function row2_lop (adj, add, mod, dat) result(stat)

         !! subroutine arguments
         integer                       :: stat
         logical, intent (in)          :: adj, add
         real, dimension (:)           :: mod, dat

         !! call row2 operator to compute M+muP
         call row2(nchconest1_lop,nchconest2_lop,adj,add,mod(1:ncoeffs1), &
              & mod(ncoeffs1+1:ncoeffs),dat,mu1) 
       
         stat=1
		
    end function row2_lop 

    !! same as row2_lop but for non-stationary filters
    function nrow2_lop (adj, add, mod, dat) result(stat)

         !! function arguments
         integer                       :: stat
         logical, intent (in)          :: adj, add
         real, dimension (:)           :: mod, dat

         !! call row2 operator to compute M+muP
         call row2(ncnhconest1_lop,ncnhconest2_lop,adj,add,mod(1:ncoeffs1), &
              & mod(ncoeffs1+1:ncoeffs),dat,mu1) 
		
         stat=1

    end function nrow2_lop 


  subroutine row2 (oper1, oper2, adj, add, mod1, mod2, dat, mu)
    interface
       function oper1 (adj, add, mod, dat) result (stat)
         integer              :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper1
       function oper2 (adj, add, mod, dat) result (stat)
         integer    :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper2
    end interface
    logical, intent (in)  :: adj, add
    real, dimension (:)   :: mod1, mod2, dat
    real                  :: mu 

    real, dimension(:), allocatable :: tmp,tmp1

    integer               :: stat1, stat2

    allocate(tmp(size(dat)))
    allocate(tmp1(size(mod1)))

    if (adj) then


		tmp1=0.

       stat1 = oper1 (T, F, tmp1, dat)

		tmp1=mu*tmp1

		if (add) then
       	mod1  = mod1 + tmp1
		else
			mod1=tmp1
		end if

       stat2 = oper2 (T, add, mod2, dat)

    else

		tmp = 0.

       stat1 = oper1 (F, F, mod1, tmp)

		tmp=mu*tmp

       if (add) then 
		dat = tmp+dat
       else
         dat = tmp
       end if

       stat2 = oper2 (F,   T, mod2, dat)

    end if

    deallocate(tmp,tmp1)

  end subroutine row2


  subroutine row3 (oper1, oper2, oper3, adj, add, mod, &
		dat, eps1, eps2, eps3, n1,n2,n3)
    interface
       function oper1 (adj, add, mod, dat) result (stat)
         integer              :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper1
       function oper2 (adj, add, mod, dat) result (stat)
         integer    :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper2
       function oper3 (adj, add, mod, dat) result (stat)
         integer    :: stat
         logical, intent (in) :: adj, add
         real, dimension (:)  :: mod, dat
       end function oper3
    end interface

    logical, intent (in)  :: adj, add
    integer, intent (in)  :: n1,n2,n3
    real, dimension (:)   :: mod
	real, dimension (:)	  :: dat
    real                  :: eps1, eps2, eps3

    real, dimension(:), allocatable :: tmp_mod1,tmp_mod2,tmp_mod3
    real, dimension(:), allocatable :: tmp_dat1,tmp_dat2,tmp_dat3

    integer               :: stat1, stat2, stat3

    allocate(tmp_mod1(size(mod)))
    allocate(tmp_mod2(size(mod)))
    allocate(tmp_mod3(size(mod)))
    allocate(tmp_dat1(n1))
    allocate(tmp_dat2(n2))
    allocate(tmp_dat3(n3))

    if (adj) then
	
		tmp_mod1=0.	
		tmp_mod2=0.	
		tmp_mod3=0.	

       stat1 = oper1 (T, F, tmp_mod1, dat(1:n1))
       stat2 = oper2 (T, F, tmp_mod2, dat(n1+1:n1+n2))
       stat3 = oper3 (T, F, tmp_mod3, dat(n1+n2+1:n1+n2+n3))

		if (add) then
			mod=mod+eps1*tmp_mod1+mu1*eps2*tmp_mod2+eps3*tmp_mod3
		else
			mod=eps1*tmp_mod1+mu1*eps2*tmp_mod2+eps3*tmp_mod3
		end if

    else

       tmp_dat1=0.

       stat1 = oper1 (F, F, mod, tmp_dat1)

       tmp_dat1=tmp_dat1*eps1 

       if (add) then
          dat(1:n1) = dat(1:n1)+tmp_dat1
       else
          dat(1:n1) = tmp_dat1
       end if

       tmp_dat2 = 0.

       stat1 = oper2 (F, F, mod, tmp_dat2)

       tmp_dat2=mu1*tmp_dat2*eps2 

       if (add) then
          dat(n1+1:n1+n2) = dat(n1+1:n1+n2)+tmp_dat2
       else
          dat(n1+1:n1+n2) = tmp_dat2
       end if

       tmp_dat3 = 0.

       stat3 = oper3 (F, F, mod, tmp_dat3)

       tmp_dat3=tmp_dat3*eps3 

       if (add) then
          dat(n1+n2+1:n1+n2+n3) = dat(n1+n2+1:n1+n2+n3)+tmp_dat3
       else
          dat(n1+n2+1:n1+n2+n3) = tmp_dat3
       end if

    end if

    deallocate(tmp_mod1,tmp_mod2,tmp_mod3)
    deallocate(tmp_dat1,tmp_dat2,tmp_dat3)

  end subroutine row3

end module row_mod
