subroutine ell_eig2chol(n, u, lam, g)
      integer n
      double precision u(n, n), lam(n), g(n, n)
!  The n x n matrix A has the eigen decomposition A = u * lam^2 * u^T
!  and the Cholesky decomposition A = G * G^T.  Given u and lam, this
!  routine returns G.

!  Method:
!     A = u * lam^2 * u^T = B * B^T, where B = u * lam
!     B = G * Q  ( LQ factorization)

      integer j
      double precision B(n, n)

! form B
      do j = 1, n ! B = u * lam
        B(:,j) = u(:,j) * lam(j)
      end do

      call ell_bbt2chol(n, B, g)

      end