Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!

proper use of dgesvd?

aicoder
Beginner
126 Views
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....

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

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

0 Kudos
2 Replies
mecej4
Black Belt
126 Views
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
Moderator
126 Views
that's right.
Benson - please see C based DGESVD examples.
--Gennady
Reply