! Programul principal pentru metoda Cholesky
! Cu triunghiul superior al matricii A
!
! Apeluri:
! call Cholesky_S(n,a)
! call back(n,s,b)
! call forward(n,st,b)
!
! Datele se citesc din fisier de intrare (sau fisier de intrare + fisier specificat):
! Fisier de date (intrare/iesire):
! 1) Titlu
! 2) kodInput
! 3) n, nb
! 4) [A] ... pentru kodInput =1
! 5) [b] ... pentru kodInput =1
! 6) kod_Inv, kod_PrintA, kod_PrintS
! 7) check_text
!
! Fisier extern (intrare):				
! 1) Titlu =linie text
! 2) Matricea A
! 3) vectorul/vectorii b

! Semnificatia:
! Titlu: text <= 80 caractere
!
! 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 (nxn): 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_PrintA: Cod tiparire matrici. Se tiparesc:
!	kod_PrintA: >=1: Matricea A
!	kod_PrintA: >=2: termenii liberi b
!
! kod_PrintS: cod_Tiparire. Se tiparesc:
!   kod_PrintS: >=1: matricea S
!   kod_PrintS: >=2: determinantul.
!   kod_PrintS: >=3: proba A =ST*S.
!
! 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";
!
! Stocare matrice A:
! A(i,j) -> a(k): k =(j-1)*j/2 +i =loca(i,j)

Program Cholesky_Method

real, allocatable ::a(:),b(:,:)

character(260) fname, fname2, titlu*80
character(4) nx, text*5

integer width

logical print_A, print_b, print_S, print_Det, print_STS

character(20) check_text
logical do_check, full_check

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

character(80) dummy
real max_dif
integer status

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

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

! Citeste datele:
read(1,'(a80)') titlu
read(1,*) kodInput
! Verifica kodInput:
if(Input_Check(kodInput) <0) goto 510

read(1,*) n, nb

allocate (a(n*(n+1)/2), b(n, nb))

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

elseif (kodInput ==1) then
! Citeste A (triunghi superior):
do i =1,n
read(1,*) (a(loca(i,j)), j =i,n)
enddo

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

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

	read(2,*) dummy
	! A:
	do i =1,n
	read(2,*) (a(loca(i,j)), j =i,n)
	enddo

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

	close(2)

endif


read(1,*) kod_Inv, kod_PrintA, kod_PrintS

print_A =kod_PrintA >=1
print_b =kod_PrintA >=2

print_S   =kod_PrintS >=1
print_Det =kod_PrintS >=2
print_STS =kod_PrintS >=3

read(1, '(a)') 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

! Inversare elem. A:
if(kod_Inv /=0) then
a =1./a
endif

IF(do_check) THEN
allocate(a0(n,n), b0(n,nb))
b0 =b

do i =1,n
do j =i,n
lij =loca(i,j)
a0(i,j) =a(lij)
enddo
enddo

do i=2,n
do j=1,i-1
a0(i,j) =a0(j,i)
enddo
enddo
ENDIF

!*
Ab_IF: if(print_A .or. print_b) then
write(1, '(/, a)') ' Matricea A'
do i=1,n
write(1,*) a0(i,:)
enddo

if(print_b) then
write(1,'(/,a)') 'Termenii liberi b'

	do i =1,n
	write(1,*) (b(i,j), j= 1,nb)
	enddo
endif
endif Ab_IF
!*

call Cholesky_S(n, a, kod)
if(kod <0) then
GOTO 500
endif

write(1,*)
if(Print_S) then
! Tipareste matricea S:
write(1,'(a)') ' Matricea S'
do i =1,n
write(1,*) (a(loca(i,j)), j =i,n)
enddo
endif

do k=1,nb
call forward(n, a, b(:,k))
call back(n, a, b(:,k))
enddo

if(print_Det) then
det1 =1.
do i=1,n
ii =loca(i,i)
det1 =det1*a(ii)
enddo
det =det1*det1
write(1,'(/, a, 1pg14.7)') ' Determinant: ', det
endif

! Solutia:
iwnb =width(nb)

nb_do: do k =1,nb
write(1, '(/, a, I<iwnb>)') ' * Solutia pentru termenii liberi nr. ', k

	n_do: do i=1,n
	iwi =width(i)
	write(1, '(a, I<iwi>, a, G15.7)') ' x(', i, ') = ', B(i,k)
	enddo n_do
enddo nb_do

! Proba:
CHECK_IF: IF(do_check) THEN

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

ENDIF CHECK_IF


! Check ST*S =A:
if(print_STS) then
allocate(a1(n,n))
call compact_mul(n, a, a1)

write(1, '(/, a)') ' ST*S -A:'
do i=1,n
write(1,*) a1(i,:) -a0(i,:)
enddo

endif

500 deallocate(a,b)
if(allocated(a0)) deallocate(a0)
if(allocated(b0)) deallocate(b0)

510 continue
END Program Cholesky_Method
