module solver
use precision
implicit none

integer,private :: bw !! half-bandwidth, tridiag has bw=1

private :: computebw

contains

subroutine gauss_lu(A)
! performs an LU factorization of the matrix A using
! Gaussian reduction.
! A is the matrix to be factored.
! L is a lower triangular factor with 1's on the diagonal
! U is an upper triangular factor.
! A = L * U
! L is stored as the lower triangle without diagonal of A, U is
!      stored as the upper triangle of A
! from the Matlab
! will abort if pivoting is needed

! mms 3/4/13
real(PREC) :: A(:,:)
integer :: n,Irow,Jcol,i,i4max(1),colup
real(PREC) :: mult, maxPivotRatio=0.0_PREC
logical  :: pivot=.false.

call computebw(A)
!print *,"bw=",bw

n=size(A,1);

do Jcol=1,n
  ! check if pivoting needed.  Hope not!
  i4max = maxloc ( abs ( A(Jcol:min(Jcol+bw,n), Jcol ) ) );
  if (i4max(1) /= 1) then
    pivot=.true.
    maxPivotRatio=max(maxPivotRatio,abs(A(Jcol-1+i4max(1),Jcol)/A(Jcol,Jcol)))
    if (abs(A(Jcol-1+i4max(1),Jcol)) > abs(A(Jcol,Jcol))*10._PREC) &
      & stop 'gauss_lu: Pivoting needed!'
  endif
  do Irow=Jcol+1,min(Jcol+bw,n)
    ! compute Irow multiplier and save in L(Irow,Jcol)
    mult=A(Irow,Jcol)/A(Jcol,Jcol);
  
    ! multiply row "Jcol" by L(Irow,Jcol) and subtract from row "Irow"
    A(Irow,Jcol:min(Jcol+bw,n))=A(Irow,Jcol:min(Jcol+bw,n)) &
      &                         -mult*A(Jcol,Jcol:min(Jcol+bw,n));

    ! Now A(Irow,Jcol) should be zero
    if(abs(A(Irow,Jcol))>abs(mult)*epsilon(mult)*1.0e5_PREC)then
      print *,'Irow=',Irow,' Jcol=',Jcol
      print *,'A(Irow,Jcol)=',A(Irow,Jcol),' was=',mult
      stop 'gauss_lu: did not get zero below diagonal'
    endif
    ! put multiplier there
    A(Irow,Jcol)=mult;
  enddo
enddo
if(pivot) print '(a,f8.4)',"Maximum ratio (pivot)/(diag)=",maxPivotRatio

end subroutine gauss_lu

subroutine l_solve(A,z)
  ! L is the lower-triangular piece of A, with diagonal assumed=1
  ! z is the right hand side
  ! the solution of L*y=z is returned in z
  ! from the Matlab

  real(PREC),intent(in)    :: A(:,:)
  real(PREC),intent(inout) :: z(:)
  integer                  :: Irow,Jcol
  real(PREC)               :: y(size(z))
  
  ! first row is really an easy one, especially since the diagonal
  ! entries of L are equal to 1
  Irow=1;
  
  do Irow=2,size(z)
    do Jcol=max(1,Irow-bw),Irow-1
      z(Irow) = z(Irow) - A(Irow,Jcol)*z(Jcol);
    enddo
  enddo
  
end subroutine l_solve

subroutine u_solve(A,y)
  ! U is the upper-triangular part of matrix A
  ! y is the right hand side
  ! x is the solution of U*x=y
  ! from the Matlab
  
  real(PREC),intent(in)    :: A(:,:)
  real(PREC),intent(inout) :: y(:)
  integer                  :: Irow,Jcol
  real(PREC)               :: x(size(y))
  
  ! last row is the easy one
  Irow=size(y);
  y(Irow)=y(Irow)/A(Irow,Irow);
  
  do Irow=size(y)-1,1,-1
    do Jcol=Irow+1,min(Irow+bw,size(y))
      y(Irow) = y(Irow) - A(Irow,Jcol)*y(Jcol);
    enddo
    y(Irow) = y(Irow)/A(Irow,Irow);
  enddo
end subroutine u_solve

subroutine computebw(A)
real(PREC),intent(in) :: A(:,:)
integer :: i,j,ubw,lbw

bw=0
do i=1,size(A,1)
  ! upper half-band for this row
  ubw=0
  do j=i+1,size(A,2)
    if(A(i,j) /= 0.0_PREC)ubw=j-i
  enddo

  ! lower half-band for this row
  lbw=0
  do j=1,i-1
    if(A(i,j) /= 0.0_PREC)lbw=i-j
  enddo
  bw=max(bw,lbw,ubw)
enddo

end subroutine computebw

end module solver