Program Gauss
! Program principal pt. ELIM
! Cu matricile: "A(n,n)" si "b(n,nb)"
! A. b: generate; citite din fisierul de intrare; citite din din fisier de date.
!
! Ian.2013: 
! Afiseaza si (v. kod_Print):
! - matricile L,U
! - determinantul lui A
! - termenii liberi g
!
! Selectia fisierului de intrare - prin fereastra de dialog standard.
!
! Fisier de date (intrare):
! 1) Titlu
! 2) n, nb, kodAb
! 3) [A]
! 4) [b]
! 5) kod_Inv, kod_PrintA, kod_Print
! 6) check_text
!
! Semnificatia:
! Titlu: text <= 80 caractere
!
! n, nb: nr. de ecuatii; nr. de termeni liberi.
! kodAb: Cod de intare a matricilor "A" si "b":
!		kodAb <=0: matrici A,b, generate - de subrutina GenAB;
!		kodAb =1: matrici citite din fisierul de intrare
!
!		kodAb =2: matrice A generala, si b, citite dintr-un fisier de date;
!		kodAb >=3:matrice A simetrica, si b, citite dintr-un fisier de date
!				Prima linie in fisierul specificat e o linie text (max. 80 caractere).
!
! A: Matricea coeficientilor n x n
!    Se scrie asa cum e data - pe linii si coloane
!    Pentru kod_Ab >=3: se introduce triunghiul superior (pe linii).
!
! b: Termenii liberi
!    Se scriu pe linii: "nb" linii a "n" termeni liberi.
!
! kod_Inv: Cod_inversare_coeficienti
!	      kod_Inv: #0: a(i,j) -> 1./a(i,j)
! kodPrintA: cod tiparire matrici
!	kod_PrintA: #0: Se tiparesc matricile A si b
! kod_PrintP: cod_Tiparire
!	kod_PrintP: #0: Se tiparesc permutarile si pivotii.
!   kod_PrintP: >=2: Se tiparesc matricile "L" si "U"; si Determinantul.
!   kod_PrintP: >=3: Se tiparesc si termenii liberi "g". 
!   kod_PrintP: >=4: Se face proba: LU =A, sau A-permutat.
!
! Check_text: text <= 20 caractere. Primul caracter e semnificativ:
!	'C', 'c', 'P', 'p', 'F', 'f': se face proba "Ax-b" (pentru fiecare termen liber).
!   C,c, P, p: se tipareste diferenta maxima in modul "Ax-b"
!
!   F,f: - se tiparesc diferentele "Ax-b";
!
real, allocatable :: A(:,:), B(:,:)
character(260) fname, fname2, titlu*80

real prag /1E-8/    !7/	! Curent /1E-6/; Pentru test:/1E-7/; Prag pentru Pivot ~0.
integer width

real, allocatable :: a0(:,:), b0(:,:)

character(20) check_text
logical do_check, full_check 

character(80) dummy
character(40) text

real max_dif
integer status
real, allocatable :: xLU(:), row(:), col(:)

! Selecteza fisierul de intrare:
call GetFile(fname, status, 'Select Input File'C)

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

! citeste:
read(1,'(a)') TITLU

read(1,*) n, nb, kodAb

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

! citeste/genereaza A si b:
if (kodAb <=0) then
call GenAB(n, nb, a, b)

elseif (kodAb ==1) then
	do i =1,n
	read(1,*) (a(i,j), j= 1,n)
	enddo

	do j =1,nb
	read(1,*) (b(i,j), i= 1,n)
	enddo

else	! kodAb >=2
	call GetFile(fname2, status, 'Select Ab-Matrices File'C)
	open(2, file =fname2, status ='old', form ='formatted')

	read(2,*) dummy
		sim_if: if(kodAb ==2) then
! A:
		do i =1,n
		read(2,*) (a(i,j), j= 1,n)
		enddo

		else	! matrice simetrica: triunghiul superior
		do i =1,n
		read(2,*) ( a(i,j), j=i,n )
		enddo

		! completeaza simetric:
		do i =1,n
		do j =1,i
		a(i,j) =a(j,i)
		enddo
		enddo

		endif sim_if

! b:
	do j =1,nb
	read(2,*) (b(i,j), i= 1,n)
	enddo

	close(2)

endif

! citeste coduri:
read(1,*) kod_Inv, kod_PrintA, kod_PrintP

! citeste "Proba"
read(1,*) check_text

full_check =check_text(1:1) =='F' .or. check_text(1:1) =='f'

do_check =check_text(1:1) =='C' .or. check_text(1:1) =='c' &
	 .or. check_text(1:1) =='P' .or. check_text(1:1) =='p' &
	 .or. full_check

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

! 
if(kod_PrintA /=0) then
! write A:
write(1,'(/,a)') 'Matricea A'
	do i =1,n
	write(1,*) (a(i,j), j= 1,n)
	enddo

