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; And
• x_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.