Program LU_complex
! A = Complexa.
! Programul principal pentru LU (DOOLITLLE)
! Apel: LUdecomp; LUsolve; APVT.
! Cu Proba: (1) Az-b ~0; (2) LU =PA.
!
! Intrare = Fisier de date:
! 1) Titlu: 80 caractere
! 2) n, nb: nr. de ecuatii; nr. de termeni liberi.
! 3) Matricea coeficientilor "A": n x 2n (Se scrie pe linii. Input real: real(a_ij), imag(a_ij)).
! 4) Termenii liberi "b": Se scriu pe linii: "nb" linii a "n" termeni liberi = numere complexe.
! 5) Cod_Tiparire, Cod_inversare coeficienti: 0/#0. 
!	 - Cod_Tiparire #0: Se tiparesc: Solutia - la display; Matricile L si U - in fisierul de date.
!                       kodPrint =2: se tipareste si vectorul "order", la fiecare pivotare. 
!    - Cod_Inversare #0: a(i,j) -> 1./a(i,j)
! 6) Cod pivotare (logic), kod_pivot: 
!		Cod_pivotare: T =Da; #T =Nu 
!		kodPivot: #0 = pivot maxim; 0 = primul pivot #0
! 7) Text_proba (primul caracter e semnificativ): C,c, P,p = Proba (1); F,f =Proba (1) si (2).
 
complex, allocatable :: A(:,:), b(:,:),x1(:)
integer, allocatable :: order(:)
character(260) fname, titlu*80
integer width

complex, allocatable :: l(:,:), u(:,:), c(:,:)
complex dl

common no_interchange, kodPrint
complex diag

real, allocatable :: ar(:,:)
complex, allocatable :: a0(:,:), b0(:,:)

complex sum

character(1) check_text; logical do_check, full_check;

logical do_pvt
real, allocatable :: be(:)

! Determinant:
complex(8) det


print '(a,\)', ' Fisier de intrare: '
read '(a)', fname
open(1,file=fname,status='old', form ='formatted')

read(1,'(a)') TITLU

! Citeste n, nb:
read(1,*) n, nb

allocate (A(n,n), b(n, nb), x1(n))
allocate(order(n))

! citeste matricea coeficientilor:
!read(1,*) ((A(i,j),j=1,n),i=1,n)

allocate(ar(n, 2*n))

read(1,*) (ar(i,:), i=1,n)

do i=1,n
	do j=1,n
	j2 =j*2
	a(i,j) =cmplx(ar(i,j2-1), ar(i,j2))
	enddo
enddo

deallocate(ar)

! citeste vectorii termenilor liberi - in numar de NB:

read(1,*) ((b(i,kb),i=1,n), kb=1,nb)

! citeste 'kod':
read(1,*) kodPrint, kodInv

read(1, *) do_pvt, kodPivot

read(1, '(a1)') check_text

do_check =check_text =='C' .or.check_text =='c' .or. &
		check_text =='P' .or. check_text =='p' .or. &
		check_text =='F' .or. check_text =='f'

full_check = check_text =='F' .or. check_text =='f'

IF(kodInv /=0) THEN
do i=1,n
	do j=1,n
	a(i,j) =(1.,0.)/a(i,j)
	enddo
enddo
ENDIF

! Pt. Proba 1-2:
if(do_check .or. full_check) then
allocate(a0(n,n), b0(n,nb))
a0 =a; b0 =b;
endif



! Descompunera LU:
! dl =(1,0)	- in LUdecomp
no_interchange =0
kodLU =0

call LUdecomp_z(A, n, order, do_pvt, kodPivot, kodLU)
if(kodLU <0) then
	write(*,*) 'U-diagonal < Threshold. LU decomposition not feasible.'
	write(1,*) 'U-diagonal < Threshold. LU decomposition not feasible.'
call EXIT
endif

! Solutia pentru "nb" termeni liberi:
do kb =1, nb
x1 =0.
call LUsolve_z(A, b(1:n,kb), n, x1, order)
b(1:n, kb) =x1

	if(kb ==1 .and. kodPrint /=0) then
	iinter =width(no_interchange)+1
	write(1, '(/, a, I<iinter>)') 'Numar de schimbari de linii =', no_interchange
	endif
enddo

Print_it: IF(kodPrint /=0) THEN
iwn =width(n)
diag =(1.,0.)


! Write L:
write(1, '(/, a)') 'Matricea L'
do i=1,n
write(1, '(<i>(g15.7,g15.7))') (a(i,j), j =1, i-1), diag
enddo

! Write U:
write(1, '(/, a)') 'Matricea U'
do i=1,n
!it =(i-1)*30
write(1,'(<n-i+1>(g15.7, g15.7))') (a(i,j), j =i,n)
enddo

ENDIF Print_it

! Calculeaza determinantul - pentru kod =2:
IF(kodPrint ==2) THEN
det =(1D0,0D0)
do i=1,n
det =det*a(i,i)
enddo


if(mod(no_interchange,2) /=0) then
det =-det
endif

write(1, '(/, a, (1pg15.7, 1pg15.7))') 'Determinant(A) =', det
ENDIF

! Scrie solutia:
iwnb =width(nb)

nb_do: do kb =1,nb
write(1, '(/, a, I<iwnb>)') ' * Solutia pentru termenii liberi nr. ', kb
		if(kodPrint /=0) print '(/, a, I<iwnb>)', ' * Solutia pentru termenii liberi nr. ', kb
	n_do: do i=1,n
	iwi =width(i)
	write(1, '(a, I<iwi>, a, (g15.7, g15.7))') ' x(', i, ') = ', b(i,kb)
		if(kodPrint /=0) print '(a, I<iwi>, a, (g15.7, g15.7))', ' x(', i, ') = ', b(i,kb)
	enddo n_do
enddo nb_do

! Proba-1: Solutia
Check_If: IF(do_check) THEN

write(1, '(/,a)') '* Proba Ax =b:'

DO kb=1,nb
iwk =width(kb)
if(nb >1) then
write(1, '(/, a, I<iwk>, a)') 'Solutia nr. ', kb, ': A*x -b'
else
write(1, '(/, a)') ' A*x -b'
endif

x1 =b(1:n,kb)

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

write(1,*) i, sum -b0(i,kb)
enddo

ENDDO

ENDIF Check_If

! Proba-2: LU =PA
LU_check: IF(full_check) THEN
allocate(l(n,n), u(n,n), c(n,n))

l =(0.,0.); u =(0.,0.); c=(0.,0.);
do i=1,n
	do j=1, i-1
	l(i,j) =a(i,j)
	enddo
l(i,i) =diag
enddo

do i=1,n
	do j=i,n
	u(i,j) =a(i,j)
	enddo
enddo


c =matmul(L,U)
write(1, '(/,a)') '* Proba LU =PA: L*U'
do i =1,n
if(kodInv ==0) then
write(1,*) (c(i,j), j=1,n);
else
write(1,*) (1./c(i,j), j=1,n)
endif
enddo

ENDIF LU_check

!deallocate(A, b, x, order, c, l, u)
deallocate(A)
deallocate(b)
deallocate(x1)
deallocate(order)

if(allocated(c)) deallocate(c)
if(allocated(l)) deallocate(l)
if(allocated(u)) deallocate(u)

END Program LU_complex

