! Program principal pt. GAUSS-SEIDEL
! Cu matricile: "A(n,n)" & "b(n, nb)"
!
! Fisier de date (intrare/iesire):
! 1) Titlu
! 2) tip, kodInput
! 3) n, nb
! 4) [A] ... pentru kodA =1
! 5) [b] ... pentru kodA =1
! 6) kod_Inv, kod_PrintA, kod_Print
! 7) check_text
! 8) eps, nlim
!
! Fisier extern (intrare):				
! 1) Titlu =linie text
! 2) Matricea A
! 3) vectorul/vectorii b

! Semnificatia:
! Titlu: text <= 80 caractere
!
! tip: text <= 10 caractere; primul e semnificativ: s, S =simetrica; alt caracter =generala;
! kodInput:     Matricie A si b:
!           0 = generate (de Gen_Ab);
!           1 = citite din fisier de intrare;
!           2 = citita din fiser extern (altul decat fisierul de intrare)
!           alta valoare = invalida

! n:  nr. de ecuatii;
! nb: nr. de termeni liberi.
!
! A: Matricea coeficientilor n x n
!    Se scrie asa cum e data - pe linii si coloane
!    Pentru kodA#0: se introduce triunghiul superior (pe linii).
!
! b: Termenii liberi
!    Se scriu pe linii: "nb" linii a "n" termeni liberi.
!
! kod_PrintA: #0/0: Se tiparesc/Nu se tiparesc A si b
! kod_Print:  #0/0: Se tiparesc/Nu se tiparesc iteratiile.
! 
! 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".
!
! eps, nlim: toleranta, numarul limita de iteratii
!
Program Gauss_Seidel
use dfwin

real, allocatable :: a(:,:), b(:,:), b1(:)
real, allocatable :: x(:,:), x0(:), x1(:)

character(max_path) fname, fname2, dummy*80
real norm

logical Print_It

character(1) check_text
logical do_check, full_check 

character(30) text(-2:1)
integer width
integer status

character(10) tip
logical simetric, general
real prag/1E-6/
character(3) char3


text(-2) =' Iteratie stationara'
text(-1) =' Divergenta'
text(0) =' Toleranta atinsa.'
text(1) =' Numar de iteratii depasit.'

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

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

read(1,'(a80)') dummy	! Titlu

! Matricea A:
read(1,*) tip, kodInput
! Verifica kodInput:
if(Input_Check(kodInput) <0) goto 500

read(1,*) n, nb
allocate (A(n,n), b(n,nb))
A =0.; b =0.;

simetric =tip(1:1) =='s' .or. tip(1:1) =='S'
general =.not. simetric

If_2: IF(kodInput ==2) THEN
	call GetFile(fname2, status, 'Select Ab-Matrices File'C)
	open(2, file =fname2, status ='old', form ='formatted')
ENDIF If_2

!if(kod_Ab /=1) then     ! Generata de gen_Ab; sau citita din fisier exterior;
allocate(x(n,nb))
x =0.
!endif

allocate (b1(n))
allocate (x0(n), x1(n))
b1 =0;
x0 =0; x1 =0;

! Citeste/Genereaza:
call Read_Ab(n, nb, A, b, general, kodInput)

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

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

Print_It =kod_Print /=0

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

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

read(1,*) eps, nlim

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,jb), jb= 1,nb)
	enddo

endif


! Initializeaza "guess" vector:
!x0 =0.

n8 =min(n,8)

! Pivotare:
! Check for a(i,i) #0:
n_sch_linii =0
pvt_DO: do i=1,n
if( abs(a(i,i)) < prag) then
call pivotare(n, nb, a, b, prag, i, n_sch_linii, kod_pivot)
    if(kod_pivot <0) then
    write(1,*) 'Pivot < prag, linia:', i
    write(*,*) 'Pivot < prag, linia:', i
    goto 500
    endif
endif
enddo pvt_DO

iwn =width(n)
write(1, '(/,a, I<iwn>)') 'Numar de schimbari de linii: ', n_sch_linii

DO_NB: DO ib =1, nb
b1 =b(:,ib)
! Initializeaza "guess" vector:
x0 =0.

if(Print_It) then
itmax =n8*15 +6
write(1,'(/, a, t<itmax>,a)') ' Iter	x', 'max|x - x_ant|'
write(1, '(a, t6, (<n8>g15.7))') '  0', x0
endif

dif_ant =1E10
n_div =0
! Iteratie:
ITER_DO: DO ITER =1, NLIM

DO1: DO i=1,n

sum =0.
	do2: do j=1,n
	
	if(j<=i-1) then
	xx =x1(j)

	elseif( j==i) then
	cycle do2

	else	! j>i
	xx =x0(j)
	endif

	sum =sum +a(i,j)*xx
	enddo do2

x1(i) =(b1(i) - sum)/a(i,i)
enddo DO1


! Test eroare:
dif =norm(x1 -x0, n)

! Scrie iteratiile:
iwi =width(i)+2
if(Print_It) write(1, '(I<iwi>, t6, (<n8+1>g15.7))') iter, x1, dif

if( dif <=eps) then
kod =0
exit ITER_DO

elseif(dif < dif_ant) then
kod =1
x0 =x1
dif_ant =dif

elseif(dif >= dif_ant) then
!DIVERGENTA:
n_div =n_div+1

if(dif ==dif_ant) then
write(1,*) '* Iteratie stationara', iter
kod =-2
exit ITER_DO

else
    if(Print_It) then
    write(1,*) '* Divergenta'
    else
    write(1,*) '* Divergenta la iteratia', iter
    endif
endif

    if(n_div ==2) then
    KOD =-1
    exit ITER_DO
    else
    x0 =x1
    dif_ant =dif
    endif
endif

ENDDO ITER_DO

if(nb>1) then
iwib =width(ib)
write(1,'(/, a, i<iwib>)') ' Solutia - Termenii liberi nr. ', ib
else
write(1,'(/, a)') ' Solutia:'
endif

! Incheierea iteratiei cf. "kod":
if(kod >=0) then
write(1, '(/, a, a, 1pg12.6, a)') trim(text(kod)), ' (', dif, ')'
write(*, '(/, a)') trim(text(kod))

else        !
write(*, '(/, a)') trim(text(kod))
    if(kod ==-1) then
    char3 =' > '
    else
    char3 =' = '
    endif
write(1, '(/, a, a, 2(1pg12.6, a))') trim(text(kod)), ' (', dif, char3, dif_ant, ')'
endif

! Scrie solutia:

if(kod ==0) then
iwiter =width(iter)
write(1, '(a, I<iwiter>,/)') ' Numar iteratii: ', iter
write(*, '(a, I<iwiter>, /)') ' Numar iteratii: ', iter

!write(1,'(/, a)') ' Solutia:'

else
write(*,*)
write(1, '(/, a)') ' Ultima aproximatie:'
endif

do i=1,n
iwi =width(i)
write(1, '(a, i<iwi>, a, g15.7)') ' x(', i, ') =', x1(i)
enddo

x(:, ib) =x1
ENDDO DO_NB

! Proba: Solutia
CHECK_IF: IF(do_check) THEN

call Proba(n, nb, a, b, x, full_check)

ENDIF CHECK_IF

500 continue
END Program Gauss_Seidel

