! Inversa matr. HILBERT cu utilizare rutine IMSL
!
!  Pentru matricea HILBERT: Se pune: a(i,j)=1./a(i,j)
!
!Usage

!CALL LSARG (N, A, LDA, B, IPATH, X)

!Arguments

!N = Number of equations.   (Input)

!A = N by N matrix containing the coefficients of the linear system.   (Input)

!LDA = Leading dimension of A exactly as specified in the dimension statement of the calling program.   (Input)
!B = Vector of length N containing the right-hand side of the linear system.   (Input)
!IPATH = Path indicator.   (Input) 
!IPATH = 1 means the system AX = B is solved. 
!IPATH = 2 means the system ATX = B is solved.

!X = Vector of length N containing the solution to the linear system.   (Output)
use numerical_libraries

character*128 fname, charn*2
real, allocatable :: A(:,:), b(:,:), x(:,:)

print '(a, \)', 'Inversa matricii HILBERT de ordinul:'
read *, n
print '(a, \)', 'Iterative refinment (0/1):'
read *, kodRefine

if(n <=9) then
in =1
else
in =2
endif
write(charn, '(I<in>)') n

if(kodRefine ==0) then
fname ='InvHilbert'//trim(charn)//'.rez'
else
fname ='InvHilbert'//trim(charn)//'-refin'//'.rez'
endif
open(1, file =fname, status ='unknown')


allocate (A(n, n), b(n,n), x(n,n))

! Generate Hilbert matrix:
ido: do i =1,n
jdo: do j =1,n
a(i,j) =1./(i+j-1)
enddo jdo
enddo ido

! Genereaza termenii liberi:
do1: do i =1, n
do2: do j=1, n
if (i==j) then
b(i, j) =1.
else
b(i, j) =0.
endif
enddo do2
enddo do1


ndim =n
ipath =1
!	call ELIM(AB, N, NP, NDIM, kod)

solve: DO nb =1,n
if(kodRefine ==0) then
CALL LSLRG (n, a, ndim, b(:, nb), IPATH, x(:,nb))	! without iterative refinement
else
CALL LSARG (n, a, ndim, b(:, nb), IPATH, x(:,nb))	! with iterative refinement
endif
ENDDO solve
	
! Condition number:


! Write results
write(1,'(a, I3)') 'Inversa (IMSL) matricii HILBERT de ordinul', n
if(kodRefine ==0) then
write(1,'(a, /)')	'Fara rafinare iterativa' 
else
write(1,'(a, /)')	'Cu rafinare iterativa'
endif

do i =1, n
write(1,*) (x(i, j), j =1,n)
enddo

deallocate(A, b, x)
END


