Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Beginner
515 Views

Dear support team,

I met some problem when trying to run the ?unmbr() function from MKL, the bidiagonal terms could be calculated correctly, but I was unable to mutiply them back. Could you help me with it?

Thanks a lot!

```float x_real[6 * 4] = { -0.57, -1.28, -0.39,  0.25,
-1.93,  1.08, -0.31, -2.14,
2.30,  0.24,  0.40, -0.35,
1.93,  0.64, -0.66,  0.08,
0.15,  0.30,  0.15, -2.13,
0.02,  1.03, -1.43,  0.50};

int length[2] = { 4, 6 };

P = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*6*6, 32;
Q = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*6*6, 32);
TP = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8*6, 32);
TQ = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8*6, 32);

for (int j = 0; j < length[0] * length[1]; j++)
{
Spec.real = x_real;
Spec.imag = 0;
}

LAPACKE_cgebrd(LAPACK_ROW_MAJOR, length[1], length[0], Spec, length[0], s1, s2, TP, TQ); //here is correct

for (j = 0; j < length[1] * length[0] - 1; j++)
{
P.real = Spec.real;
P.imag = Spec.imag;
Q.real = Spec.real;
Q.imag = Spec.imag;

}

LAPACKE_cunmbr(LAPACK_ROW_MAJOR, 'Q', 'L', 'N', length[1], length[0], length[0], Q, length[0], TQ, S, length[1]);

LAPACKE_cunmbr(LAPACK_ROW_MAJOR, 'P', 'R', 'C', length[1], length[0], length[1], P, length[0], TP, S, length[1]);```

5 Replies
Beginner
515 Views

When I use

( 0.96,-0.81) (-0.03, 0.96) (-0.91, 2.06) (-0.05, 0.41)
(-0.98, 1.98) (-1.20, 0.19) (-0.66, 0.42) (-0.81, 0.56)
( 0.62,-0.46) ( 1.01, 0.02) ( 0.63,-0.17) (-1.11, 0.60)
(-0.37, 0.38) ( 0.19,-0.54) (-0.98,-0.36) ( 0.22,-0.20)
( 0.83, 0.51) ( 0.20, 0.01) (-0.17,-0.46) ( 1.47, 1.59)
( 1.08,-0.28) ( 0.20,-0.12) (-0.07, 1.23) ( 0.26, 0.26)

as input, the output bidiagonal matrix

-3.08701   2.11257   0                 0
0               2.06604   1.26281      0
0               0               1.87313     -1.61263
0               0                0                2.00218

is also incorrect, the correct output should be

Diagonal
-3.0870 -2.0660 -1.8731 -2.0022
Super-diagonal
2.1126 -1.2628 1.6126

Employee
515 Views

Hi Zhonghan,

I am afraid your code implementation is wrong. The bidiagonal matrix B just return as float array to s1 & s2 which s1 store diagonal elements, and s2 store super-diagonal. However, the ?unmbr is to calculate complex matrix C multiply with Q or P. I did not find any place you setting s1 &s2 to the real part of complex matrix S to computer S=S*Q, S=S*PH.

Best regards,
Fiona

Beginner
515 Views

Fiona Z. (Intel) wrote:

Hi Zhonghan,

I am afraid your code implementation is wrong. The bidiagonal matrix B just return as float array to s1 & s2 which s1 store diagonal elements, and s2 store super-diagonal. However, the ?unmbr is to calculate complex matrix C multiply with Q or P. I did not find any place you setting s1 &s2 to the real part of complex matrix S to computer S=S*Q, S=S*PH.

Best regards,
Fiona

Hi Fiona,

Thanks for your response! I didn't attach whole part of the code, I did set s1 and s2, sorry for the confuse. As you said, the s1 and s2 consist the bidiagonal term, and in my another comment it also unable to yield the correct bidiagonal result when comes with the complex number. Also there was a little typo in my first post(which won't affect the output bidiagonal terms), the order of TP and TQ is mistaken, should be

LAPACKE_cgebrd(LAPACK_ROW_MAJOR, length[1], length[0], Spec, length[0], s1, s2, TQ, TP);

Best regards,

Zhonghan

Employee
515 Views

Hi Zhonghan,

I do believe the result of MKL function is correct. Please follow below program:

```#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>

int main()
{
float x_real[6 * 4] = {
0.96, -0.03, -0.91, -0.05,
-0.98, -1.20, -0.66, -0.81,
0.62, 1.01, 0.63, -1.11,
-0.37, 0.19, -0.98, 0.22,
0.83, 0.20,-0.17, 1.47,
1.08, 0.20, -0.07, 0.26};
float x_imag[6*4] = {
-0.81, 0.96, 2.06, 0.41,
1.98, 0.19, 0.42, 0.56,
-0.46, 0.02,-0.17, 0.60,
0.38,-0.54,-0.36,-0.20,
0.51, 0.01,-0.46, 1.59,
-0.28,-0.12, 1.23, 0.26};

int length[2] = { 4, 6 };

MKL_Complex8* P = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*4*4, 32);
MKL_Complex8* Q = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*6*6, 32);
MKL_Complex8* TP = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*4, 32);
MKL_Complex8* TQ = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*4, 32);
MKL_Complex8* Spec = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*4*6,32);
MKL_Complex8* S = (MKL_Complex8*)MKL_malloc(sizeof(MKL_Complex8)*6*4,32);

float s1[4];
float s2[3];

for (int j = 0; j < length[0] * length[1]; j++)
{
Spec.real = x_real;
Spec.imag = x_imag;
}

LAPACKE_cgebrd(LAPACK_ROW_MAJOR, length[1], length[0], Spec, length[0], s1, s2, TQ, TP); //here is correct

printf("Matrix C:\n");
for(int i=0;i<length[1];i++)
{
for(int j=0;j<length[0];j++)
{
S[i*length[0]+j].real=0;
S[i*length[0]+j].imag=0;
if(i==j)
{
S[i*length[0]+j].real=s1;
}
else if(j-i==1)
{
S[i*length[0]+j].real=s2[j-1];
}
printf("%f+%fi\t",S[i*length[0]+j].real,S[i*length[0]+j].imag);
}
printf("\n");
}

LAPACKE_cunmbr(LAPACK_ROW_MAJOR, 'Q', 'L', 'N', length[1], length[0], length[0], Spec, length[0], TQ, S, length[0]);

LAPACKE_cunmbr(LAPACK_ROW_MAJOR, 'P', 'R', 'C', length[1], length[0], length[1], Spec, length[0], TP, S, length[0]);
printf("real part of Q*C*P':\n");
for(int i=0;i<length[1];i++)
{
for(int k=0;k<length[0];k++)
{
printf("%f\t", S[i*length[0]+k].real);
}
printf("\n");
}
printf("\nimag part of Q*C*P':\n");
for(int i=0;i<length[1];i++)
{
for(int k=0;k<length[0];k++)
{
printf("%f\t", S[i*length[0]+k].imag);
}
printf("\n");
}

return 0;
}```

The matrix A can get by multiplied back by using ?unmber function.

Best regards,
Fiona

Employee
515 Views

Another point is about the result of bidiagonal matrix B. The result of factorization would not be the only. For example, 1 can be factorized as 1*1, also can be -1*-1, it is depends on the algorithm of factorization. If you get the result with different sign, that means the elements in Q or PH maybe also different. It would be fine the results only with different sign. You could use ?ungbr to generate Q and P to check.