Subroutine LUdecomp_z(A, n, order, do_pvt, kodPivot, kod)
! Descompunerea LU a matricii A. 
!* L are (1,1) pe diagonala (DOOLITLLE). 
! Stocare compacta: L, U in A. 
! Pivotare partiala - prin apel APVT
! Ordinea liniilor (dupa pivotare) e stocata in vectorul ORDER:
! Ex.: ORDER(1) = 3: linia initiala "3" este in pozitia 1.
! ORDER se utilizeaza pentru ordonarea membrilor II - la rezolvare (in LUsolve)
!
! Parametri:
! A = matricea sistemului
! n = numarul de ecuatii
! Order = vector intreg
! kod = cod de test pivot: kod =0: o.k.; kod <0: pivot <prag
! Indici: i* =linie; j*, k* =coloana
! Produsele in dubla precizie

complex A(n, n)
integer order(n)
real prag /1E-7/	! prag pentru Pivot
complex(8) sum
integer width

common no_interchange, kodPrint
logical do_pvt


kod =0

! Initializeaza ORDER:
do i =1, n
order(i) =i
end do

!no_interchange =0 - in Main

! Pivoteaza pentru col. 1:
Pvt1_If: if(do_pvt) then
call APVT_z(A, n, order, 1, kodPivot, prag)
endif Pvt1_If

	if(kodPrint >=1) then
	iwn =width(n)
	write(1, '(/, a, <n>(I<iwn>,1x))') 'Pas 1 - ORDER: ', order
		if(kodPrint >=2) write(1,*) a(1,1)
	endif

! Mesaj pentru Pivot < prag:
if (abs(a(1, 1)) < prag) then
write (*,*) ' u(1) =', a(1,1), ' <', prag
write (1,*) ' u(1) =', a(1,1), ' <', prag

kod =-1
return
endif

! Linia 1 - U: Cf. L(1,1)=1, linia 1, j>=1, ramane nemodificata U(1,j) =a(1,j):
! Coloana 1 - L:
do irow =2,n
a(irow, 1) =a(irow, 1)/a(1,1)
enddo


! Calculeaza L si U: 
MAIN_DO: DO irow =2, n
! U - pe linii:
! Cautare pivot si test prag -inainte de calculul liniei -U:
! Nu - pentru ultima linie:
IF(irow <n) THEN

Pvt_If: if(do_pvt) then
call APVT_z(a, n, order, irow, kodPivot, prag)
endif Pvt_If

ENDIF


! Linia "irow" a lui U:
i1 = irow-1
do jcol =irow, n
sum =(0D0, 0D0)
	do k =1, i1
	sum =sum +dcmplx(a(irow, k))*dcmplx(a(k, jcol))
	enddo
a(irow, jcol) =a(irow, jcol) -sum		! Cf.: L(i,i) =1.
enddo

	if(kodPrint >=1) then
	iwn =width(n); iwr =width(irow);
	write(1, '(/, a, I<iwr>, a, <n>(I<iwn>, 1x))') 'Pas ' , irow, ' - ORDER: ', order(irow:n)
		if(kodPrint >=2) write(1,*) a(irow,irow)
	endif

! Check for u(i,i) =0:
if(abs(a(irow, irow)) <prag) then

iwi =width(irow)
write(1, '(/,a, I<iwi>, a, I<iwi>, a, 1pg15.7, 1pg15.7)') ' u(', irow, ',', irow, ') = ', a(irow, irow)
kod =-1
RETURN
endif

! Calculeaza coloana "irow" a lui L: irow =2,...,n:
i2 = irow +1
	do iL =i2, n
sum =(0D0, 0D0)
		do k =1, i1
		sum =sum +dcmplx(a(iL, k))*dcmplx(a(k, irow))
		enddo
	a(iL, irow) =(a(iL, irow) -sum)/a(irow, irow)
	enddo

ENDDO Main_DO

END

