! Genereaza matricea R =S^(-T)*K*S^(-1)
! Cu Cholesky_S: M=S^(T)*S
!
SUBROUTINE gen_R(kod_Print)
use work

implicit real(8) (a-h,o-z)

real(8) YM(ngl*(ngl+1)/2)	!, S_1(ngl,ngl)
logical diagonal_m

real(8) e(ngl)


R =0.

! Check if Mass is diagonal:
diagonal_M =.true.

i_do: do i=1,ngl
j_do: do j=i+1,ngl

if(m(i,j) /=0) then
diagonal_m =.false.
exit i_do
endif

enddo j_do
enddo i_do

! Compute S =Mass^(1/2):
! 1) Diagonal:
Diag_IF: IF(diagonal_m) THEN

! Triunghiul sup.:
do i=1,ngl
do j=i,ngl
R(i,j) =K(i,j)/dsqrt( M(i,i)*M(j,j) )
enddo
enddo

! Simetric:
do i =2,ngl
do j =1,i
R(i,j) =R(j,i)
enddo
enddo

! Defineste S_1:
do i =1,ngl
S_1(i,i) =1./sqrt(M(i,i))
enddo

RETURN

ELSE

! 2) Non-Diagonal:
! Assign Mass to YM:
nn =ngl*(ngl+1)/2
YM =0.

! Store UPPER triangle - as required by Cholesky_S:
do i =1,ngl
do j =i,ngl	! upper triangle
YM(loca(i,j)) =M(i,j)
enddo

enddo


! Cholesky for Mass:
call Cholesky_S(ngl, YM, kod)		! YM =L
if(kod<0) RETURN

if(kod_Print /=0) then
write(1, '(/,a)') ' S:'
do i =1,ngl
write(1,*) (YM(loca(i,j)), j=i,ngl)
enddo
endif


! Now, YM stores the upper triangle S (of Mass).
! Compute S^(-1) by S*x =e, and Assign to S_1:
DO I =1,ngl
e =0.; e(I) =1.;
call Back(ngl, YM, e)

! Assign:
S_1(:,I) =e
ENDDO


! Compute R =S^(-T)*K*S^(-1):
R =matmul(K, S_1)
R =matmul(transpose(S_1), R)

ENDIF Diag_IF


END subroutine gen_R
