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
Comments