subroutine qrdcmp(a,c,d,sing)
real(8), dimension(:,:), intent(inout) :: a
real(8), dimension(:), intent(out) :: c,d
logical, intent(out) :: sing
integer :: k,n
real(8) :: scale,sigma
n=assert_eq4(size(a,1),size(a,2),size(c),size(d),'qrdcmp')
sing=.false.
do k=1,n-1
scale=maxval(abs(a(k:n,k)))
if (scale == 0.0) then
sing=.true.
c(k)=0.0
d(k)=0.0
else
a(k:n,k)=a(k:n,k)/scale
sigma=sign(vabs(a(k:n,k)),a(k,k))
a(k,k)=a(k,k)+sigma
c(k)=sigma*a(k,k)
d(k)=-scale*sigma
a(k:n,k+1:n)=a(k:n,k+1:n)-outerprod(a(k:n,k),&
matmul(a(k:n,k),a(k:n,k+1:n)))/c(k)
end if
end do
d(n)=a(n,n)
if (d(n) == 0.0) sing=.true.
end subroutine qrdcmp