Program Gauss_AB
! Program principal pt. ELIM
! Versiunea cu matricea extinsa AB
! Cu proba.
!
! Intrare = Fisier de date:
! 1) Titlu: 80 caractere
! 2) n, nb: nr. de ecuatii; nr. de termeni liberi.
! 3) Matricea coeficientilor "A": n x n (Se scrie asa cum e data - pe linii si coloane).
! 4) Termenii liberi "b": Se scriu in linii: "nb" linii a "n" termeni liberi.
! 5) Cod inversare coeficienti: 0/#0. Pentru cod #0: a(i,j) -> 1./a(i,j)
! 6) Cod Proba: text<= 20 caractere. Pentru proba, primul caracter este: "C","c" sau "P","p".

real, allocatable :: AB(:,:)
character(260) fname, titlu*80
real prag /1E-7/        !-8/	! Prag pentru Pivot ~0.
integer width

real, allocatable :: ab0(:,:)
character(20) check_text; logical do_check; 
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,nb

ncol =n+nb
allocate (AB(n, ncol))

! citeste matricea coeficientilor:
read(1,*) ((AB(i,j),j=1,n),i=1,n)

! citeste vectorii termenilor liberi - in numar de NB:
read(1,*) ((AB(i,n+j),i=1,n),j=1,nb)

! citeste 'kod':
read(1,*) kodInv

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

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'


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

! Pt. Proba:
if(do_check) then
allocate(ab0(n,ncol))
ab0 =ab
endif

!	NCol =n+nb
	call ELIM(AB, N, NCol, prag, kod, pivot)

if(kod <0) then
write(1, '(a, I3, a,1pE8.2, a,1pE8.2,a, /,a)') '* Linia', -kod,': Pivot =', pivot, ' <', prag,'!',&
'* Sistemul nu se poate rezolva.'
goto 500
endif

iwnb =width(nb)

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, ') = ', AB(i,n+k) 
	enddo n_do
enddo nb_do

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

DO kb =1,nb
iwk =width(kb)
if(nb >1) then
write(1, '(/, a, I<iwk>, a)') 'Solutia nr. ', kb, ': A*x -b'
else
if(nb >1) write(1, '(/, a)') ': A*x -b'
endif

!x =b(1:n,kb)

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

write(1,*) i, sum -ab0(i, n+kb)
enddo

ENDDO
ENDIF Check_if

500 deallocate(AB)
if(allocated(ab0)) deallocate(ab0)


END Program Gauss_AB

