topic Quote:Chunyu W. wrote:I think in Intel® oneAPI Math Kernel Library

LAPACKE_zgesvd has a bug

Chunyu_W_ — Sat, 06 Sep 2014 20:59:03 GMT

I think LAPACKE_zgesvd has a bug.

When I use this function to calculate the SVD of a 2*1 matrix [1;3], the result seems to be wrong.

Thanks.

Please show us the result,

mecej4 — Sat, 06 Sep 2014 21:41:27 GMT

Please show us the result, and tell us why it "seems to be wrong".

I tried the same 2 X 1 matrix example with MKL 11.1 and 11.2, and the results agreed with those from Matlab.

Given that your matrix is real, why do you call the complex SVD routine?

Thank you for your reply!

Chunyu_W_ — Sun, 07 Sep 2014 14:25:04 GMT

Thank you for your reply!

Do you use LAPACKE_zgesvd to calculate?

So you mean with different data type, we must use corresponding function. I can not call complex SVD routine when I use real matrix.

Oh, I think that's the problem. I always treat the real number as a special case of complex number and I just let the imaginary part to be 0 to call the complex routine.

Quote:Chunyu W. wrote:Do you

mecej4 — Sun, 07 Sep 2014 18:29:41 GMT

Chunyu W. wrote:
Do you use LAPACKE_zgesvd to calculate?

Yes. I used a modified version of the file lapacke_zgesvd_row.c that is distributed in the examples subdirectory of MKL.

So you mean with different data type, we must use corresponding function. I can not call complex SVD routine when I use real matrix.

You should know the properties of the three matrices that the SVD produces. If the original matrix is real, the U, S and V^T matrices are all real. When the routines to handle this case (sgesvd and dgesvd) are provided, why would you resort to using the more general routine for complex matrices? All the imaginary parts in the results will be zero, and the calculation will be much slower, without any purpose.

I still think that complex

Chunyu_W_ — Sun, 07 Sep 2014 23:19:47 GMT

I still think that complex SVD routine can calculate real matrix and it should be compatible with real matrix.

I use LAPACKE_zgesvd because sometimes I get complex matrix and sometimes I get real matrix, so I call this function overall, you know it is troublesome to judge what matrix I get.

I think the problem might have something to do with the superb parameter.

Below is an explanation about superb from the documentation(https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/)

On exit, superb(0:min(m,n)-2) contains the unconverged superdiagonal elements of an upper bidiagonal matrix B whose diagonal is in s (not necessarily sorted). B satisfies A = u*B*vt, so it has the same singular values as A, and singular vectors related by u and vt.

So my question is that when min(m,n) = 1, what should be the size of superb? This is just the case when I caculate [1;3].

Quote:Chunyu W. wrote:I think

mecej4 — Mon, 08 Sep 2014 01:25:14 GMT

Chunyu W. wrote:
I think the problem might have something to do with the superb parameter.

Since you have not yet stated what exactly the problem is, I cannot comment on that. The matrix B is of the same size as S. However, S is diagonal whereas B is upper bi-diagonal.

In your example, S is 1 X 1, so the size of the first upper diagonal is zero. However, C does not allow you to declare an array of size zero, so you could declare "double superb[max(1,min(M,N)-1)]. Alternatively, pass NULL as the last argument to LAPACKE_zgesvd.

Thank you for your patience!

Chunyu_W_ — Mon, 08 Sep 2014 02:40:11 GMT

Thank you for your patience!

I might paste my code here as well.

/**********************code*******************************/

void svd(int m, int n, FemasDComplex * a, FemasDComplex * u, FemasDComplex * vt, double * s, double * superb)
{
(void)LAPACKE_zgesvd(LAPACK_ROW_MAJOR, 'A', 'A', m, n, (MKL_Complex16 *)a, m, s, (MKL_Complex16 *)u, m, (MKL_Complex16 *)vt, n, superb);

}

int _tmain(int argc, _TCHAR* argv[])
{

FemasDComplex a[2], u[4], vt[1];

  double s[1];
   a[0].real = 1;
   a[0].imag = 0;
   a[1].real = 3;
   a[1].imag = 0;

for(int i=0; i<2; i++)
cout<<a.real<<'+'<<a.imag<<"i\n";
cout<<endl;

svd(2, 1, a, u, vt, s, NULL);

//print the result

cout<<s[0]<<endl;

   for(int i=0; i<4; i++)
       cout<<u.real<<'+'<<u.imag<<"i\n";
   cout<<endl;
   for(int i=0; i<1; i++)
       cout<<vt.real<<'+'<<vt.imag<<"i\n";
   cout<<endl;

while(1);

return 0;

}

/****************************code*******************************************/

above is c++ result

above is matlab result.

FemasDComplex is a structure

Chunyu_W_ — Mon, 08 Sep 2014 02:45:00 GMT

FemasDComplex is a structure the same as MKL_Complex16.

typedef struct FemasDComplex {
double real;
double imag;
} FemasDComplex;

I test some other matrix, they give the same result as matlab, just this one show a different result. [1 ; 3]

Since you are using the

mecej4 — Mon, 08 Sep 2014 08:36:44 GMT

Since you are using the LAPACK_ROW_MAJOR convention, the value passed for ldA should be n, not m. Many C users make this mistake and I think that the imprecise description in the MKL documentation contributes by using the term "leading dimension" in an inappropriate way.

Thank you. That's exactly the

Chunyu_W_ — Mon, 08 Sep 2014 23:17:04 GMT

Thank you. That's exactly the reason.

I also notice that svd is not unique actually, so it is sometimes possible that the result may not be the same as matlab.