Help with QL factorization

[sguazt]

Hello,
I am a beginner LAPACK user.
I am playing with QL factorization.
I'd have problem to obtain the full Q matrix through the ORGQL subroutine.
Specifically, given an m-by-n matrix A, I would like to build the "full" m-by-m matrix Q such that A=QL, for any value of m and n (i.e., both when m >= n or n < m).
To date, when m > n, I was able to successfully obtain only the "reduced" m-by-n matrix Q, which corresponds to the n trailing columns of the "full" Q matrix.

Is there a way in LAPACK to obtain the "full" Q matrix?
Or simply I compute the "reduced" Q and set the m-n heading columns to zero (when m>n)?

As a test case, I'm using the following 6-by-4 input matrix A:

Code: Select all

  [ -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 ]


I perform the following steps:
1. Call DGEQLF with LWORK=-1 (in order to get the optimal workspace size)
DGEQLF(6, 4, A, 6, TAU, WORK, -1, INFO)

2. Then I call DGEQLF again but this time with workspace properly resized (LWORK=INFO):
DGEQLF(6, 4, A, 6, TAU, WORK, LWORK, INFO)
After this call I get the following result:

Code: Select all

A = [ -1.942940   1.506660   0.933778  -0.096466;
         0.899115   1.403620  -1.277710   0.825749;
         2.894850   0.455552  -1.038550   0.135052;
         0.504110   1.581270   1.283920  -0.030869;
         1.921280  -1.053210   1.692760   0.821891;
         0.872950  -0.901803  -0.213936   3.091590 ]

TAU = [ 0.35820;
             0.36713;
             0.32095;
             0.83827 ]


3. Then I call DORGQL in order to obtain the matrix Q (as before I call it two times, where the first one is for getting the optimal workspace size). Here is the problem

3.a. If I call: DORGQL(6, 4, min(6,4), A, 6, TAU, WORK, LWORK, INFO)
I get the trailing 4 columns of Q matrix

Code: Select all

Q = [ 0.0833107  -0.9100041  -0.2201730   0.0808647;
        -0.2971665   0.1079856  -0.2706159  -0.6922014;
         0.6403848   0.2350681   0.2219915  -0.1132102;
        -0.4461129   0.1619876  -0.3866238   0.0258767;
         0.2937612  -0.2021707   0.0015402  -0.6889674;
         0.4575411   0.1945615  -0.8243342   0.1617290 ]


4.a While if I call: DORGQL(6, 6, min(6,4), A, 6, TAU, WORK, LWORK, INFO)
I get a totally wrong matrix.

If I compute the Q matrix manually (i.e., by computing H(i), such that Q=H(4)*H(3)*H(2)*H(1)) I get:

Code: Select all

Q = [ -0.1280591  -0.3058246   0.0833107  -0.9100041  -0.2201730   0.0808647;
         -0.5697119  -0.1519170  -0.2971665   0.1079856  -0.2706159  -0.6922014;
         -0.0751564  -0.6833014   0.6403848   0.2350681   0.2219915  -0.1132102;
          0.5633922  -0.5542341  -0.4461129   0.1619876  -0.3866238   0.0258767;
          0.5762694   0.2570164   0.2937612  -0.2021707   0.0015402  -0.6889674;
         -0.0621795   0.2079659   0.4575411   0.1945615  -0.8243342   0.1617290 ]


which successfully satisfies A=QL.

Thank you very much for helping me!!

Best,

-- Marco

[Julien Langou]

Hi Marco,

All you are writing makes sense. The header of DORGQL explains how the routine is supposed to be called but, yep, it's not the same as DORGQR and so this leads to some misuse. A good hint on how to use the subroutine is to go in the TESTING directory and in your case, you want to look at: TESTING/dqlt01.f (for the square Q). (See TESTING/dqlt02.f for rectangular Q.) So the header says:

Code: Select all

*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the (n-k+i)-th column must contain the vector which
*          defines the elementary reflector H(i), for i = 1,2,...,k, as
*          returned by DGEQLF in the last k columns of its array
*          argument A.
*          On exit, the M-by-N matrix Q.


This means that on INPUT of DORGQL, " A " needs to be shifted on the left.

So: in the case, A is M-by-N with M > N, and you want to compute the full M-by-M Q (as opposed to the reduced M-b-N Q): you need to start your M-by-M array (that will hold Q) by leaving M-N irrelevant columns, then filling the last N relevant columns with the input matrix A, then you call DGEQLF on "A(1,M-N+1)" and call DORGQL on "A(1,1)". So you were almost there. See code below.

The LAPACK TESTING/dqlt01.f performs the "same". It starts with an M-by-N matrix A in "standard form" (so not shifted on the left by M-N+1), calls DGEQLF on A, then performs a "copy details of Q" which performs the shift and then calls DORGQL.

Please note that your "reduced Q factor" for the QL factorization is not correct as well. (I think, and I think you implied it was not correct in your email, but I was not sure.)

