! Descompunerea Cholesky: A=ST*S
! Matricea A = Banda simetrica + pozitiv definita
! A este generata de functia AA.
! Daca A este data numeric: AA citeste din tabloul 'a' ~ Y(NY0+1:NY).
! Matricea S in locatiile matricii A (stocare compacta)
! Cod scris cu tablouri; Cod unitar pentru: J=1,N
! Y(NY0+n*LIM)
! CholoeskyBand3 + Dubla precizie pentru "sum".

subroutine Cholesky_for_Band(Y, kod)
common n, LIM, NY0
real(4) Y(NY0 + n*LIM)
!real(4) w1(LIM-1), w2(LIM-1)	! Tablouri de lucru
real(8) w1(LIM-1), w2(LIM-1)
real(8) sum, si, s

integer width


kod =0
NY = NY0 + n*LIM

! PASII J =1,...,LIM
! Genereaza primul triunghi, pe linii (locale = generale):
DO i =1, LIM
	do j = i, LIM		! pana la a(i,LIM)
	ly = LOF(i,j)
	Y(ly) = AA(i,j, Y)
	enddo
ENDDO

! Test:
!print *, 'Y generat la pasul I'
!print *, Y(1:NY0)
!print *
!write(1,*)' Y generat la pasul I'
!write(1,*) Y(1:NY0)
!write(1,*)

! PROCESEAZA coloanele: J =1,..,LIM; LIM+1,...,N
! Element (1,1):
ly1 = LOF(1,1)
s1 = sqrt(Y(ly1))
Y(ly1) = s1

DO_J:	DO J = 2, N
if(J  <=LIM) then
iLIM =J
else
iLIM =LIM
endif

! Skip pentru J <=LIM:
if (J <=LIM) goto 100

! Atribuire: liniile 2-LIM -> 1-(LIM-1)
	do il =2, LIM
		do jl =il, LIM
		ly1 = LOF(il,jl)
		ly2 = LOF(il-1,jl-1)
		Y(ly2) =Y(ly1)
		enddo
	enddo

! Generare coloana LIM:
	do il =1, LIM
	ly =LOF(il, LIM)
	Y(ly) =AA(J-LIM+il, J, Y)
	enddo

100 continue
! Calculeaza coloana locala iLIM ~ coloana 'J'
! Primul element, linia il =1:
ly  = LOF(1,iLIM)		! jl = il +LIM -1
ly1 = LOF(1,1)			! elementul 1 din 'il'
Y(ly) =Y(ly)/Y(ly1)

! Elementele din liniile locale 2,...,LIM-1:
DO_il:	do il =2, iLIM-1
	lyi =LOF(il,il)
	si =Y(lyi)			! numitor = (il,il)

	ly =LOF(il,iLIM)		! element curent
	s =Y(ly)			! s =a(i,j)

! a(i,j): se inmultesc coloanele i, LIM - din liniile locale 1,...,(i-1)
!	k1 =1; 	k2 =il-1
		do k =1, il-1
		lyki =LOF(k,il)
		lykj =LOF(k,iLIM)
		w1(k) =Y(lyki)
		w2(k) =Y(lykj)
		enddo
	sum =dot_product(w1(1:il-1),w2(1:il-1))
	s =(s-sum)/si

	Y(ly) =s
Enddo DO_il

! Elementul diagonal, il=iLIM:
ly = LOF(iLIM,iLIM)
s = Y(ly)

!
w1 =0.

! sumeaza elementele supra-diagonale: coloana de 'LIM-1' elemente
	do k =1, iLIM-1
	lyk = LOF(k,iLIM)
	w1(k) =Y(lyk)
	enddo

sum =dot_product(w1,w1)
s =s -sum
if(s <=0) then
iwj =width(J)
print '(/, a)', ' * Matricea nu este pozitiv definita !'
print '(a, I<iwj>, a, g15.7 )', ' Linia ', J,':	DIAG =', s

write(1, '(/, a)') ' * Matricea nu este pozitiv definita !'
write(1, '(a, I<iwj>, a, g15.7 )') ' Linia ', J, ':	DIAG =', s

kod =-1
RETURN

else
Y(ly) =dsqrt(s)
endif

! Scrie linia procesata: pentru J =LIM,...,N-1, in 'a' (linia=1,...,N-LIM+1):
! La J = N (ultimul triunghi) scrie TOATE liniile:
IF(J>= LIM .and. J<N) THEN
I =J -LIM +1
ly0 =(I-1)*LIM +NY0
do jl =1, LIM
Y(ly0 +jl) =Y(jl)
enddo

ELSEIF (J==N) THEN
I0 =N -LIM		! J=N-1 => I0 =N-LIM
do il =1, LIM
I = I0 + il
ly0 =(I-1)*LIM	+NY0
	do ll =1, LIM
	jl = il +ll -1		! coloana locala
	if(jl >LIM) cycle	! Y e initializat cu 0.
	ly =LOF(il,jl)
	Y(ly0 +ll) =Y(ly)
	enddo
enddo
ENDIF

ENDDO DO_J

END subroutine Cholesky_for_Band
