- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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]);
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page