Julie and Julien.

Code: Select all

      PROGRAM QLFACT
     
      INTEGER LWORK
      DOUBLE PRECISION   ROGUE
      PARAMETER ( LWORK = 128 )
      PARAMETER ( ROGUE = -1.0D+10 )
      DOUBLE PRECISION A(6,6), TAU(4)
      DOUBLE PRECISION WORK(LWORK)
      INTEGER INFO, I, J
      EXTERNAL DGEQLF, DORGQL

      A(1,1) = ROGUE
      A(1,2) = ROGUE
      A(2,1) = ROGUE
      A(3,1) = ROGUE
      A(4,1) = ROGUE
      A(5,1) = ROGUE
      A(6,1) = ROGUE

      A(2,2) = ROGUE
      A(3,2) = ROGUE
      A(4,2) = ROGUE
      A(5,2) = ROGUE
      A(6,2) = ROGUE

      A(1,3) = -0.57
      A(2,3) = -1.93
      A(3,3) =  2.30
      A(4,3) = -1.93
      A(5,3) =  0.15
      A(6,3) = -0.02

      A(1,4) = -1.28
      A(2,4) =  1.08
      A(3,4) =  0.24
      A(4,4) =  0.64
      A(5,4) =  0.30
      A(6,4) =  1.03

      A(1,5) = -0.39
      A(2,5) = -0.31
      A(3,5) =  0.40
      A(4,5) = -0.66
      A(5,5) =  0.15
      A(6,5) = -1.43

      A(1,6) =  0.25
      A(2,6) = -2.14;
      A(3,6) = -0.35
      A(4,6) =  0.08
      A(5,6) = -2.13
      A(6,6) =  0.50
     
      DO I=1,6
         DO J=3,6
            WRITE(*,*) 'A(',I,',',J-2,')=',A(I,J),';'
         END DO
      END DO

      CALL DGEQLF(6, 4, A(1,3), 6, TAU, WORK, LWORK, INFO)
      CALL DORGQL(6, 6, min(6,4), A, 6, TAU, WORK, LWORK, INFO)
     
      DO I=1,6
         DO J=1,6
            WRITE(*,*) 'Q(',I,',',J,')=',A(I,J),';'
         END DO
      END DO

      RETURN
      END PROGRAM QLFACT


[sguazt]

Hi Marco,

All you are writing makes sense. The header of DORGQL explains how the routine is supposed to be called but, yep, it's not the same as DORGQR and so this leads to some misuse. A good hint on how to use the subroutine is to go in the TESTING directory and in your case, you want to look at: TESTING/dqlt01.f (for the square Q). (See TESTING/dqlt02.f for rectangular Q.)

Hi Julie & Julien,

Sorry for the late reply, but I was out of office for a few days.

I just want to thank you for your help!!

Please note that your "reduced Q factor" for the QL factorization is not correct as well. (I think, and I think you implied it was not correct in your email, but I was not sure.)

Please can you tell me what is it wrong?

I have rearranged your sample program qlfact in order to compute the "reduced" Q matrix (see qlfact2 program below) and the resulting Q matrix is the same as the one you can find in the first post.

Code: Select all

      PROGRAM QLFACT2

      INTEGER LWORK
      DOUBLE PRECISION   ROGUE
      PARAMETER ( LWORK = 128 )
      PARAMETER ( ROGUE = -1.0D+10 )
      DOUBLE PRECISION A(6,4), TAU(4)
      DOUBLE PRECISION WORK(LWORK)
      INTEGER INFO, I, J
      EXTERNAL DGEQLF, DORGQL

      A(1,1) = -0.57
      A(2,1) = -1.93
      A(3,1) =  2.30
      A(4,1) = -1.93
      A(5,1) =  0.15
      A(6,1) = -0.02

      A(1,2) = -1.28
      A(2,2) =  1.08
      A(3,2) =  0.24
      A(4,2) =  0.64
      A(5,2) =  0.30
      A(6,2) =  1.03

      A(1,3) = -0.39
      A(2,3) = -0.31
      A(3,3) =  0.40
      A(4,3) = -0.66
      A(5,3) =  0.15
      A(6,3) = -1.43

      A(1,4) =  0.25
      A(2,4) = -2.14;
      A(3,4) = -0.35
      A(4,4) =  0.08
      A(5,4) = -2.13
      A(6,4) =  0.50

      WRITE(*,*) 'A = '
      DO I=1,6
         WRITE(*,*) ( A(I,J), J = 3,6 )
      END DO

      CALL DGEQLF(6, 4, A, 6, TAU, WORK, LWORK, INFO)
      CALL DORGQL(6, 4, min(6,4), A, 6, TAU, WORK, LWORK, INFO)

      WRITE(*,*) 'Q = '
      DO I=1,6
         WRITE(*,*) ( A(I,J), J = 1,4 )
      END DO

      RETURN
      END PROGRAM QLFACT2


Thank you very much!!

Cheers,

-- Marco