How to efficiently generate GSL matrices drawn from gaussian distributions?

Question

What is the most efficient way in which one can generate a gsl matrix with entries drawn from Gaussian distributions?

I know about the function called gsl_ran_gaussian but, as far as I know, it only generates scalar values. I tried filling a matrix using that function in a double for loop, but it is very slow for my intended application.

I would like to generate 3xN matrices: All the entries of the 1st 2 rows should be filled with samples from the same zero-mean Gaussian distribution of variance V1, while the entries of the last row should be filled with samples from a zero-mean Gaussian of variance V2.

EDIT:

This is how I am currently doing it:

gsl_matrix * M;
M = gsl_matrix_calloc(3, size2);
const gsl_rng_type * T;
gsl_rng * r;

gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc(T);

for (uint16_t k=0; k<size2; k++) {
    for (uint8_t h=0; h<2; h++)
        gsl_matrix_set(M, h, k, gsl_ran_gaussian(r,S1));
    gsl_matrix_set(M, 2, k, gsl_ran_gaussian(r,S2));
}