Program Crout_complex
! Programul principal pentru LU
! Apel: LUdecomp_Crout_z; LUsolve_Crout_z; APVT_z.
!
! Intrare = Fisier de date:
! 1) Titlu: 80 caractere
! 2) n, nb: nr. de ecuatii; nr. de termeni liberi.
! 3) Matricea coeficientilor "A": n x n (Se scrie asa cum e data - pe linii si coloane).
! 4) Termenii liberi "b": Se scriu pe linii: "nb" linii a "n" termeni liberi.
! 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.
!									 Se calculeaza si determinantul lui A. 
!    - Cod_Inversare #0: a(i,j) -> 1./a(i,j)
! 6) Cod pivotare (logic), kod_pivot: 
!		Cod_pivotare (T/F): T =Da; F =Nu 
!		kodPivot: #0 = pivot maxim; 0 = primul pivot #0
! 7) Check_text (20 car.; primul caracter este semnificativ): 
!    "C", "c", "P","p"/ alt caracter: Verifica/Nu verifica - solutia.
!    "F", "f": Se verifica solutia si LU =PA
!

complex, allocatable :: A(:,:), b(:,:),x(:)
integer, allocatable :: order(:)
character(260) fname, titlu*80
integer width

logical do_pvt

common no_interchange, kodPrint

real, allocatable :: ar(:,:)

complex, allocatable :: a0(:,:), b0(:,:)
character(20) check_text; logical do_check, full_check; 

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

complex(8) det


print '(a,\)', ' Nume 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), x(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

! citeste "Proba"
read(1,*) check_text

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' &
		.or. check_text(1:1) =='F' .or. check_text(1:1) =='f'
full_check =check_text(1:1) =='F' .or. check_text(1:1) =='f'

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

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


! Descompunerea LU:
no_interchange =0

call LUdecomp_Crout_z(A, n, order, do_pvt, kodPivot, kodLU)

if(kodLU /=0) then
	write(*,*) 'L-diagonal < Prag. Descompunerea LU nu e fezabila.'
	write(1, '(/, a)') ' L-diagonal < Prag. Descompunerea LU nu e fezabila.'
call EXIT
endif

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

	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.)

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

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

ENDIF Print_it

! Calculeaza determinantul - pentru kod =2:
IF(kodPrint ==2) THEN
det =(1.,0.)
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, '(t16, a)') 'A*x -b'
endif

x =b(1:n,kb)

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

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

ENDDO
ENDIF Check_if

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

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

do i=1,n
	do j=i+1,n
	u(i,j) =a(i,j)
	enddo
u(i,i) =1.
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 Full_if

deallocate(A, b, x, order)
if(allocated(a0)) deallocate(a0)
if(allocated(b0)) deallocate(b0)

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

END Program Crout_complex
