Program LU_D
! DUBLA Precizie
!
! Programul principal pentru LU (DOOLITLLE)
! Cu matricile: "A(n,n)" si "b(n,nb)"
! A. b: generate; citite din fisierul de intrare; citite din din fisier de date.
!
! Apel: LUdecomp; LUsolve; APVT.
!
! Selectia fisierului de intrare - prin fereastra de dialog standard.
!
! Fisier de date (intrare):
! 1) Titlu
! 2) n, nb, kod_Ab
! 3) [A - pentru kod_Ab =1]
! 4) [b - pentru kod_Ab =1]
! 5) kod_Print, kod_Inv
! 6) do_pvt, kod_Pivot
! 7) check_text
!
! Semnificatia:
! Titlu: text <= 80 caractere
!
! n, nb: nr. de ecuatii; nr. de termeni liberi.
! kod_Ab: Cod de intare a matricilor "A" si "b":
!		kodA_b <=0: matrici A,b, generate - de subrutina GenAB;
!		kod_Ab  =1: matrici citite din fisierul de intrare
!
!		kod_Ab  =2: matrice A generala, si b, citite dintr-un fisier de date;
!		kod_Ab >=3:matrice A simetrica, si b, citite dintr-un fisier de date
!				Prima linie in fisierul specificat e o linie text (max. 80 caractere).
!
! A: Matricea coeficientilor n x n
!    Se scrie asa cum e data - pe linii si coloane
!    Pentru kod_Ab >=3: se introduce triunghiul superior (pe linii).
!
! b: Termenii liberi
!    Se scriu pe linii: "nb" linii a "n" termeni liberi.
!
! kod_Inv: Cod_inversare_coeficienti
!	      kod_Inv: #0: a(i,j) -> 1./a(i,j)
! kod_Print: cod tiparire matrici
!	kod_Print >=1 : Se tiparesc: Matricile L si U.
!	kod_Print >=2: Se tipareste si vectorul "order", la fiecare pivotare.
!	               Se calculeaza si determinantul lui A.
!				   Daca |u(i,i,)| <prag: Se tipareste si linia "i" a matricii A (la pasul resoectiv).
!
! do_pvt: Cod pivotare (logic) 
!		do_pvt: T/F: Da/Nu - pivoteaza
! kodPivot: #0 = pivot maxim; 0 = primul pivot #0
!
! Check_text: text <= 20 caractere. Primul caracter e semnificativ:
!	'C', 'c', 'P', 'p', 'F', 'f': se face proba "Ax-b" (pentru fiecare termen liber).
!   C,c, P, p: se tipareste diferenta maxima in modul "Ax-b"
!   F,f: se tiparesc diferentele "Ax-b". Se verifica LU =PA
!
!
implicit real(8) (a-h, o-z)

real(8), allocatable :: A(:,:), b(:,:), x1(:)
integer, allocatable :: order(:)
character(260) fname, fname2, titlu*80
integer width

logical do_pvt
common no_interchange, kod_Print

real(8), allocatable :: a0(:,:), b0(:,:)

character(20) check_text;
logical do_check, full_check; 

real(8), allocatable :: l(:,:), u(:,:), c(:,:)

real(8) det
real(8) max_dif

character(80) dummy
character(100) title
integer status

! Selecteza fisierul de intrare:
title ='Select Input File'C
call GetFile(fname, status, title)

open(1, file =fname, status ='old', form ='formatted')

read(1,'(a)') TITLU

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

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

! citeste/genereaza A si b:
if (kod_Ab <=0) then
call GenAB(n, nb, a, b)

elseif (kod_Ab ==1) then
	do i =1,n
	read(1,*) (a(i,j), j= 1,n)
	enddo

	do j =1,nb
	read(1,*) (b(i,j), i= 1,n)
	enddo

else	! kod_Ab >=2
title ='Select Ab-Matrices File'C
	call GetFile(fname2, status, title)
	open(2, file =fname2, status ='old', form ='formatted')

	read(2,*) dummy
		sim_if: if(kod_Ab ==2) then
! A:
		do i =1,n
		read(2,*) (a(i,j), j= 1,n)
		enddo

		else	! matrice simetrica: triunghiul superior
		do i =1,n
		read(2,*) ( a(i,j), j=i,n )
		enddo

		! completeaza simetric:
		do i =1,n
		do j =1,i
		a(i,j) =a(j,i)
		enddo
		enddo

		endif sim_if

! b:
	do j =1,nb
	read(2,*) (b(i,j), i= 1,n)
	enddo

	close(2)

endif


! Citeste coduri:
read(1,*) kod_Print, kod_Inv

read(1, *) do_pvt, kod_Pivot

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

full_check =check_text(1:1) =='F' .or. check_text(1:1) =='f'

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

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

a =1./a
ENDIF

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


! Descompunera LU:
no_interchange =0

call LUdecomp(A, n, order, do_pvt, kod_Pivot, kodLU)
if(kodLU <0) then
	write(*,*) 'U-diagonal < Prag. Descompunerea LU nu e fezabila.'
	write(1, '(/, a)') ' U-diagonal < Prag. Descompunerea LU nu e fezabila.'

call EXIT
endif

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

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

Print_it: IF(kod_Print >=1) THEN
iwn =width(n)
diag =1.


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

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

ENDIF Print_it

! Calculeaza determinantul - pentru kod =2:
IF(kod_Print >=2) THEN
det =1.
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)') '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

	i_do: do i=1,n
	iwi =width(i)
	write(1, '(a, I<iwi>, a, G15.7)') ' x(', i, ') = ', b(i,kb)
	enddo i_do

enddo nb_do

! Proba: Solutia
Check_if: IF(do_check) THEN
write(1, '(//,a)') ' * Proba'

if(full_check) then
write(1, '(/,a)') ' Diferente (A*x-b)'

else
write(1, '(/,a)') ' Diferenta (A*x-b), maxima in modul'
endif


kb_do: DO kb=1,nb
iwk =width(kb)

if(nb >1) then
write(1, '(/, a, I<iwk>, a)') ' Solutia nr. ', kb, ':'
endif

!x =b(1:n,kb)

max_dif =0.
imax =0

ip_do: do i=1,n
sum =0.
	do k=1,n
	sum =sum +a0(i,k)*b(k,kb)		!x(k)
	enddo

dif =sum -b0(i,kb)
if(full_check) write(1,*) i, dif

!max_dif =max(abs(dif), max_dif)
adif =abs(dif)
if(adif >max_dif)then
max_dif =adif
dif_max =dif
imax =i
endif

enddo ip_do

if(full_check) then
write(1, '(/, a, \)') ' Diferenta (A*x-b), maxima in modul:'
endif

	if(imax >0) then
	iwm =width(imax)+1
	write(1, '(1pg15.7, a, i<iwm>)') dif_max, '; ecuatia', imax

	else
	write(1,*) dif_max
	endif


ENDDO kb_do


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-1
	l(i,j) =a(i,j)
	enddo
l(i,i) =1.
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(kod_Inv ==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, x1, 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 LU_D
