proper use of dgesvd?

[aicoder]

Hello,

I'm trying to compute an SVD, but I'm not sure if I'm using dgesvd properly. I copied the example code from here:

http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/dgesvd_ex.c.htm

and modified a few things to use the array from Wikipedia to try to verify the numbers (I verified the numbers in Wikipedia with Matlab):

http://en.wikipedia.org/wiki/Singular_value_decomposition#Example

Here is the code I used (I note below the code the changes I made):

#define M 5
#define N 4
#define LDA M
#define LDU M
#define LDVT N

...
...

int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
double wkopt;
double* work;
/* Local arrays */
double s, u[LDU*M], vt[LDVT*N];
double a[LDA*N] = {
1, 0, 0, 0, 2,
0, 0, 3, 0, 0,
0, 0, 0, 0, 0,
0, 4, 0, 0, 0
};
// Executable statements */
printf( " DGESVD Example Program Results\\n" );
/* Query and allocate the optimal workspace */
lwork = -1;
dgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork,
&info );
lwork = (int)wkopt;
work = (double*)malloc( lwork*sizeof(double) );
/* Compute SVD */
dgesvd( "All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,
&info );
/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm computing SVD failed to converge.\\n" );
exit( 1 );
}
/* Print singular values */
print_matrix( "Singular values", 1, n, s, 1 );
/* Print left singular vectors */
print_matrix( "Left singular vectors (stored columnwise)", m, m, u, ldu );
/* Print right singular vectors */
print_matrix( "Right singular vectors (stored rowwise)", n, n, vt, ldvt );
/* Free workspace */
free( (void*)work );

CHANGES:
* I changed the array to be the Wikipedia example
* I changed the print range of the "Left singular vectors" to be m, m (instead of m, n) [is that an error in the example?]
* I changed the #define to be M == 5 and N == 4

Here are my results:

Singular values
4.00 3.00 2.24 0.00

Left singular vectors (stored columnwise)
0.00 0.00 -0.45 0.00 -0.89
-1.00 0.00 -0.00 0.00 0.00
0.00 -1.00 -0.00 0.00 0.00
-0.00 -0.00 -0.00 1.00 0.00
0.00 -0.00 -0.89 -0.00 0.45

Right singular vectors (stored rowwise)
-0.00 -0.00 -0.00 -1.00
-0.00 -1.00 0.00 0.00
-1.00 -0.00 -0.00 -0.00
-0.00 -0.00 -1.00 0.00


The singular values are fine. The left and right singular vectors appear to be swapped from U and V on the Wikipedia entry. But here's my main problem: Notice the #define for M and N. The matrix on Wikipedia is a 4 x 5 matrix. The documentation for 'm' and 'n' from:

http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/cpp/win/mkl/refman/lse/functn_gesvd.html

says:

m INTEGER. The number of rows of the matrix A (m 0).

n INTEGER. The number of columns in A (n 0).

BUT

if I put the "right" values for M and N and change the #define to M == 4 and N == 5, then here are my results for the singular values:

Singular values
4.85 2.55 0.00 0.00 2.36

(which does not match Wikipedia)

What's going on?

Thanks!
Benson

[mecej4]

Please read the MKL documentation and observe that Lapack routines with dense matrix arguments (such as xGESVD) expect arrays of rank-2 (and higher, if needed) to be stored in column-major order (Fortran convention) instead of row-major order (C convention).

As a result of the mismatched storage order, you are passing M^T rather than M, and misinterpreting the matrices U and V^T after they are returned to you.

[Gennady_F_Intel]

that's right.

Benson - please see C based DGESVD examples.

--Gennady