I am trying to solve problems of the form Ax = b
by using the function sgesvx
. However, it is returning a wrong solution for one test case: - 0.33333 * x_1 - 0.06667 * x_2 = - 0.133333
, which is returning:
x_1 = 0.01346
; Andx_2 = 0
.
That, when applied over the expression - 0.33333 * x_1 - 0.06667 * x_2 = - 0.133333
, gives - 0.0044866218 = - 0.133333
, which is invalid. Below we have the code.
program Main
implicit none
external :: sgesv, sgesvx
! SCALAR
integer, parameter :: n=2,m=1
integer :: i,j
! ARRAY
real(kind=8) :: A(m,n),b(m),x(n)
character(len=15) :: subname
! LAPACK
integer :: IPIV(n),IWORK(n),INFO
character :: EQUED='N'
real(kind=8) :: R(n),C(n),RCOND,FERR(1),BERR(1),WORK(4*n),AF(m,n)
!
x(1) = 0
x(2) = 0
A(1,1) = -0.33333
A(1,2) = -0.06667
b(1) = -0.133333
!
do i = 1,m
print *,"Constraint ",i,": "
do j = 1,n
print "(f10.5)",A(i,j)
end do
print *," = ",b(i)
end do
!
call sgesvx( &
'E', &
'N', &
n, &
1, &
A, &
n, &
AF, &
n, &
IPIV, &
EQUED, &
R, &
C, &
b, &
n, &
x, &
n, &
RCOND, &
FERR, &
BERR, &
WORK, &
IWORK, &
INFO &
)
if (INFO .lt. 0) then
print *,INFO,"-th argument of sgesv is invalid"
else if (INFO .gt. 0) then
print *,INFO
end if
do j = 1,n
print "(f10.5)",x(j)
end do
end program Main
For compiling, simply type gfortran -O3 linearsystem.f90 -o main -lblas -llapack
.
13