linear algebra – GSL gsl_linalg_SV_decomp returning U and V reversed

I have the following C code:

``````void calculate_svd_example() {
int m = 4;
int n = 5;
double M(4)(5) = {{1.0, 0.0, 0.0, 0.0, 2.0},
{0.0, 0.0, 3.0, 0.0, 0.0},
{0.0, 0.0, 0.0, 0.0, 0.0},
{0.0, 2.0, 0.0, 0.0, 0.0}};

gsl_matrix *mat = gsl_matrix_alloc(m, n);
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
double x = M(i)(j);
gsl_matrix_set(mat, i, j, x);
}
}
printf("M = n");
pretty_print(mat);
run_svd(mat);
}
``````

``````#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <time.h>

#include "../include/run_svd.h"

/*
gsl_matrix_printf prints a matrix as a column vector.  This function
prints a matrix in block form.
*/
void pretty_print(const gsl_matrix * M)
{
// Get the dimension of the matrix.
int rows = M->size1;
int cols = M->size2;
// Now print out the data in a square format.
for(int i = 0; i < rows; i++){
for(int j = 0; j < cols; j++){
printf("%f ", gsl_matrix_get(M, i, j));
}
printf("n");
}
printf("n");
}

void pretty_print_vector(const gsl_vector * M)
{
int cols = M->size;
for(int j = 0; j < cols; j++){
printf("%f ", gsl_vector_get(M, j));
}
printf("n");
}

int run_svd(const gsl_matrix * a) {
gsl_matrix *aa;
int m = a->size1;
int n = a->size2;
if (m >= n) {
aa = gsl_matrix_alloc(m, n);
gsl_matrix_memcpy(aa, a);
} else {
aa = gsl_matrix_alloc(n, m);
// Need to transpose the input
gsl_matrix_transpose_memcpy(aa, a);
}

m = aa->size2;
gsl_matrix * V = gsl_matrix_alloc(m, m);
gsl_vector * S = gsl_vector_alloc(m);
gsl_vector * work = gsl_vector_alloc(m);

/**
* On output the matrix A is replaced by U. The diagonal elements of
* the singular value matrix S are stored in the vector S. The
* singular values are non-negative and form a non-increasing sequence
* from S_1 to S_N. The matrix V contains the elements of V in
* untransposed form. To form the product U S V^T it is necessary to
* take the transpose of V. A workspace of length N is required in
* work.
*/
gsl_linalg_SV_decomp(aa, V, S, work);
printf("U:n");
pretty_print(aa);
printf("S:n");
pretty_print_vector(S);
printf("V:n");
pretty_print(V);

gsl_matrix_free(V);
gsl_vector_free(S);
gsl_vector_free(work);

return 0;
}
``````

It is giving me the following results:

``````U:
-0.000000 -0.447214  0.000000 -0.000000
0.000000 -0.000000 -1.000000 -0.000000
-1.000000 -0.000000  0.000000  0.000000
-0.000000 -0.000000  0.000000  1.000000
-0.000000 -0.894427  0.000000  0.000000

S:
3.000000 2.236068 2.000000 0.000000

V:
-0.000000 -1.000000 -0.000000 0.000000
-1.000000 -0.000000 -0.000000 0.000000
-0.000000 -0.000000 -0.000000 1.000000
-0.000000 -0.000000 -1.000000 0.000000
``````

Aren’t `U` and `V` reversed here? Is this a problem with the GSL code, the GSL documentation, or something I’m doing wrong?