- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I have a matrix A=QR and the QR decomposition is computed using LAPACKE_dgeqrf. To be more concrete A and the output QR are included below. At this point I append a new column to A and it becomes A1, A1 is included below as well. At this point, I try to apply an update which requires computing first Q'*A1 and therefore I use the function LAPACKE_dormqr. This below is the documentation and the way I am invoking it for my given matrices:
matrix A of dimensions 4 x 2:
===================
7.500000e+01 4.000000e+00
4.000000e+01 2.800000e+01
6.600000e+01 5.000000e+00
1.800000e+01 1.000000e+01
matrix QR of dimensions 4 x 2:
===================
-1.091100e+02 -1.768857e+01
2.172614e-01 -2.474095e+01
3.584813e-01 -5.777660e-02
9.776762e-02 1.640589e-01
matrix A1 of dimensions 4 x 3:
===================
7.500000e+01 4.000000e+00 3.900000e+01
4.000000e+01 2.800000e+01 7.700000e+01
6.600000e+01 5.000000e+00 8.000000e+01
1.800000e+01 1.000000e+01 1.900000e+01
mklman.pdf MKL 10.3 Update 9:
===================
C:
lapack_int LAPACKE_ormqr( int matrix_order, char side, char trans, lapack_int m,
lapack_int n, lapack_int k, const* a, lapack_int lda, const *
tau,* c, lapack_int ldc );
int l_order = LAPACK_COL_MAJOR;
lapack_int m = 4; // C->rows;
lapack_int n = 3; // C->cols;
lapack_int k = 3; // B->QR->rows;
lapack_int lda = m;
lapack_int ldc = m;
lapack_int info = 0;
lapack_int info = LAPACKE_dormqr(l_order, 'L', 'T', m, n, k, B->QR->data, lda, B->tau->data, C->data, ldc);
So far so good, but now the matrix coming out Q'*A1 is:
matrix Q'*A1 of dimensions 4 x 3
===================
-1.091100e+02 -1.768857e+01 -1.065621e+02
0.000000e+00 -2.474095e+01 -4.110871e+01
0.000000e+00 0.000000e+00 3.281543e+01
0.000000e+00 8.881784e-16 -9.419678e+00
Cross checking with matlab turns out that the elements in column Q'*A1(:,3) are incorrect ... what am I doing wrong here?
This is the same but done using matlab:
>> A = [7.500000e+01 4.000000e+00;
4.000000e+01 2.800000e+01;
6.600000e+01 5.000000e+00;
1.800000e+01 1.000000e+01]
A =
75 4
40 28
66 5
18 10
>> A1 = [A [39 77 80 19]']
A1 =
75 4 39
40 28 77
66 5 80
18 10 19
>> [Q R] = qr(A)
Q =
0.6874 -0.3298 -0.5647 0.3161
0.3666 0.8696 -0.2063 -0.2585
0.6049 -0.2304 0.6340 -0.4232
0.1650 0.2862 0.4865 0.8088
R =
109.1100 17.6886
0 24.7410
0 0
0 0
>> Q'*A1
ans =
109.1100 17.6886 106.5621
0.0000 24.7410 41.1087
0.0000 0.0000 22.0579
-0.0000 -0.0000 -26.0582
The final matlab answer and the final MKL LAPACK answer are different in the last column:
>> Matlab=Q'*A1
Matlab =
109.1100 17.6886 106.5621
0.0000 24.7410 41.1087
0.0000 0.0000 22.0579
-0.0000 -0.0000 -26.0582
>> MKL=[
-1.091100e+02 -1.768857e+01 -1.065621e+02;
0.000000e+00 -2.474095e+01 -4.110871e+01;
0.000000e+00 0.000000e+00 3.281543e+01;
0.000000e+00 8.881784e-16 -9.419678e+00]
MKL =
-109.1100 -17.6886 -106.5621
0 -24.7410 -41.1087
0 0 32.8154
0 0.0000 -9.4197
>> abs(Matlab) - abs(MKL)
ans =
0.0000 -0.0000 0.0000
0.0000 0.0000 -0.0000
0.0000 0.0000 -10.7575
0.0000 0.0000 16.6385
Why is it so different, what am I doing wrong when calling LAPACKE_dormqr?
Thanks in advance,
Best regards,
Giovanni
I have a matrix A=QR and the QR decomposition is computed using LAPACKE_dgeqrf. To be more concrete A and the output QR are included below. At this point I append a new column to A and it becomes A1, A1 is included below as well. At this point, I try to apply an update which requires computing first Q'*A1 and therefore I use the function LAPACKE_dormqr. This below is the documentation and the way I am invoking it for my given matrices:
matrix A of dimensions 4 x 2:
===================
7.500000e+01 4.000000e+00
4.000000e+01 2.800000e+01
6.600000e+01 5.000000e+00
1.800000e+01 1.000000e+01
matrix QR of dimensions 4 x 2:
===================
-1.091100e+02 -1.768857e+01
2.172614e-01 -2.474095e+01
3.584813e-01 -5.777660e-02
9.776762e-02 1.640589e-01
matrix A1 of dimensions 4 x 3:
===================
7.500000e+01 4.000000e+00 3.900000e+01
4.000000e+01 2.800000e+01 7.700000e+01
6.600000e+01 5.000000e+00 8.000000e+01
1.800000e+01 1.000000e+01 1.900000e+01
mklman.pdf MKL 10.3 Update 9:
===================
C:
lapack_int LAPACKE_ormqr( int matrix_order, char side, char trans, lapack_int m,
lapack_int n, lapack_int k, const
tau,
int l_order = LAPACK_COL_MAJOR;
lapack_int m = 4; // C->rows;
lapack_int n = 3; // C->cols;
lapack_int k = 3; // B->QR->rows;
lapack_int lda = m;
lapack_int ldc = m;
lapack_int info = 0;
lapack_int info = LAPACKE_dormqr(l_order, 'L', 'T', m, n, k, B->QR->data, lda, B->tau->data, C->data, ldc);
So far so good, but now the matrix coming out Q'*A1 is:
matrix Q'*A1 of dimensions 4 x 3
===================
-1.091100e+02 -1.768857e+01 -1.065621e+02
0.000000e+00 -2.474095e+01 -4.110871e+01
0.000000e+00 0.000000e+00 3.281543e+01
0.000000e+00 8.881784e-16 -9.419678e+00
Cross checking with matlab turns out that the elements in column Q'*A1(:,3) are incorrect ... what am I doing wrong here?
This is the same but done using matlab:
>> A = [7.500000e+01 4.000000e+00;
4.000000e+01 2.800000e+01;
6.600000e+01 5.000000e+00;
1.800000e+01 1.000000e+01]
A =
75 4
40 28
66 5
18 10
>> A1 = [A [39 77 80 19]']
A1 =
75 4 39
40 28 77
66 5 80
18 10 19
>> [Q R] = qr(A)
Q =
0.6874 -0.3298 -0.5647 0.3161
0.3666 0.8696 -0.2063 -0.2585
0.6049 -0.2304 0.6340 -0.4232
0.1650 0.2862 0.4865 0.8088
R =
109.1100 17.6886
0 24.7410
0 0
0 0
>> Q'*A1
ans =
109.1100 17.6886 106.5621
0.0000 24.7410 41.1087
0.0000 0.0000 22.0579
-0.0000 -0.0000 -26.0582
The final matlab answer and the final MKL LAPACK answer are different in the last column:
>> Matlab=Q'*A1
Matlab =
109.1100 17.6886 106.5621
0.0000 24.7410 41.1087
0.0000 0.0000 22.0579
-0.0000 -0.0000 -26.0582
>> MKL=[
-1.091100e+02 -1.768857e+01 -1.065621e+02;
0.000000e+00 -2.474095e+01 -4.110871e+01;
0.000000e+00 0.000000e+00 3.281543e+01;
0.000000e+00 8.881784e-16 -9.419678e+00]
MKL =
-109.1100 -17.6886 -106.5621
0 -24.7410 -41.1087
0 0 32.8154
0 0.0000 -9.4197
>> abs(Matlab) - abs(MKL)
ans =
0.0000 -0.0000 0.0000
0.0000 0.0000 -0.0000
0.0000 0.0000 -10.7575
0.0000 0.0000 16.6385
Why is it so different, what am I doing wrong when calling LAPACKE_dormqr?
Thanks in advance,
Best regards,
Giovanni
1 Solution
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Giovanni,
Your parameter setting should be correct. But the different result between MKL and Matlab should be the QR implementation itself.
for example, the below showA=QR. Theelements in Greenaresignificant. But
Best Regards,
Ying
Your parameter setting should be correct. But the different result between MKL and Matlab should be the QR implementation itself.
for example, the below showA=QR. Theelements in Greenaresignificant. But
Red values could take any value with only one condition Q*Q^T=E.
As a result, MKL and Matlab have different last twocolumns, then different Q^T *A1 in last two row.

Best Regards,
Ying
Link Copied
7 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
A quick scan of your post showed me that you may be overlooking the fact that if A = Q R is the Q-R factorization, so also is (-Q) (-R), unless some additional restrictions are placed on R. Therefore, you should not compare only the Q-s or only the R-s given by Matlab and MKL. You should look at both.
Since Matlab itself uses MKL for internal operations, I would not expect to see any significant discrepancies.
Since Matlab itself uses MKL for internal operations, I would not expect to see any significant discrepancies.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello Mecej4,
I know what you mean but that's why I compare taking abs (see the end of my OP). It is also possible that the Q come different by permuting or changing signs for some row elements for a given column. But AFAIK Q'*C should give the same result specially like as you say matlab is just invoking MKL under the hood.
I guess this might have something to do with the way LAPACKE_dormqr works? Note that I am always using the compact QR representation upper trapezoidal R with lower min(M,N) householder reflectors and tau coeffients vector.
Best regards,
Giovanni
I know what you mean but that's why I compare taking abs (see the end of my OP). It is also possible that the Q come different by permuting or changing signs for some row elements for a given column. But AFAIK Q'*C should give the same result specially like as you say matlab is just invoking MKL under the hood.
I guess this might have something to do with the way LAPACKE_dormqr works? Note that I am always using the compact QR representation upper trapezoidal R with lower min(M,N) householder reflectors and tau coeffients vector.
Best regards,
Giovanni
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
> But AFAIK Q'*C should give the same result specially like as you say matlab > is just invoking MKL under the hood.
Two points:
1) The version of MKL used by Matlab need not be the same as the one used by your Fortran/C compiler.
2) If Matlab gives you Q1 and Fortran/C gives you Q2 = -Q1, Q1'C will be equal to -Q2'C, since you are supplying the C after the decomposition was performed.
Two points:
1) The version of MKL used by Matlab need not be the same as the one used by your Fortran/C compiler.
2) If Matlab gives you Q1 and Fortran/C gives you Q2 = -Q1, Q1'C will be equal to -Q2'C, since you are supplying the C after the decomposition was performed.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Mecej4,
Ok point taken :)
Could you confirm this is the right way to call the function?
QR is computed using DGEQRF over matrix B (m x (n - 1)).
C is a copy of matrix B with one extra column C (m x n).
C (m x n)
B->QR (m x (n - 1)) as returned by DGEQRF having MIN(m, n - 1) reflectors in the lower triangular.
int l_order = LAPACK_COL_MAJOR;
lapack_int m = C->rows;
lapack_int n = C->cols;
lapack_int k = MIN(B->QR->rows, B->QR->cols);
lapack_int lda = m;
lapack_int ldc = m;
lapack_int info = 0;
lapack_int info = LAPACKE_dormqr(l_order, 'L', 'T', m, n, k, B->QR->data, lda, B->tau->data, C->data, ldc);
Many thanks in advance,
Best regards,
Giovanni
Ok point taken :)
Could you confirm this is the right way to call the function?
QR is computed using DGEQRF over matrix B (m x (n - 1)).
C is a copy of matrix B with one extra column C (m x n).
C (m x n)
B->QR (m x (n - 1)) as returned by DGEQRF having MIN(m, n - 1) reflectors in the lower triangular.
int l_order = LAPACK_COL_MAJOR;
lapack_int m = C->rows;
lapack_int n = C->cols;
lapack_int k = MIN(B->QR->rows, B->QR->cols);
lapack_int lda = m;
lapack_int ldc = m;
lapack_int info = 0;
lapack_int info = LAPACKE_dormqr(l_order, 'L', 'T', m, n, k, B->QR->data, lda, B->tau->data, C->data, ldc);
Many thanks in advance,
Best regards,
Giovanni
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Rather than my attempting to play "compiler", it will be more effective if you post your code. Assuming that it is not more than a couple of 100 lines, I'll be happy to take a look at it.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Giovanni,
Your parameter setting should be correct. But the different result between MKL and Matlab should be the QR implementation itself.
for example, the below showA=QR. Theelements in Greenaresignificant. But
Best Regards,
Ying
Your parameter setting should be correct. But the different result between MKL and Matlab should be the QR implementation itself.
for example, the below showA=QR. Theelements in Greenaresignificant. But
Red values could take any value with only one condition Q*Q^T=E.
As a result, MKL and Matlab have different last twocolumns, then different Q^T *A1 in last two row.

Best Regards,
Ying
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Ying,
Thank you for following up on this :) I kinda of forgot this Thread already :) but is good to know. I got this bit working long ago.
Best!
Giovanni
Thank you for following up on this :) I kinda of forgot this Thread already :) but is good to know. I got this bit working long ago.
Best!
Giovanni

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