Program LU_2015
! 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.
!
! 2013: kod_PrintP >=3: Proba LU=A
!
! Apel: LUdecomp; LUsolve; APVT.
!
! Selectia fisierului de intrare - prin fereastra de dialog standard.
!
! Fisier de date (intrare/iesire):
! 1) Titlu
! 2) tip, kodInput
! 3) n, nb
! 4) [A] ... pentru kodA =1
! 5) [b] ... pentru kodA =1
! 6) kod_Inv, kod_PrintA, kod_PrintP
! 7) check_text
!
! Fisier extern (intrare):				
! 1) Titlu =linie text
! 2) Matricea A
! 3) vectorul/vectorii b

! Semnificatia:
! Titlu: text <= 80 caractere
!
! tip: text <= 10 caractere; primul e semnificativ: s, S =simetrica; alt caracter =generala;
! kodInput:     Matricie A si b:
!           0 = generate (de Gen_Ab);
!           1 = citite din fisier de intrare;
!           2 = citita din fiser extern (altul decat fisierul de intrare)
!           alta valoare = invalida

! n:  nr. de ecuatii;
! nb: nr. de termeni liberi.
!
! A: Matricea coeficientilor n x n
!    Tip "genarala": Se scrie asa cum e data - pe linii si coloane
!    Tip "simetrica": 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)
!
! kodPrintA: cod tiparire matrici
!	kod_PrintA: #0: Se tiparesc matricile A si b
!
! kod_PrintP: cod_Tiparire
!	kod_PrintP >=1 : Se tiparesc: Matricile L si U.
!	kod_PrintP >=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).
!	kod_PrintP >=3: Se face proba LU =A sau A-permutat.
!
! 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".
!
!* Pentru studiu pivotare, modifica valorile (linia 118):
! do_pvt = Cod pivotare (logic): T/F: Da/Nu - pivoteaza
! kodPivot: #0 = pivot maxim; 0 = primul pivot #0


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

logical do_pvt
common no_interchange, kod_PrintP

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

character(20) check_text;
logical do_check, full_check; 

real, allocatable :: xLU(:), row(:), col(:)

real(8) det
real max_dif

character(80) dummy
character(100) title
integer status
character(10) tip   ! gen/sim
logical simetric, general


! 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

! Matricea A:
read(1,*) tip, kodInput
! Verifica kodInput:
if(Input_Check(kodInput) <0) goto 510

read(1,*) n, nb
allocate (A(n,n), b(n, nb), x1(n))
allocate(order(n))
A =0.; b =0.; x1 =0;

simetric =tip(1:1) =='s' .or. tip(1:1) =='S'
general =.not. simetric

If_2: IF(kodInput ==2) THEN
	call GetFile(fname2, status, 'Select Ab-Matrices File'C)
	open(2, file =fname2, status ='old', form ='formatted')
ENDIF If_2

! Citeste/Genereaza:
call Read_Ab(n, nb, A, b, general, kodInput)

! Citeste coduri:
read(1,*) kod_Inv, kod_PrintA, kod_PrintP

do_pvt =.true.
kod_Pivot =1
!* Studiu: Fara pivotare:
!do_pvt =.false.
!* 
! 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

if(kod_PrintA /=0) then
! write A:
write(1,'(/,a)') 'Matricea A'
	do i =1,n
	write(1,*) (a(i,j), j= 1,n)
	enddo

! write b:
write(1,'(/,a)') 'Termenii liberi b'

	do i =1,n
	write(1,*) (b(i,j), j= 1,nb)
	enddo
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_PrintP >=2) then
	iinter =width(no_interchange)+1
	write(1, '(/, a, I<iinter>)') 'Numar de schimbari de linii =', no_interchange
	endif
enddo

LU_Print: IF(kod_PrintP >=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 LU_Print

! Calculeaza determinantul - pentru kod =2:
Det_Print: IF(kod_PrintP >=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, 1pE14.6)') 'Determinant(A) =', det
ENDIF Det_Print

! Proba LU=A
LU_A: IF(kod_PrintP >=3) THEN
if(no_interchange ==0) then
write(1, '(/,a)') '* Proba LU =A:'

else
write(1, '(/,a)') '* Proba LU =A-permutat:'
endif

if(kod_Inv ==0) then
write(1, '(a)') '  Matricea L*U'
else
write(1, '(a)') '  Matricea 1./L*U'
endif

allocate(xLU(n), row(n), col(n))
xLU =0.; row =0.; col =0.;

do i =1,n
if(i >=2) row(1:i-1) =a(i,1:i-1); 
row(i) =1.; row(i+1:n) =0;

xLU =0.
    do j=1,n
    col(1:j) =a(1:j,j); col(j+1:n) =0.;

    kij =min(i,j)
    xLU(j) =DOT_PRODUCT(row(:kij), col(:kij))
    enddo

if(kod_Inv ==0) then
write(1,*) xLU(:)    	
else
write(1,*) 1./xLU(:)    	
endif

enddo
ENDIF LU_A

! Scrie solutia:
iwnb =width(nb)
write(1, '(a)')

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

call Proba(n, nb, a0, b0, b, full_check)

ENDIF CHECK_IF


deallocate(A, b, x1, order)

if(allocated(a0)) deallocate(a0)
if(allocated(b0)) deallocate(b0)

if(allocated(xLU)) deallocate(xLU)
if(allocated(row)) deallocate(row)
if(allocated(col)) deallocate(col)

510 continue
END Program LU_2015
