!  Program principal pt. ELIM+
!  Matricea HILBERT: generata
! - Inversa
! + Inversa Exacta
! + Eroarea
! + Determinantul

character*128 fname, charn*2
real, allocatable :: A(:,:), B(:,:)
integer width

real, allocatable :: A0(:,:), C(:,:)
logical do_check

real, allocatable :: Hinv(:,:), Rel_Error(:,:)
real Norm1
integer result(2)

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

print '(a, \)', ' Verificare (T/F): '
read(*,'(L)') do_check

iwn =width(n)
write(charn, '(I<iwn>)') n

fname ='Inv-H'//trim(charn)//'.rez'
open(1, file =fname, status ='unknown')

allocate (A(n, n), B(n,n))

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

if(do_check) then
allocate (A0(n, n), C(n,n))
A0 =A; C=0.;
endif

! Genereaza termenii liberi:
b =0.
do i =1, n
b(i,i) =1.
enddo

nb =n
prag =1E-10

!call ELIM2 (n, nb, A, B, prag, kod)
kodPrint =0
call ELIM (n, nb, A, B, prag, kod, kodPrint, n_sch_linii)

! Rezultate:
!iwn =width(n)
write(1,'(a, I<iwn>)') 'Inversa GAUSS a matricii HILBERT de ordinul ', n

if (kod <0) then
write(1,'(/, a)') '** Pivot < prag: Matrice aproape singulara.'
goto 500
endif

nn =min(8,n)

do i =1, n
write(1, '(<nn>(1pg16.7))') b(i,1:n)
enddo


!Inversa Exacta:
allocate(Hinv(n,n))
call HInverse_Exact(n, Hinv)

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

do i =1, n
!write(1,*) (Hinv(i,j), j =1,n)
write(1, '(<nn>(1pg16.7))') Hinv(i,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

! Eroarea Relativa maxima:
allocate(Rel_Error(n,n))
Rel_Error =0

rel_err_max =0

do i=1,n
do j=1,n
Rel_Error(i,j) =(Hinv(i,j) -b(i,j))/Hinv(i,j)
Rel_Error_1 =Rel_Error(i,j)

if(abs(Rel_Error_1) > rel_err_max) then
rel_err_max =Rel_Error_1
endif

enddo
enddo

!rel_err_max =maxval(abs(Rel_Error))
!result =maxloc(abs(Rel_Error))
write(*,*) rel_err_max
write(1,'(/, a, 1pE10.3)') 'Eroarea Relativa maxima (in modul) =', rel_err_max

do i=1,n
do j=1,n
Rel_Error_1 =Rel_Error(i,j)

if( abs(Rel_Error_1) == abs(rel_err_max) ) then
write(1,'(a, I1,a,I1, a)') 'Element: (', i,',',j, ')'
write(*,'(a, I1,a,I1, a)') 'Element: (', i,',',j, ')'
endif

enddo
enddo

! Determinant-Gauss:
det =1.

do i=1,n
det =det*a(i,i)
enddo

det =det*(-1)**n_sch_linii

write(1,'(/, a, \)') 'Determinant-Gauss: '
write(1,*) det

! Proba:
if(do_check) then
C =matmul(A0,B)

write(1, '(/, a)') ' Proba (inversa Gauss):'
do i =1, n
write(1, '(<nn>(1pg16.7))') C(i, 1:n)
enddo

endif

500 deallocate(A, B)
if(allocated(A0)) deallocate(A0)
if(allocated(C)) deallocate(C)

deallocate(Hinv)        !, Rel_Error)

write(*,'(/, a, a )') '* Fisier de iesire: ', trim(fname)

END

