Program Nr_Conditie
!  Program principal pt. Numarul de conditie
! Matricea A = matricea data (n,n); Matricea B = matricea inversa (n,n).
! Norma-infinit; Norma-1. ; Norma-2; Cond(A)*.
!
! Fisier de date (intrare/iesire)
! 1) Titlu: 80 caractere
! 2) n: nr. de ecuatii
! 3) Matricea coeficientilor "A": n x n (Se scrie asa cum e data - pe linii si coloane)
!
! 4) Cod_inversare coeficienti, Cod_Tiparire: 0/#0. 
!    - Cod_Inversare coeficienti: #0 => a(i,j) -> 1./a(i,j)
!    - Cod_Tiparire: #0: se tipareste si A^(-1)
! 5) Cod_Norma
!	 - Cod_Norma: 0 = norma-infinit (norma liniilor)
!				  1 = Norma-1 (norma coloanelor)
!				  2 = Norma-2 (norma euclidiana).
!				 <0 = Cond(A)*
!       Alta valoare: Se ia norma-infinit
!
real, allocatable :: A(:,:), B(:,:), A0(:,:), ATxA(:,:)
character(80) titlu, fname*260
real Norm1, Norm$, Norm2
real prag /1e-8/        !7/    !/1E-6/

complex eigen_max, eigen_max_1
integer status

!print '(a,\)', ' Nume fisier de intrare: '
!read '(a)', fname
! Selecteza fisierul de intrare:
call GetFile(fname, status, 'Select Input File'C)

open(1, file =fname, status ='old', form ='formatted')

read(1,'(a)') TITLU
! Citeste n, nb:
read(1,*) n

allocate (A(n,n), B(n,n), A0(n,n))
A =0.0; B =0.0; A0 =0.

! Citeste matricea coeficientilor (pe linii):
read(1,*) ((A(i,j),j=1,n),i=1,n)
A0 =A

read(1,*) kod_inv, kod_Print	

read(1,*) kod_norma		! 0 =infinit; 1 =norma-1; 2 =norma-2; <0: Cond()*

if(kod_Inv /=0) then
do i=1,n
do j=1,n
a(i,j) =1./a(i,j)
enddo
enddo
endif

! Atribuie vectorii termenilor liberi = vectorii unitate
do i =1,n
B(i, i) =1.0
enddo

! Calculeaza norma lui A (inainte de eliminare):
select case(kod_norma)
case(:-1)
call Eigen_Power(n, A, eigen_max, iter_flag, test_val)

if(iter_flag <0) then
write(1, '(/,a)') ' A: Toleranta 1E-5 la "lambda", nu s-a atins in 1000 de iteratii (metoda Puterii).'
write(1, '(a, 1pg11.4)') ' Valoarea de test =', test_val
write(1, '(a)') ' Calculeaza valorile proprii ale lui A cu alta metoda !'
endif

C =abs(eigen_max)

case(0)
C = Norm$(n, A)

case(1)
C = Norm1(n, A)

case(2)
allocate(ATxA(n,n))
ATxA =matmul(transpose(A),A)

call Eigen_Power(n, ATxA, eigen_max, iter_flag_2, test_val)

if(iter_flag_1 <0) then
write(1, '(/,a)') ' A_transpus*A: Toleranta 1E-5 la "lambda", nu s-a atins in 1000 de iteratii (metoda Puterii).'
write(1, '(a, 1pg11.4)') ' Valoarea de test =', test_val
write(1, '(a)') ' Calculeaza valorile proprii ale lui "A_transpus*A" cu alta metoda !'
endif

C =sqrt( abs(eigen_max) )

case default
C = Norm$(n, A)
end select

! nb =n
	call ELIM2(n, n, A, B, prag, kod)

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

! Tipareste matricea inversa:
if(kod_Print /=0) then
write(1, '(/, a)') ' Matricea inversa:'

do i =1,n
write(1, '(<n>G15.7)') B(i,:) 
enddo
endif

! Norma inversei:
write(1,*)
select case(kod_norma)
Case(:-1)
call Eigen_Power(n, B, eigen_max_1, iter_flag_1, test_val)

