! 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)

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(loca(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(loca(i,j))/dsqrt( M(loca(i,i))*M(loca(j,j)) )
enddo
enddo

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

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

RETURN

ELSE

! 2) Non-Diagonal:
! Assign Mass to YM:

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


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



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

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

deallocate(M)

ENDIF Diag_IF


! Compute R =S^(-T)*K*S^(-1):
!R =matmul(KK, S_1)
call compact_mul(ngl, K, S_1, R)

!*
!write(1,'(//, a)') ' Matricea R1 =K*S^(-1)'

!do i =1,ngl
!write(1, '(<ngl>(1pg24.16))') R(i,:)
!enddo
!*

R =matmul(transpose(S_1), R)

deallocate(K)

END subroutine gen_R
