top of page
  • Writer's pictureAdisorn O.

LAPACK Package Testing for LU Solver Ax = b

To solve the linear system Ax = b, LAPACK provide a function dgesv() see more detail

dgesv employs fast LU decomposition which is one of the most efficient and robust direct solver for a linear system.

F90 code:

program solve_linear_system

implicit none

integer :: n, lda, ldb, nrhs, info

integer, dimension(:), allocatable :: ipiv

real(8), dimension(:,:), allocatable :: a, b

! Define the order of the matrix A

n = 3

nrhs = 1

lda = n

ldb = n

! Allocate arrays

allocate(a(n, n), b(n, nrhs), ipiv(n))

! Initialize matrix A (note that A is stored by column so that it's equivalent to A' in Matlab)

a = reshape([3.0d0, 1.0d0, 2.0d0, &

6.0d0, 3.0d0, 4.0d0, &

3.0d0, 1.0d0, 5.0d0], shape(a))

! Initialize right-hand side B

b = reshape([1.0d0, 2.0d0, 3.0d0], shape(b))

! Print the matrix A and vector b before solving

print *, "Matrix A:"

print *, a

print *, "Vector b:"

print *, b

! Call LAPACK routine to solve the system A*x = b

call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

! Check if the solution is successful

if (info == 0) then

print *, "The solution is:"

print *, b


print *, "The solution could not be computed. Info =", info

end if

! Deallocate arrays

deallocate(a, b, ipiv)

end program solve_linear_system


Matrix A:

3.0000000000000000 1.0000000000000000 2.0000000000000000

6.0000000000000000 3.0000000000000000 4.0000000000000000

3.0000000000000000 1.0000000000000000 5.0000000000000000

Vector b:

1.0000000000000000 2.0000000000000000 3.0000000000000000

The solution is:

-3.7777777777777781 1.6666666666666667 0.77777777777777779

For this test example, I use BLAS/ LAPACK package installed from Simply Fortran 3.0 via the package manager. The package library must be linked to the project before building.

It should be noted that the matrix A is stored by column in F90 (against by row in MATLAB). So [A] is equal to [A]' in MATLAB



bottom of page