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