program solve9
use precision
use fem,only: makemat,makerhs,Mat,rhs,mesh
use solver, only: gauss_lu,l_solve,u_solve
implicit none

integer :: Ntest=7
integer :: i,j
real(PREC),allocatable :: x(:),A(:,:)
real(PREC) :: err
character(len=100) :: line

! CALL GET_COMMAND_ARGUMENT(NUMBER [, VALUE, LENGTH, STATUS]) 
call get_command_argument(number=1,value=line,length=j,status=i)
if (i < 0)stop 'problem reading command line'
if (j>=100)stop 'command line buffer too short'
if (len_trim(line) > 0) then
  if (len_trim(line) /= j) stop 'bad result from get_command_argument'
  read(line,*)Ntest  !! this format might be illegal
endif
print *,'Using Ntest=',Ntest

allocate(x(Ntest),A(Ntest,Ntest))

call makemat(Ntest)
A=Mat

call makerhs(frhs9)
x=rhs

! factor A, put factors inside A
call gauss_lu(A)
! L-solve, put back into x
call l_solve(A,x)
! U-solve, put back into x
call u_solve(A,x)

!print '(10e12.4)',mesh
!print '(10e12.4)',x
!print '(10e12.4)',exact1(mesh)

err=norm(exact9(mesh)-x)/norm(exact9(mesh))

print *,'N=',Ntest,' err=',err

contains
function frhs1(k,x)
  real(PREC),intent(in) :: x
  real(PREC)            :: frhs1
  integer,intent(in)    :: k
! y''+y'+y=rhs
! y  = x*(1-x) = x-x^2
! y' = 1-2*x
! y''=-2
! rhs=-1-x-x^2
  frhs1=-1.0_PREC-1.0_PREC*x-x**2
end function frhs1

function exact1(x)
  real(PREC),intent(in) :: x(:)
  real(PREC)            :: exact1(size(x))
! y''+y'+y=rhs
! y  = x*(1-x) = x-x^2
! y' = 1-2*x
! y''=-2
! rhs=-1-x-x^2
  exact1=x*(1.0_PREC-x)
end function exact1

function frhs9(k,x)
  real(PREC),intent(in) :: x
  real(PREC)            :: frhs9
  integer,intent(in)    :: k
  frhs9=1.0_PREC+x
end function frhs9

function exact9(x)
  real(PREC),intent(in) :: x(:)
  real(PREC)            :: exact9(size(x))
  real(PREC)            :: omega,one=1.0_PREC,two=2.0_PREC,three=3.0_PREC
  omega=sqrt(three)/two;
  exact9=x-exp(-(x-one)/two)*sin(omega*x)/sin(omega)
end function exact9

function norm(x)
  real(PREC) :: x(:),norm
  norm = sqrt(sum(x*x))
end function norm


end program solve9