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

issue with LAPACKE_dormqr

Azua_Garcia__Giovann
523 Views
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
0 Kudos
1 Solution
Ying_H_Intel
Employee
523 Views
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

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

View solution in original post

0 Kudos
7 Replies
mecej4
Honored Contributor III
523 Views
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.
0 Kudos
Azua_Garcia__Giovann
523 Views
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
0 Kudos
mecej4
Honored Contributor III
523 Views
> 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.
0 Kudos
Azua_Garcia__Giovann
523 Views
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
0 Kudos
mecej4
Honored Contributor III
523 Views
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.
0 Kudos
Ying_H_Intel
Employee
524 Views
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

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
0 Kudos
Azua_Garcia__Giovann
523 Views
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
0 Kudos
Reply