! Metoda GAUSS-SEIDEL cu accelerare SOR
! Cu scalarea matricii A(n,n) si vectorului b(n).
!
! Fisier de date (intrare):
! 1) Titlu
! 2) n, kod_Ab
! 3) [A]
! 4) [b]
! 5) kod_PrintA, kod_Print
! 6) check_text
!
! Titlu: text <= 80 caractere
! n: nr. de ecuatii
! kod_Ab: Cod de intare a matricilor "A" si "b":
!		kod_Ab <=0: matrici A,b, generate - de subrutina GenAB;
!		kod_Ab =1: matrici citite din fisierul de intrare
!
!		kod_Ab =2: matrice A generala, si b, citite dintr-un fisier de date;
!		kod_Ab >=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 "A": n x n (Se scrie pe linii si coloane).
! b: Termenii liberi "b": "n" termeni liberi (pe linie).
! eps, nlim: toleranta, numarul limita de iteratii
!
! 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".
!

Program SOR

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

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

logical Print_It

character(20) check_text
logical do_check, full_check 

character(30) text(0:1)

integer width
real max_dif


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	! Comments

read(1,*) n, kod_Ab
allocate (a(n,n), b(n), x0(n), x1(n))

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

elseif (kod_Ab ==1) then
	do i =1,n
	read(1,*) a(i,:)
	enddo

	read(1,*) b


else	! kod_Ab >=2
	call GetFile(fname2, status, 'Select Ab_Matrices File'C)
	open(2, file =fname2, status ='old', form ='formatted')

	read(2,*) dummy
		sim_if: if(kod_Ab ==2) then
! A:
		do i =1,n
		read(2,*) a(i,:)
		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:
	read(2,*) b

	close(2)

endif


read(1,*) eps, nlim
read(1,*) kod_PrintA, kod_Print

read(1,*) omega
omega_1 =1 -omega

read(1,*) check_text

Print_It =kod_Print /=0

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

endif


if(Print_It) write(1,'(/, a)') 'Iter	x' 

! Scale A and b:
do i =1,n
b(i) =b(i)/a(i,i)
	j_do: do j=1,n
	if(j==i) cycle j_do
	a(i,j) =a(i,j)/a(i,i)
	enddo j_do
a(i,i) =1.
enddo

! Init "guess" vector:
x0 =0.
kod =1

n8 =min(n,8)

! Begin iteration:
ITER_DO: DO ITER =1, NLIM

i_DO: DO i=1,n
sum1 =0.; sum2 =0.
	j1_do: do j=1,i-1
	sum1 =sum1 +a(i,j)*x1(j)
	enddo j1_do

	j2_do: do j=i+1,n
	sum2 =sum2 +a(i,j)*x0(j)
	enddo j2_do

xGS =b(i) - sum1 -sum2					!(b(i) - sum)/a(i,i)
if(i ==1) then
x1(i) =xGS
else
x1(i) =omega*xGS +omega_1*x0(i)
endif
enddo i_DO

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

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

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

else
kod =1
x0 =x1
endif

ENDDO ITER_DO

! Incheierea iteratiei cf. "kod":
write(1, '(/, a, a, 1pg12.6, a)') trim(text(kod)), ' (', dif, ')'

! Scrie solutia:

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

else
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

! 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



max_dif =0.

do i=1,n
sum =0.
	do k=1,n
	sum =sum +a(i,k)*x1(k)
	enddo

dif =sum -b(i)
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


ENDIF Check_if


END Program SOR

