topic issue with LAPACKE_dormqr in Intel® oneAPI Math Kernel Library

issue with LAPACKE_dormqr

Azua_Garcia__Giovann — Thu, 29 Mar 2012 12:25:39 GMT

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

issue with LAPACKE_dormqr

mecej4 — Thu, 29 Mar 2012 13:23:25 GMT

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.

issue with LAPACKE_dormqr

Azua_Garcia__Giovann — Thu, 29 Mar 2012 14:53:35 GMT

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

issue with LAPACKE_dormqr

mecej4 — Thu, 29 Mar 2012 15:58:00 GMT

> 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 Q₁ and Fortran/C gives you Q₂ = -Q₁, Q₁'C will be equal to -Q₂'C, since you are supplying the C after the decomposition was performed.

issue with LAPACKE_dormqr

Azua_Garcia__Giovann — Thu, 29 Mar 2012 16:07:28 GMT

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

issue with LAPACKE_dormqr

mecej4 — Thu, 29 Mar 2012 16:15:29 GMT

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.

issue with LAPACKE_dormqr

Ying_H_Intel — Thu, 07 Jun 2012 08:07:35 GMT

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

issue with LAPACKE_dormqr

Azua_Garcia__Giovann — Fri, 08 Jun 2012 18:12:00 GMT

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