program testblas3
  use RealKind
  use blocklu
  ! compares times for lu1blas2 and lu1blas3 on random well-conditioned
  ! test matrices
  integer :: n
  real (kind=rk),dimension(:,:),allocatable :: A,A1
  real (kind=rk),dimension(:),allocatable :: b,x,x_e
  real :: t1, t2
  integer :: isolve
  integer :: bsize,ierror
  print*,'type n'
  read*, n
  print*, 'value of n =',  n
  ! allocate storage
  allocate(A(n,n))
  allocate(A1(n,n))
  allocate(x(n))
  allocate(x_e(n))
  allocate(b(n))
  ! create a random well-conditioned symmetric positive definite matrix
  call rsymm(A,n)
  A1 = A                     ! save as A1 (only needed for the error check)
  ! do i = 1,n
  ! print*,'matrix A'
  ! print*,A(i,:)
  ! end do 
  !  random solution vector
  call random_number(x_e)
  ! assemble RHS: so that Ax = b has the exact solution x_e
  b  = matmul(A,x_e)
  ! now solve
  isolve = 0
  do
    if (isolve==0) then
      A = A1             ! reload A
      print*,'which solver?  2  for lu1blas2, 3 for lu1blas3 '
      read*,isolve
      call cpu_time(t1)    ! start the clock
      if (isolve == 2) then
        print*,'using solver blu1blas2'
        call blu1blas2(A,n,n,1)
      elseif (isolve == 3) then
        print*,'using solver lu1blas3'
        print*,'type block size'
        read*,bsize
        call lu1blas3(A,n,bsize,ierror)
      endif
      call cpu_time(t2)    ! stop the clock
      call bs1(A,b,n,x)     ! backsolve
      print*, 'error in solve'   ! check that the solver has worked
      print*, sqrt(dot_product(x-x_e,x-x_e))
      print*, 'Time taken for solver', t2 - t1, 'seconds   '
                                 ! print the time 
      print*,'Type 0 for another solve, 1 to stop'
      read*, isolve
    else
      stop
    endif
  enddo
  contains
    subroutine rsymm(A,n)
    ! assembles a random well-conditioned n times n symmetric matrix A
    ! It is just a small symmetric random perturbation of the identity
    ! and so should be well-conditioned
      integer, intent(in) :: n
      real (kind=rk), dimension(n,n), intent(out) :: A
      integer :: i
      call random_number(A)
      A = 0.01_rk*(A*transpose(A))    ! small symmetric matrix
      do i = 1,n
        A(i,i) = A(i,i) + 1.0_rk
      end do
    end subroutine rsymm
end program testblas3