'Draw random samples from a matrix normal distribution.
Parameters
%(_matnorm_doc_default_callparams)s
size : integer, optional
Number of samples to draw (default 1).
%(_doc_random_state)s
Returns
rvs : ndarray or scalar
Random variates of size (`size`, `dims`), where `dims` is the
dimension of the random matrices.
Notes
%(_matnorm_doc_callparams_note)s'

def rvs(self, mean=None, rowcov=1, colcov=1, size=1, random_state=None):

    size = int(size)
    (dims, mean, rowcov, colcov) = self._process_parameters(mean, rowcov, colcov)
    rowchol = scipy.linalg.cholesky(rowcov, lower=True)
    colchol = scipy.linalg.cholesky(colcov, lower=True)
    random_state = self._get_random_state(random_state)
    std_norm = random_state.standard_normal(size=(dims[1], size, dims[0]))
    roll_rvs = np.tensordot(colchol, np.dot(std_norm, rowchol.T), 1)
    out = (np.rollaxis(roll_rvs.T, axis=1, start=0) + mean[np.newaxis, :, :])
    if (size == 1):
        out = out.reshape(mean.shape)
    return out
