! Inversa exacta a matricii Hilbert de ordinul "n"
! + Numarul de conditie Cond1(A)

real, allocatable :: Hinv(:,:)
character*20 fname, charn*2
real numer
real norm1

print '(a, \)', 'Inversa matricii HILBERT de ordinul:'
read *, n

allocate(Hinv(n,n))
if(n <=9) then
in =1
else
in =2
endif
write(charn, '(I<in>)') n
fname ='InvHilbert'//trim(charn)//'.rez'
open(1, file =fname, status ='unknown')

ido: do i=1, n
fi_1 =fact(i-1)
fn_i =fact(n-i)
fni_1=fact(n+i-1)

jdo: do j=1, n
if(mod(i+j,2) ==0) then
dij =1
else
dij =-1
endif

fnj_1 =fact(n+j-1)
fj_1 =fact(j-1)
fn_j =fact(n-j)	

numer =dij*fni_1*fnj_1
denom =(i+j-1)*(fi_1*fj_1)*(fi_1*fj_1)*fn_i*fn_j

Hinv(i,j) = numer/denom 
enddo jdo
enddo ido

! Write results
write(1,'(a, I3, /)') 'Inversa (EXACTA) matricii HILBERT de ordinul', n

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

! Numarul de conditie:
! Sumeaza prima linie din H = norma1(H)
sum =1.
do j =2,n
sum =sum +1./j
enddo

cond1 =norm1(Hinv,n)*sum

write(1,'(/, a, 1pE10.2)') 'Numarul de conditie: Cond1 =', cond1

deallocate(Hinv)
END

! Factorial
recursive function fact(n)	!integer

fact =1
if(n==0) return

do i =1, n
fact =fact*i
enddo
END

Real function norm1(A,n)
! 1-norma matricii A(n*n)= norma liniilor
real a(n,n), aline(n)

do i=1,n
aline =SUM(abs(A), dim =2)
enddo

norm1 =maxval(aline)

END function norm1
