About MKL bidiagonal decomposition

[Zhonghan_L_]

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]);

[Zhonghan_L_]

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

[Zhen_Z_Intel]

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*P^H.  

Best regards,
Fiona

[Zhonghan_L_]

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*P^H.  

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

[Zhen_Z_Intel]

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

[Zhen_Z_Intel]

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 P^H maybe also different. It would be fine the results only with different sign. You could use ?ungbr to generate Q and P to check.