! Programul principal pentru metoda Cholesky - Matrici BANDA simetrica
! Se lucreaza cu semi-banda SUPERIOARA.
! Apeluri:
! call Cholesky(a, kod); call forward(l,b); call back(lt,b)
!
! Datele se citesc din fisier de date:
! 1) Titlu = 80 chars.
! 2) n, LIM1, nb		
!    n =nr.ecuatii; LIM =semi-latime de banda (inclusiv elem. diagonal); nb =nr.termeni liberi
! 3) a: SEMI-BANDA SUPERIOARA (inclusiv elem. diagonal).
! 4) b: Se scriu pe linii; "nb" linii a "n" termeni liberi.
! 5) kod: cod de tiparire: 
!    kod =0/#0: Se tipareste matricea S; kod >=2: se tipareste proba ST*S =A
!
! 6) check_text
! Check_text (20 char; Primul caracter e semnificativ):
! "C", "c", "P', "p" : se verifica solutia
!
! Stocare matrice A:
! A -> a(n*LIM): A(i,j) -> a(k): k =(i-1)*LIM +j-i+1 = loca(i,j)
! Tabloul Y = dimensionat ca: F + a

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

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


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, LIM, nb

! Alocari: (a(n*LIM), b(n, nb))
NY0 = LIM*(LIM+1)/2
NY = NY0 +n*LIM
allocate (Y(NY))
allocate (b(n, nb))
Y(:) = 0

! Citeste A (banda superioara):
do i=1,n
j2 =min(n, i+LIM-1)
read(1,*) (Y(loca(i,j)), j =i, j2)		!i+LIM-1
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

a0 =0
do i=1,n
do j=i,n
!lij =loca(i,j)
a0(i,j) =AA(i, j, Y)
enddo
enddo

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

ENDIF

! Tipareste A:
print '(/, a)', ' Matricea A:'
write(1,'(/, a)') ' Matricea A:'
do i =1,n
j1 =i; j2 =min(i+LIM-1, n); nj =j2-j1+1
write(*, '(<nj>g15.7)') (Y(loca(i,j)), j= j1,j2)
write(1, '(<nj>g15.7)') (Y(loca(i,j)), j= j1,j2)
enddo

call Cholesky_for_Band(Y, kod)
if(kod <0) then
GOTO 500
endif

write(1,*)
IF(print_it) THEN	! Tipareste S:
print '(/, a)', ' Matricea S:'
write(1,'(/, a)') ' Matricea S:'

do i =1,n
j1 =i; j2 =min(i+LIM-1, n); nj =j2-j1+1
write(*, '(<nj>g15.7)') (Y(loca(i,j)), j= j1,j2)
write(1, '(<nj>g15.7)') (Y(loca(i,j)), j= j1,j2)

enddo

ENDIF

! Proba-1: ST*S =A
IF(print_it .and. kod0 >=2) THEN
write(1, '(/, a)') ' Proba A =ST*S:'

DO_i: DO i=1,n
	DO_j: do j=1,n
	aij =0.
			do k=1,n
			aij =aij +AA(k,i, Y)*AA(k,j, Y)
			enddo
	if(j <=n-1) then
	write(1,'(g15.7, \)') aij
	else
	write(1,'(g15.7)') aij
	endif
	ENDDO DO_j
ENDDO DO_i

ENDIF

do k=1,nb
call solve(Y, b(:,k))
enddo


! Scrie 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 2 - Solutia:
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.D0,0.D0)
	do k=1,n
	sum =sum +a0(i,k)*b(k,kb)		! ok
	enddo

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

ENDDO
endif Check_if

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

END Program Cholesky_Band
