program sgefs_main

implicit none
integer, parameter :: n=3
integer :: i, j, itask, ind
integer, dimension(n) ::   iwork
real :: rcond
real, dimension(n) :: b, work
real, dimension(n,n) :: A

! Set up problem:
itask  =  1
A = reshape( (/10.0,-3.0,5.0,-7.0,2.0,-1.0,0.0,6.0,5.0/),(/3,3/) )
b = (/7.0,4.0,6.0/)

! Print problem data:
write(*,*) 'Augmented matrix [A,b]:'
do i=1,n
   write(*,10) (A(i,j), j = 1,n), b(i)
10 format(4f12.6)
enddo

! Solve linear system:
call sgefs(A,n,n,b,itask,ind,work,iwork,rcond)

! Solution:
write(*,*)
write(*,*) 'Sgefs results:'
if (ind==-10) then
   write(*,*) 'Error code =', ind
elseif (ind<0) then
   write(*,*) 'Error code =', ind
   stop
else
   write(*,*) 'Number of accurate digits =', ind
endif

write(*,*) 'Solution:'
write(*,20) 'x = ',(b(j), j = 1,n)
20 format(a,4x,3f12.6)

stop
end program sgefs_main