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 https://netlib.org/lapack/explore-html-3.6.1/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html


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

else

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

end if


! Deallocate arrays

deallocate(a, b, ipiv)

end program solve_linear_system


Output:


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




3 views

Comentarios


bottom of page