I want to compute the SVD of two fairly small square matrices (of size 32x32 or 16x16). I use the following code to compute SVD:
int iRows = m1; //=32 int iCols = m2; //=16 int INFO; char JOBU = 'A'; char JOBVT = 'A'; float* A_Z1 = new float[iRows*iRows]; float* U_Z1 = new float[iRows*iRows]; float* VT_Z1 = new float[iRows*iRows]; float* S_Z1 = new float[iRows]; float* A_Z2 = new float[iCols*iCols]; float* U_Z2 = new float[iCols*iCols]; float* VT_Z2 = new float[iCols*iCols]; float* S_Z2 = new float[iCols]; int LWORK_Z1 = max(1,5*iRows); float* WORK_Z1 = new float[LWORK_Z1]; int LWORK_Z2 = max(1,5*iCols); float* WORK_Z2 = new float[LWORK_Z2]; //Compute SVD of vZ1a[a]................................................ for(int i = 0; i < iRows; ++i) for(int j = 0; j < iRows; ++j) A_Z1[j+i*iRows] = vZ1a[a](j,i); sgesvd(&JOBU, &JOBVT, &iRows, &iRows, A_Z1, &iRows, S_Z1, U_Z1, &iRows, VT_Z1, &iRows, WORK_Z1, &LWORK_Z1, &INFO); //Compute SVD of vZ2a[a]................................................ for(uint32_t i = 0; i < iCols; ++i) for(uint32_t j = 0; j < iCols; ++j) A_Z2[j+i*iCols] = vZ2a[a](j,i); sgesvd(&JOBU, &JOBVT, &iCols, &iCols, A_Z2, &iCols, S_Z2, U_Z2, &iCols, VT_Z2, &iCols, WORK_Z2, &LWORK_Z2, &INFO);
Later in my application I use left and right singular vectors for a machine learning algorithm and I have algorithms for measuring the accuracy of results.
I linked the same program with ATLAS+CLAPACK and MKL.
Surprisingly, I get very inaccurate results using MKL (i.e. PSNR using MKL is 53 and using CLAPACK is 73!).
Note that I am using the same code for both of them. All I do is to link to different libraries.
Whats wrong here? I don't see a difference between lapack function declarations (obviously there shouldn't be any). So my function call should be correct. Also I take care of the fact that I should pass the transpose of my input matrix since sgesvd needs column-major ordering
I compiled with both gcc and intel compiler but same results. I use the sequential libraries just to be safe.
I hope this is one of those stupid mistakes that I make because it has been bugging me for quite some time :)