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