Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Zhonghan_L_
Beginner
63 Views

About MKL bidiagonal decomposition

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

 

0 Kudos
5 Replies
Zhonghan_L_
Beginner
63 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

Zhen_Z_Intel
Employee
63 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

Zhonghan_L_
Beginner
63 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

Zhen_Z_Intel
Employee
63 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

Zhen_Z_Intel
Employee
63 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. 

Reply