I find that MKL and OpenBLAS give slightly different SVD decompositions.
I have this declaration of dgesvd
,
extern "C" {
void dgesvd_(char const& jobu, char const& jobvt, int const& mm, int const& nn, double* aa, int const& lda, double* s, double* u, int const& ldu, double* vt, int const& ldvt, double* work, int const& lwork, int& info); // NOLINT_INT(readability-identifier-naming)
}
This declaration works in principle with generic implementations of BLAS/LAPACK, such as OpenBLAS, and also with Intel MKL.
When I pass this array as input
0.5 1
2 2.5
I get for s
, vt
, and u
outputs, respectively as
3.38391 0.221637
-0.324536 0.945873
-0.945873 -0.324536
-0.606994 -0.794707
-0.794707 0.606994
When I link with generic LAPACK, I reconstruct the original array as vt^T.Diag(s).u
.
However when I link to MKL, (same code), I get this instead
3.38391 0.221637
0.324536 0.945873
0.945873 -0.324536
0.606994 0.794707
-0.794707 0.606994
I have to reconstruct the original array in a different way as vt.Diag(s).u
, that is there is an extra transposition.
So, it seems that these implementations give different kinds of SVD decompositions. Differing on the transposition on one of the array in the reconstruction.
I am reporting these arrays in C-ordering, (that is why V
appears on the left), but this is not the source of the problem.
Is this a bug in MKL? or am I missing something? Is this expected?
(This is very inconvinient because the program cannot be linked to generic BLAS and MKL.)