! write b:
write(1,'(/,a)') 'Termenii liberi b'

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

endif

! Pt. Proba:
if(do_check) then
allocate(a0(n,n), b0(n,nb))
a0 =a; b0 =b;
endif

n_sch_linii =0
	call ELIM(n, nb, A, B, prag, kod, kod_PrintP, n_sch_linii)

if(kod <0) then
write(1, '(/, a, 1pE8.2,a, /,a)') '* Pivot <', prag,' !',&
'* Sistemul nu se poate rezolva.'

write(*, '(/, a, 1pE8.2,a, /,a)') '* Pivot <', prag,' !',&
'* Sistemul nu se poate rezolva.'

goto 500
endif

! Nr. schimb. de linii:
if(kodPrint /=0) then
i_nsl =width(n_sch_linii)
write(1,'(/,a, I<i_nsl>)') 'Nr. de pivotari: ', n_sch_linii
endif

! Matricile L si U ~ A:
Print_LU_Det: if(kod_PrintP >=2) then
write(1,'(/,a)') 'Matricea "L":'
	do i =1,n
	write(1,*) a(i,1:i-1), 1.E0
	enddo
	
write(1,'(/,a)') 'Matricea "U":'
	do i =1,n
	write(1,*) a(i,i:n)
	enddo

det =1.

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

!det =det*(-1)**n_sch_linii
if(mod(n_sch_linii,2) /=0) then
det =-det
endif
write(1,'(/,a, 1pE14.6)') 'Determinant(A) =', det

endif Print_LU_Det

! write b:
Print_g: if(kod_PrintP >=3) then
write(1,'(/,a)') 'Termenii liberi "g":'

	do i =1,n
	write(1,*) b(i,1:nb)
	enddo
endif Print_g

Print_LU: if(kod_PrintP >=4) then
allocate(xLU(n), row(n), col(n))
xLU =0; row =0.; col =0;


if(n_sch_linii ==0) then
write(1,'(/,a)') ' * Proba L*U =A:'
else
write(1,'(/,a)') ' * Proba L*U =A-permutat:'
endif

if(kod_Inv ==0) then
write(1, '(a)') '   Matricea L*U'
else
write(1, '(a)') '   Matricea 1./L*U'
endif

do i =1,n
if(i >=2) row(1:i-1) =a(i,1:i-1); 
row(i) =1.; row(i+1:n) =0;

xLU =0.
    do j=1,n
    col(1:j) =a(1:j,j); col(j+1:n) =0.;

    kij =min(i,j)
    xLU(j) =DOT_PRODUCT(row(:kij), col(:kij))
    enddo

!write(1,*) xLU(:)    	
if(kod_Inv ==0) then
write(1,*) xLU(:)    	
else
write(1,*) 1./xLU(:)    	
endif

enddo
endif Print_LU


! SOLUTIA:
iwnb =width(nb)
write(1, '(a)')

nb_do: do k =1,nb
write(1, '(/, a, I<iwnb>)') ' * Solutia pentru termenii liberi nr. ', k

	n_do: do i=1,n
	iwi =width(i)
	write(1, '(a, I<iwi>, a, G15.7)') ' x(', i, ') = ', B(i,k)
	enddo n_do
enddo nb_do


! Proba: Solutia
Check_if: IF(do_check) THEN
write(1, '(//,a)') ' * Proba'

if(full_check) then
write(1, '(/,a)') ' Diferente (A*x-b)'

else
write(1, '(/,a)') ' Diferenta (A*x-b), maxima in modul'
endif


b_do: DO kb =1,nb
iwk =width(kb)

if(nb >1) then
write(1, '(/, a, I<iwk>, a)') ' Solutia nr. ', kb, ':'
endif

!x =b(1:n,kb)

max_dif =0.
imax =1

do i =1,n
sum =0.
	do k=1,n
	sum =sum +a0(i,k)*b(k,kb)		!x(k)
	enddo

dif =sum -b0(i,kb)
if(full_check) write(1,*) i, dif

!max_dif =max(abs(dif), max_dif)
adif =abs(dif)
if(adif >max_dif)then
max_dif =adif
dif_max =dif
imax =i
endif

enddo

if(full_check) then
write(1, '(/, a, \)') ' Diferenta (A*x-b), maxima in modul:'
endif

	if(imax>0) then
	iwm =width(imax)+1
	write(1, '(1pg15.7, a, i<iwm>)') dif_max, '; ecuatia', imax

	else
	write(1,*) dif_max
	endif

ENDDO b_do

! Proba LU =A

ENDIF Check_if


500 deallocate(A, b)
if(allocated(a0)) deallocate(a0)
if(allocated(b0)) deallocate(b0)

END Program Gauss

