! Solve3 + dubla precizie: sum
subroutine SOLVE(Y,b)
common n, LIM, NY0
real Y(NY0+n*LIM), b(n)
!real w(n)		! Tablou de lucru
real(8) w(n)
real(8) sum, si

!1) Forward:
! Calculeaza y din: ST*y =b. Prin substitutie inainte.
! Matricea S - data de Cholesky = stocata in 'a'
! 'y' stocat in 'b'

w=0.
b(1) =b(1)/Y(NY0+1)

F_DO: do i =2, n
si =Y(loca(i,i))
k1 =i-LIM+1
if(k1 <=0) then
k1=1
endif

	kf_do: do k =k1,i-1			! k =coloana
	w(k) =Y(loca(k,i))		! Conform ST: s(i,k)-> s(k,i)
	enddo kf_do

sum =dot_product(w(k1:i-1),b(k1:i-1))
yi =(b(i) -sum)/si
b(i) =yi

ENDDO F_DO

!2) Back:
! Calculeaza x din: S*x =y. Prin substitutie inapoi
! Matricea L - data de Cholesky = stocata in 'a'
! Vectorul 'y'->'b' dat de Forward
! 'x' stocat in 'b'

b(n) =b(n)/Y(loca(n,n))

B_DO: do i =n-1,1,-1
si =Y(loca(i,i))
k1 =i+LIM-1
if(k1 >n) then
k1=n
endif

	kb_do: do k =k1, i+1, -1	! k = coloana
	w(k) =Y(loca(i,k))
	enddo kb_do

sum =dot_product(w(i+1:k1),b(i+1:k1))
xi =(b(i) -sum)/si
b(i) =xi

ENDDO B_DO

END subroutine SOLVE


