subroutine ell_bbt2chol(n, B, gf)
      integer n
      double precision B(n, n), gf(n, n)
!  Given an n x n matrix B, compute the lower  
!  Cholesky triangle  G:   G * G^T = B * B^T

!  The Cholesky triangle G is returned in gf

!  Method: B = L * Q;  so that B * B^T = L * L^T; and so G = L.

      integer info, lwork
      double precision tau(n), work(10+n*(10+n)), vt(n, n)

!  perform LQ factorization: B = G*Q
      lwork = 10+n*(10+n)
      gf = B

      call dgelqf(n, n, vt, n tau, work, lwork, info)

      if( info /= 0 ) then
        write(0,*)'ell_bbt2chol: dgelqf failed, info = ', info
        stop
      endif

      end