subroutine LUsolve(LU, b, n, x, order)
! Solutia lui Ax=b, dupa efectuarea desc. LU
! Pentru: DOOOLITLLE ~GAUS: L(i,i) =1.
! Elem. din "b" sunt re-aranjate cf. schimbarii de linii
! Solutia este returnata in "x"
! Parametri:
! LU = echivalentul LU a lui A (stocare compacta)
! b = vectorul term. liberi
! x = solutia
! order = vector intreg
implicit real(8) (a-h, o-z)

real(8) LU(n, n)
real(8) x(n), b(n)
integer order(n)
real(8) sum

! Re-aranjeaza elem, lui "b" - stocate in "x":
do i =1, n
j =order(i)
x(i) =b(j)
enddo


! Substitutie inainte: L*y =b	("b": stocat in "x")
! y(i) =(b(i) -SUM[j=1,i-1] l(i,j)*x(j))/l(i,i)
! x(1) =x(1), cf.: x(1) =x(1)/LU(1,1) si LU(i,i) =1.

row: do i =2, n
sum =0.D0
	do jcol =1, i-1
	sum =sum +dprod(LU(i, jcol), x(jcol))
	enddo
x(i) =(x(i) -sum)
enddo row

! Substitutie inapoi: U*x =y	("y": stocat in "x")
! x(i) =y(i) - SUM[j=i+1,n] u(i,j)*x(j)		! u(i,i) =1.
! Ultima ec.: LU(n,n)*x(n)=y(n).
x(n) =x(n)/LU(n,n)

do i = n-1, 1, -1
sum =0.D0
	do jcol =i+1, n
	sum =sum +LU(i, jcol)*x(jcol)      !dprod(LU(i, jcol), x(jcol))
	enddo
x(i) =(x(i) -sum)/LU(i,i)
enddo

END
