! 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):
! 1) Titlu = 80 chars.
! 2) n, nb, kod_Ab		
! 3) [a]
! 4) [b]
! 5) kod_Print 
! 6) check_text
!
! n =nr.ecuatii;
! nb =nr.termeni liberi
! kod_Ab =cod pentru matricile A si b.
!	kod_Ab <=0: matrici A, b = generate de Gen_Ab
!	kod_Ab  =1: matrici citite din fisiserul de intrare
!	kod_Ab >=2: matrici citite din fisiser de date specificat
!
! a: semi-banda superioara (inclusiv elem. diagonal).
!    Numai pt. kod_Ab =1
!
! b: Termenii liberi. Se scriu pe linii; "nb" linii a "n" termeni liberi.
!    Numai pt. kod_Ab =1
!
! kod_Print = Cod de tiparire
!	kod_Print =1: se tiparesc matricile A si b.
!   kod_Print >=2: se tipareste matricea S.
!
! Check_text (20 char; Primul caracter e semnificativ):
!	"C", "c", "P', "p", "F". "f" : se verifica solutia.
!	"F". "f" : se tipareste si proba A =ST*S
!
! 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_Ab, print_S

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,*) n, nb, kod_Ab

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

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

elseif (kod_Ab ==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_Print
print_Ab  =kod_Print ==1
print_S = kod_Print >=2

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


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

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

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

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

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


! 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
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


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

max_dif =0
imax=0
dif_max =0

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

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 b_DO



!* Check ST*S =A:
if(full_check) 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
!*
endif Check_if

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

END Program Cholesky_Method
