- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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....

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

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/...

says:

- m INTEGER. The number of rows of the matrix
`A`().`m`0 - n INTEGER. The number of columns in
`A`().`n`0 -

- n INTEGER. The number of columns in
- 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

Link Copied

2 Replies

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page