! Matricea HILBERT/LOTKIN: generata
! Inversa Gauss
! + Inversa Exacta
! + Eroarea
! + Numarul de conditie Cond(A)1 (cu inversa exacta).

character*128 fname, charn*2
real, allocatable :: A(:,:), B(:,:) ! A =matricea Hilbert/Lotkin
integer width

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

real, allocatable :: Hinv(:,:), Rel_Error(:,:)
real Norm1
integer result(2)
logical is_Hilbert; character(1) w1;
character(7) matrix_name


print '(a, \)', ' Matricea HILBERT (H,h; Y,y) sau LOTKIN (char #H,h; Y,y): '
read(*,'(a1)') w1

is_Hilbert =(w1 =='H' .or. w1 =='h' .or. w1 =='Y' .or. w1 =='y')

if(is_Hilbert) then
matrix_name ='HILBERT'
else
matrix_name ='LOTKIN '
endif

write(*,*) matrix_name


print '(a, \)', ' Inversa matricii '//trim(matrix_name)//' de ordinul: '
read *, n

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


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

if(is_Hilbert) then
fname ='Inv-H'//trim(charn)//'.rez'
else
fname ='Inv-L'//trim(charn)//'.rez'
endif

open(1, file =fname, status ='unknown')
! Title:
iwn =width(n)
write(1,'(a,I<iwn>))') 'Matricea '//trim(matrix_name)//' de ordinul: ', n
write(1,'(a,/)')  '------------------------------'

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(.not. is_Hilbert) then
a(1,:) =1.
endif

! Pentru Numar de conditie:
ANorm1 =Norm1(A,n)
 
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)

! Determinant-Gauss:
det =1.

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

det =det*(-1)**n_sch_linii

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

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))
Hinv =0

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

if(is_Hilbert) then
call H_Inverse_Exact(n, Hinv)
else
call L_Inverse_Exact(n, Hinv)
endif

do i =1, n
!write(1,*) (Hinv(i,j), j =1,n)
write(1, '(<nn>(1pg16.7))') Hinv(i,1:n)
enddo


! 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)
enddo
enddo

!rel_err_max2 =maxval(abs(Rel_Error))
call amaxM(n, Rel_Error, rel_err_max, i, j)

write(*,*) rel_err_max      !, rel_err_max2
write(1,'(/, a, 1pE10.3)') 'Eroarea Relativa maxima (in modul) =', rel_err_max

write(1,'(a, I1,a,I1, a)') 'Element: (', i,',',j, ')'

! Numarul de conditie:
! Norma-1 (norma coloanelor):
! Cu Inversa EXACTA !:

cond1 =ANorm1*Norm1(Hinv,n)
write(1,'(/, a, 1pE10.2)') 'Numarul de conditie (cu inversa exacta) - Norma 1: Cond1 =', cond1

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

write(1,'(/, a,\)') 'Determinant (exact): Det = '

delta1 =delta(n)

if(.not. is_Hilbert) then   ! Lotkin
delta1 =delta1*(-1)**(n-1)
!write(1,*) 1./delta1

else    ! Hilbert
delta1 =delta1*n
!write(1,*) 1./delta1
endif

! Det. exact:
write(1,*) 1./delta1

write(1,'(a, \)') 'delta = 1/Det = '
write(1,*) delta1

! 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

