SUBROUTINE ELIM (n, nb, A, B, prag, kod, kodPrint)
! Rezolvitor pentru sisteme lniare, prin:
! Eliminare GAUSS, cu pivotare partiala.
! Termeni liberi multipli (numar ="nb") - in coloanele matricii "B".
! Parametri:
! n: Numarul ecuatiilor
! A(n,n): matricea coeficientilor
! B(n,nb): matricea termenilor liberi.
! Solutia este returnata in coloanele din "b".
! kod = 0/#0: O.K./ pivot ~=0 (< prag)

real A(n,n), B(n,nb), savea(n), saveb(nb)
real mik
integer width


if(kodPrint /=0) then
write(1,'(/, a)') 'Permutare, Pivot:'
endif

! REDUCERE: Linia curenta = "k"
Main_DO: DO k =1, n-1
! Indicele liniei pivotului
ipvt =k; pivot =a(k,k)
	do irow =k+1,n
	if(abs(a(irow,k)) >abs(pivot)) then
	ipvt =irow
	pivot =a(irow,k)
	endif
	enddo

! Test marime pivot: pivot =ab(ipvt, k)
	if(abs(pivot) < prag) then
	iwk =width(k)
	write(1, '(/, a, I<iwk>, a, 1pg15.7)') ' Linia ', k, ': Pivot =', pivot
	kod = -1		!k
	return
	endif

! Inter-schimba liniile "k" si "ipvt"
! Nu, daca ipvt =k:
if(ipvt /=k) then
savea =a(k,:); saveb =b(k,:)

a(k,:) =a(ipvt,:); b(k,:) =b(ipvt,:)
a(ipvt,:) =savea; b(ipvt,:) =saveb
endif


! Print:
if(kodPrint /=0) then
iwk =width(k)
	if(ipvt /=k) then
	write(1, '(I<iwk>, a, I<iwk>, 1pg15.7)') k, '<->', ipvt, a(k,k)
	else
	write(1,'(I<iwk>, 3x, 1pg15.7)') k, a(k,k)
	endif
endif

! Pivot =a(k,k)
! Reducere sub a(k,k), in coloanele kcol =k+1,...,n:
row: DO irow =k+1, n
! Nu, daca deja a(irow, k)=0:
if(a(irow, k) ==0) then
cycle row
else
mik =a(irow, k)/a(k,k)	!pivot
	acol: do jcol =k+1, n
	a(irow, jcol) =a(irow, jcol) -mik*a(k, jcol)
	enddo acol
	bcol: do jcol =1, nb
	b(irow, jcol) =b(irow, jcol) -mik*b(k, jcol)
	enddo bcol
endif
ENDDO row

Enddo Main_DO

! Test pentru a(n,n):
if(kodPrint /=0) then
iwn =width(n)
write(1, '(I<iwn>, 3x, 1pg15.7)') n, a(n,n)
endif

	if(abs(a(n,n)) < prag) then
	iwn =width(n)
	write(1, '(/, a, I<iwn>, a, 1pg15.7)') ' Linia ', n, ': Pivot =', a(n,n)
	kod = -1		!n
	return
	endif

kod =0
! BACK SUBSTITUTION:
! Ciclu pentru numarul de termeni liberi:
RHS: DO jRHS =1, nb
b(n, jRHS) =b(n, jRHS)/a(n,n)

do k = n-1, 1, -1
sum =0.
	do jcol =k+1, n
	sum =sum +a(k, jcol)*b(jcol, jRHS)
	enddo
b(k, jRHS) =(b(k,jRHS) -sum)/a(k,k)
enddo

ENDDO RHS

END subroutine Elim

