# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
!Sergey's
module divide
  use weight
  use solver_mod
  use cgstep_mod
  use tridiag
    implicit none
  interface divide_apply
    module procedure divide_apply_cg
    module procedure divide_apply_tri
  end interface
  contains
  subroutine divide_apply_tri(n1,n2,num,den,rat,eps1,eps2)
    integer :: n1,n2
    real :: num(:),den(:),eps1,rat(:),eps2
    call divide_apply_tri2(n1,n2,num,den,rat,eps1,eps2)
  end subroutine
  subroutine divide_apply_tri2(n1,n2,num,den,rat,eps1,eps2)
    integer :: i2,n1,n2
    real :: num(n1,n2),eps1
    real :: den(n1,n2),eps2,rat(n1,n2)
    real, allocatable :: d(:),a(:),b(:),c(:)
    integer :: ifirst,ilast
    allocate(d(n1),a(n1),b(n1),c(n1))
!first solve
!
! d = L m ---> num*den = den*den num/den
! 0 = e A m ---> (-1 2 -1) in 1-D
!
    ifirst=1
    do while (sum(abs(den(:,ifirst)))<.001 .and. ifirst < n2-1)
      ifirst= ifirst+1
    end do
    ilast=n2
    do while (sum(abs(den(:,ilast)))<.001 .and. ilast > 1)
      ilast= ilast-1
    end do
    rat=0
    if (ifirst <ilast) then
      call rtris_init(n1)
      d=num(:,ifirst)*den(:,ifirst)
      b=den(:,ifirst)*den(:,ifirst)+ 2 *eps1
      a=-eps1
      c=a
      b(1)=b(1)-eps1
      b(n1)= b(n1)-eps1 !boundary condition
      call rtris_back(a,b,c,d,rat(:,ifirst))
      do i2=ifirst+1,ilast
!NEXT add Extra reg
!
! m(:,i2-1) * eps2 = eps2
!
        d=num(:,i2)*den(:,i2)+eps2*rat(:,i2-1)
        b=den(:,i2)*den(:,i2)+ 2 *eps1 + eps2
        b(1)= b(1)-eps1
        b(n1)= b(n1)-eps1 !boundary condition
        call rtris_back(a,b,c,d,rat(:,i2))
      end do
      call rtris_close()
    end if
  end subroutine
  subroutine divide_apply_cg (doprec, nd, num, den, rat, eps, niter,&
    & prec, nprec,res)
    logical, intent (in) :: doprec
    integer, intent (in) :: nd, nprec, niter
    real, intent (in) :: eps
    real, dimension (nd), intent (in) :: num, den
    real, dimension (nd), intent (out) :: rat,res
    interface
      integer function prec (adj, add, x, dat)
        logical, intent (in) :: adj, add
        real, dimension (:) :: x, dat
      end function
    end interface
    real, dimension (:), pointer :: wt
    allocate (wt (nd))
    wt = den
    call weight_init (wt)
    if (doprec) then
      call solver_prec (weight_lop, cgstep, prec, nprec,x = rat, dat =&
        & num, niter = niter, eps = eps,res=res)
    else
      call solver_reg (weight_lop, cgstep, prec, nprec,x = rat, dat =&
        & num, niter = niter, eps = eps,res=res)
    end if
    call cgstep_close ()
    deallocate (wt)
  end subroutine
end module
