lutest.f90

module auxsubs

contains

  subroutine matinit(a, b)
    ! initialize matrix a and vector b with random numbers
    
    use kinds
    implicit none    
    real(kind=REAL8), dimension(:,:) :: a
    real(kind=REAL8), dimension(:)   :: b

    call random_number(a)
    call random_number(b)
   
    return
  end subroutine matinit

  subroutine check(a, b, x)
    ! compute residual ||a*x  - b||/( ||x|| * ||a|| * eps * n )
    
    use kinds
    implicit none    
    real(kind=REAL8), dimension(:,:) :: a
    real(kind=REAL8), dimension(:)   :: b, x

    real(kind=REAL8), dimension(size(a,1))  :: work
    real(kind=REAL8)                 :: anorm, xnorm, eps, resid, resnorm
    real(kind=REAL8)                 :: dlamch, dlange
    integer                          :: n, info

    n = size(a,1)

    print *, 'lapack example program'
    print *, 'solving ax=b via lu decomposition'
    print *, '   a             :', n,  ' by ', n
    print *, 'info code returned by dgesv = ', info
    
    eps = dlamch('epsilon')                  ! get machine precision
    anorm = dlange('i', n, n, a, n, work)    ! get infinity-norm of a and x
    xnorm = dlange('i', n, 1, x, n, work)
    
    ! compute a*x - b and its norm
    ! dgemm replaces b by the result
    call dgemm('n', 'n', n, 1, n, 1.0_REAL8, a, n, x, n, -1.0_REAL8, b, n)
    resnorm = dlange('i', n, 1, b, n, work)
    
    ! check residual error
    resid = resnorm / (anorm*xnorm*eps*n)
    if( resid < 10.0 ) then
       print *, 'solution is correct! residual:'
    else
       print *, 'solution is incorrect! residual:'
    end if
    print *, '   ||a*x - b|| / ( ||x||*||a||*eps*n ) = ', resid
    
    return
  end subroutine check
end module auxsubs


program lutest

  ! example program solving ax=b via lapack routine dgesv
  ! set number of threads with MLIB_NUMBER_OF_THREADS
  
  use kinds
  use auxsubs
  implicit none
  
  real(kind=REAL8), dimension(:,:), allocatable  :: a, adecomp
  real(kind=REAL8), dimension(:), allocatable    :: b, x
  integer, dimension(:), allocatable             :: ipiv
  integer                                        :: n, info

  ! get size of problem
  print*, 'enter number of equations:'
  read*, n
  allocate(a(n,n))
  allocate(adecomp(n,n))
  allocate(b(n))
  allocate(x(n))
  allocate(ipiv(n))

  ! initialize matrix and vector
  call matinit(a, b)

  ! Call the lapack routine dgesv to solve the linear system a * x = b.
  ! dgesv replaces a by its lu decomposition and b by the solution.
  ! Since we need the original data for checking, we make copies.
  adecomp = a
  x = b
  call dgesv(n, 1, adecomp, n, ipiv, x, n, info)

  ! check the result
  call check(a, b, x)
  
end program lutest