if(iter_flag_2 <0) then
write(1, '(/,a)') ' A-invers: Toleranta 1E-5 la "lambda", nu s-a atins in 1000 de iteratii (metoda Puterii).'
write(1, '(a, 1pg11.4)') ' Valoarea de test =', test_val
write(1, '(a)') ' Calculeaza valorile proprii ale lui "A-invers" cu alta metoda !'
endif

Cinv =abs(eigen_max_1)

write(1, '(/, a)') ' Cond(A)*:'

case(0)
write(1, '(/, a)') ' Norma-infinit (norma liniilor):'
Cinv= Norm$(n, B)

case(1)
write(1, '(/, a)') ' Norma-1 (norma coloanelor):'
Cinv= Norm1(n, B)

case(2)
write(1, '(/, a)') ' Norma-2 (norma euclidiana):'

ATxA =matmul(transpose(B),B)

call Eigen_Power(n, ATxA, eigen_max_1, iter_flag_3, test_val)

if(iter_flag_3 <0) then
write(1, '(/,a)') ' A-invers_transpus*A-invers: Toleranta 1E-5 la "lambda", nu s-a atins in 1000 de iteratii (metoda Puterii).'
write(1, '(a, 1pg11.4)') ' Valoarea de test =', test_val
write(1, '(a)') ' Calculeaza valorile proprii ale lui "A-invers_transpus*A-invers" cu alta metoda !'
endif

Cinv =sqrt( abs(eigen_max_1) )

case default
write(1, '(/, a)') 'Norma nedefinita: Se ia Norma-infinit (norma liniilor)'
Cinv= Norm$(n, B)
end select


! Calculeaza numarul de conditie:
Cond =C*Cinv

if(kod_norma >=0 .and. kod_norma /=2) then
write(1, '(/, a, 1pg15.7)') ' Norm(A) =', C
write(1, '(a, 1pg15.7)') ' Norm(A-invers) =', Cinv

elseif(kod_norma ==2) then

	if(iter_flag_2 <0) then
	write(1, '(/, a, 1pg15.7, a)') ' Raza_spectrala(A_transpus*A) =', abs(eigen_max),' (?)'
	else
	write(1, '(/, a, 1pg15.7)') ' Raza_spectrala(A_transpus*A) =', abs(eigen_max)
	endif

	if(iter_flag_3 <0) then
	write(1, '(a, 1pg15.7, a)')    ' Raza_spectrala(A-invers_transpus*A-invers) =', abs(eigen_max_1), " (?)"
	else
	write(1, '(a, 1pg15.7)')    ' Raza_spectrala(A-invers_transpus*A-invers) =', abs(eigen_max_1)
	endif

else
	if(iter_flag <0) then
	write(1, '(/, a, 1pg15.7, a)') ' Raza_spectrala(A):', C, ' (?)'
	else
	write(1, '(/, a, 1pg15.7)') ' Raza_spectrala(A):', C
	endif

	if(iter_flag_1 <0) then
	write(1, '(a, 1pg15.7, a)')    ' Raza_spectrala(A-invers):', Cinv, ' (?)'
	else
	write(1, '(a, 1pg15.7)')    ' Raza_spectrala(A-invers):', Cinv
	endif
endif

if(Cond <=1000.) then
assign 10 to iFMT
else
assign 11 to iFMT
endif

print '(/,a)', ' '//trim(titlu)
!print '(a, 1pG15.7)', " Numarul de conditie =", Cond
write(*, iFMT) " Numarul de conditie =", Cond


write(1, '(a)')

IF(kod_norma >=0 .and. kod_norma /=2) THEN
!write(1, 10) " Numarul de conditie =", Cond
write(1, iFMT) " Numarul de conditie =", Cond

ELSEIF(kod_norma ==2) THEN
	if(iter_flag_2 >=0 .and. iter_flag_3 >=0) then
!	write(1, 10) " Numarul de conditie =", Cond
write(1, iFMT) " Numarul de conditie =", Cond

	else
	write(1, 10) " Numarul de conditie = ?"
	endif


ELSE
	if(iter_flag >=0 .and. iter_flag_1 >=0) then
!	write(1, 10) " Numarul de conditie =", Cond
write(1, iFMT) " Numarul de conditie =", Cond

	else
	write(1, 10) " Numarul de conditie = ?"
	endif

ENDIF

10 format(a, :, 1pg15.7)
11 format(a, :, 1pE15.7)


500 deallocate(A, B)

END program Nr_Conditie

