! Metoda GAUSS-SEIDEL cu accelerare SOR
! Cu scalarea matricii A si vectorului b.
!
! Intrare = Fisier de date:
! 1) Titlu: text <= 80 caractere
! 2) n: nr. de ecuatii
! 3) Matricea coeficientilor "A": n x n (Se scrie pe linii si coloane).
! 4) Termenii liberi "b": "n" termeni liberi (pe linie).
! 5) eps, nlim: toleranta, numarul limita de iteratii
! 6) Cod_Tiparire: #0/=0: Se tiparesc/Nu se tiparesc iteratiile.
! 7) check_text: text <= 20 caractere.Primul caracter e semnificativ:
!				 C, c, P, p: se face proba.

real, allocatable :: a(:,:), b(:), x0(:), x1(:)
character(128) fname, dummy*80
real norm
logical Print_It
character(20) check_text; logical do_check;
character(30) text(0:1)
integer width


text(0) =' Toleranta atinsa.'
text(1) =' Numar de iteratii depasit.'

print '(a, \)', 'Input file: '
read '(a)', fname
open(1, file=fname, status ='old')

read(1,'(a80)') dummy	! Comments
read(1,*) n
allocate (a(n,n), b(n), x0(n), x1(n))

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

read(1,*) b
read(1,*) eps, nlim
read(1,*) omega
omega_1 =1 -omega

read(1,*) kodPrint
read(1,*) check_text

! 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

Print_It =kodPrint /=0
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(Print_It) write(1,'(/, a)') 'Iter	x' 

! 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:
if(norm(x1 -x0,n) <=eps) then
kod =0
exit ITER_DO

else
x0 =x1
endif

ENDDO ITER_DO

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

! Scrie solutia:
!if(.not. Print_It) then

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

else
write(1, '(/, a, I<iwiter>)') ' Numar iteratii: ', iter
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

if(do_check) then
write(1, '(/, a)') ' Proba "Ax-b":'

do i=1,n
e =0.
	do j=1,n
	e =e +a(i,j)*x1(j)
	enddo
e =e -b(i)

iwi =width(i)+1
write(1, '(I<iwi>, 1pg16.7)') i, e
enddo

endif

END

