function [rndmtout] = rndmt ( r, c, S, P, v )
// generate a r by c matrix of multivariate t variates
// rows are inedependent
// P is a matrix s.t. P'P = M^-1, c = dim(P)
// P can be obtained as chol(invpd(M)) (this is stupid Gauss should have a routine
// to invert triangular matrices)
// Set S = v and P = I for a standard multivariate t
z = grand(r,c,'normal')*( P * sqrt( S/v ) );
chi = grand(r,1,'gam',v/2,2/v); // not chi^2(v), v*chi is
rndmtout=z./sqrt(chi);
endfunction