! Programul principal pentru metoda Cholesky
! Apeluri:
! call Cholesky(n,a)
! call forward(n,l,b)
! call back(n,lt,b)
!
! Datele se citesc din fisier de date:
! 1) Titlu = 80 chars.
! 2) n, nb		! n=nr. de ecuatii; nb= nr.termeni liberi
! 3) a: Triunghiul inferior.
! 4) b: Se scriu pe linii; "nb" linii a "n" termeni liberi.
! 5) kod		! cod de tiparire
! 6) check_text
!
! Check_text (20 char; Primul caracter e semnificativ):
! "C", "c", "P', "p" : se verifica solutia
!
! Stocare matrice A:
! A(i,j) -> a(k): k =(i-1)*i/2 +j =loca(i,j)

Program Cholesky_Method
real, allocatable ::a(:),b(:,:)
character(260) fname, titlu*80
character(4) nx, text*5
integer width
logical Print_It
common n

character(20) check_text
logical do_check
real, allocatable ::a0(:,:),b0(:,:)


print '(a,\)', ' Nume fisier de date: '
read '(a80)', fname
open(1, file=fname, status='old', form ='formatted')
! Citeste datele:
read(1,'(a80)') titlu
read(1,*) n, nb

allocate (a(n*(n+1)/2), b(n, nb))
! Citeste A (triunghi inferior):
do i=1,n
read(1,*) (a(loca(i,j)), j=1,i)
enddo

print '(/, a)', ' Matricea A:'
do i =1,n
write(*,*) (a(loca(i,j)), j=1,i)
enddo

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

read(1,*) kod0
print_it  =kod0 /=0

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

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

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

do i=1,n
do j=i+1,n
a0(i,j) =a0(j,i)
enddo
enddo

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

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

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

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

! Solutia:
if(nb ==1) then
write(1,'(/, a, /)') ' * Solutia pentru termenii liberi nr. 1'

else
iwnb =width(nb)
write(1,'(/, a, I<iwnb>, a)') ' * Solutia pentru termenii liberi nr. 1 - ', nb, ':'
write(1,'(t9, <nb>(i<iwnb>, 14x))') (k, k=1,nb)
endif

if(Print_It) then
write(*,'(/, a)') ' * Solutia:'
endif

do i=1,n
write(nx,'(I<width(i)>)') i
text =' x'//trim(nx)//':'
write(1, '(a, t6, <nb>g15.7)') trim(text),(b(i,k), k=1,nb)
	if(Print_It) write(*,*) (b(i,k), k=1,nb)
enddo

! Proba:
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
if(nb >1) write(1, '(/, a)') ': A*x -b'
endif

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

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

ENDDO
endif Check_if

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

END Program Cholesky_Method
