LAPACK 3.3.1 Linear Algebra PACKage

# VARIANTS/qr/LL/dgeqrf.f

Go to the documentation of this file.
```00001       SUBROUTINE DGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     March 2008
00006 *
00007 *     .. Scalar Arguments ..
00008       INTEGER            INFO, LDA, LWORK, M, N
00009 *     ..
00010 *     .. Array Arguments ..
00011       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
00012 *     ..
00013 *
00014 *  Purpose
00015 *  =======
00016 *
00017 *  DGEQRF computes a QR factorization of a real M-by-N matrix A:
00018 *  A = Q * R.
00019 *
00020 *  This is the left-looking Level 3 BLAS version of the algorithm.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  M       (input) INTEGER
00026 *          The number of rows of the matrix A.  M >= 0.
00027 *
00028 *  N       (input) INTEGER
00029 *          The number of columns of the matrix A.  N >= 0.
00030 *
00031 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00032 *          On entry, the M-by-N matrix A.
00033 *          On exit, the elements on and above the diagonal of the array
00034 *          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
00035 *          upper triangular if m >= n); the elements below the diagonal,
00036 *          with the array TAU, represent the orthogonal matrix Q as a
00037 *          product of min(m,n) elementary reflectors (see Further
00038 *          Details).
00039 *
00040 *  LDA     (input) INTEGER
00041 *          The leading dimension of the array A.  LDA >= max(1,M).
00042 *
00043 *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
00044 *          The scalar factors of the elementary reflectors (see Further
00045 *          Details).
00046 *
00047 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00048 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00049 *
00050 *  LWORK   (input) INTEGER
00051 *
00052 *          The dimension of the array WORK. The dimension can be divided into three parts.
00053 *
00054 *          1) The part for the triangular factor T. If the very last T is not bigger
00055 *             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
00056 *             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
00057 *
00058 *          2) The part for the very last T when T is bigger than any of the rest T.
00059 *             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
00060 *             where K = min(M,N), NX is calculated by
00061 *                   NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
00062 *
00063 *          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
00064 *
00065 *          So LWORK = part1 + part2 + part3
00066 *
00067 *          If LWORK = -1, then a workspace query is assumed; the routine
00068 *          only calculates the optimal size of the WORK array, returns
00069 *          this value as the first entry of the WORK array, and no error
00070 *          message related to LWORK is issued by XERBLA.
00071 *
00072 *  INFO    (output) INTEGER
00073 *          = 0:  successful exit
00074 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00075 *
00076 *  Further Details
00077 *  ===============
00078 *
00079 *  The matrix Q is represented as a product of elementary reflectors
00080 *
00081 *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
00082 *
00083 *  Each H(i) has the form
00084 *
00085 *     H(i) = I - tau * v * v'
00086 *
00087 *  where tau is a real scalar, and v is a real vector with
00088 *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
00089 *  and tau in TAU(i).
00090 *
00091 *  =====================================================================
00092 *
00093 *     .. Local Scalars ..
00094       LOGICAL            LQUERY
00095       INTEGER            I, IB, IINFO, IWS, J, K, LWKOPT, NB,
00096      \$                   NBMIN, NX, LBWORK, NT, LLWORK
00097 *     ..
00098 *     .. External Subroutines ..
00099       EXTERNAL           DGEQR2, DLARFB, DLARFT, XERBLA
00100 *     ..
00101 *     .. Intrinsic Functions ..
00102       INTRINSIC          MAX, MIN
00103 *     ..
00104 *     .. External Functions ..
00105       INTEGER            ILAENV
00106       REAL               SCEIL
00107       EXTERNAL           ILAENV, SCEIL
00108 *     ..
00109 *     .. Executable Statements ..
00110
00111       INFO = 0
00112       NBMIN = 2
00113       NX = 0
00114       IWS = N
00115       K = MIN( M, N )
00116       NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00117
00118       IF( NB.GT.1 .AND. NB.LT.K ) THEN
00119 *
00120 *        Determine when to cross over from blocked to unblocked code.
00121 *
00122          NX = MAX( 0, ILAENV( 3, 'DGEQRF', ' ', M, N, -1, -1 ) )
00123       END IF
00124 *
00125 *     Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
00126 *
00127 *            NB=3     2NB=6       K=10
00128 *            |        |           |
00129 *      1--2--3--4--5--6--7--8--9--10
00130 *                  |     \________/
00131 *               K-NX=5      NT=4
00132 *
00133 *     So here 4 x 4 is the last T stored in the workspace
00134 *
00135       NT = K-SCEIL(REAL(K-NX)/REAL(NB))*NB
00136
00137 *
00138 *     optimal workspace = space for dlarfb + space for normal T's + space for the last T
00139 *
00140       LLWORK = MAX (MAX((N-M)*K, (N-M)*NB), MAX(K*NB, NB*NB))
00141       LLWORK = SCEIL(REAL(LLWORK)/REAL(NB))
00142
00143       IF ( NT.GT.NB ) THEN
00144
00145           LBWORK = K-NT
00146 *
00147 *         Optimal workspace for dlarfb = MAX(1,N)*NT
00148 *
00149           LWKOPT = (LBWORK+LLWORK)*NB
00150           WORK( 1 ) = (LWKOPT+NT*NT)
00151
00152       ELSE
00153
00154           LBWORK = SCEIL(REAL(K)/REAL(NB))*NB
00155           LWKOPT = (LBWORK+LLWORK-NB)*NB
00156           WORK( 1 ) = LWKOPT
00157
00158       END IF
00159
00160 *
00161 *     Test the input arguments
00162 *
00163       LQUERY = ( LWORK.EQ.-1 )
00164       IF( M.LT.0 ) THEN
00165          INFO = -1
00166       ELSE IF( N.LT.0 ) THEN
00167          INFO = -2
00168       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00169          INFO = -4
00170       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
00171          INFO = -7
00172       END IF
00173       IF( INFO.NE.0 ) THEN
00174          CALL XERBLA( 'DGEQRF', -INFO )
00175          RETURN
00176       ELSE IF( LQUERY ) THEN
00177          RETURN
00178       END IF
00179 *
00180 *     Quick return if possible
00181 *
00182       IF( K.EQ.0 ) THEN
00183          WORK( 1 ) = 1
00184          RETURN
00185       END IF
00186 *
00187       IF( NB.GT.1 .AND. NB.LT.K ) THEN
00188
00189          IF( NX.LT.K ) THEN
00190 *
00191 *           Determine if workspace is large enough for blocked code.
00192 *
00193             IF ( NT.LE.NB ) THEN
00194                 IWS = (LBWORK+LLWORK-NB)*NB
00195             ELSE
00196                 IWS = (LBWORK+LLWORK)*NB+NT*NT
00197             END IF
00198
00199             IF( LWORK.LT.IWS ) THEN
00200 *
00201 *              Not enough workspace to use optimal NB:  reduce NB and
00202 *              determine the minimum value of NB.
00203 *
00204                IF ( NT.LE.NB ) THEN
00205                     NB = LWORK / (LLWORK+(LBWORK-NB))
00206                ELSE
00207                     NB = (LWORK-NT*NT)/(LBWORK+LLWORK)
00208                END IF
00209
00210                NBMIN = MAX( 2, ILAENV( 2, 'DGEQRF', ' ', M, N, -1,
00211      \$                 -1 ) )
00212             END IF
00213          END IF
00214       END IF
00215 *
00216       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
00217 *
00218 *        Use blocked code initially
00219 *
00220          DO 10 I = 1, K - NX, NB
00221             IB = MIN( K-I+1, NB )
00222 *
00223 *           Update the current column using old T's
00224 *
00225             DO 20 J = 1, I - NB, NB
00226 *
00227 *              Apply H' to A(J:M,I:I+IB-1) from the left
00228 *
00229                CALL DLARFB( 'Left', 'Transpose', 'Forward',
00230      \$                      'Columnwise', M-J+1, IB, NB,
00231      \$                      A( J, J ), LDA, WORK(J), LBWORK,
00232      \$                      A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1),
00233      \$                      IB)
00234
00235 20          CONTINUE
00236 *
00237 *           Compute the QR factorization of the current block
00238 *           A(I:M,I:I+IB-1)
00239 *
00240             CALL DGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ),
00241      \$                        WORK(LBWORK*NB+NT*NT+1), IINFO )
00242
00243             IF( I+IB.LE.N ) THEN
00244 *
00245 *              Form the triangular factor of the block reflector
00246 *              H = H(i) H(i+1) . . . H(i+ib-1)
00247 *
00248                CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB,
00249      \$                      A( I, I ), LDA, TAU( I ),
00250      \$                      WORK(I), LBWORK )
00251 *
00252             END IF
00253    10    CONTINUE
00254       ELSE
00255          I = 1
00256       END IF
00257 *
00258 *     Use unblocked code to factor the last or only block.
00259 *
00260       IF( I.LE.K ) THEN
00261
00262          IF ( I .NE. 1 )   THEN
00263
00264              DO 30 J = 1, I - NB, NB
00265 *
00266 *                Apply H' to A(J:M,I:K) from the left
00267 *
00268                  CALL DLARFB( 'Left', 'Transpose', 'Forward',
00269      \$                       'Columnwise', M-J+1, K-I+1, NB,
00270      \$                       A( J, J ), LDA, WORK(J), LBWORK,
00271      \$                       A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1),
00272      \$                       K-I+1)
00273 30           CONTINUE
00274
00275              CALL DGEQR2( M-I+1, K-I+1, A( I, I ), LDA, TAU( I ),
00276      \$                   WORK(LBWORK*NB+NT*NT+1),IINFO )
00277
00278          ELSE
00279 *
00280 *        Use unblocked code to factor the last or only block.
00281 *
00282          CALL DGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ),
00283      \$               WORK,IINFO )
00284
00285          END IF
00286       END IF
00287
00288
00289 *
00290 *     Apply update to the column M+1:N when N > M
00291 *
00292       IF ( M.LT.N .AND. I.NE.1) THEN
00293 *
00294 *         Form the last triangular factor of the block reflector
00295 *         H = H(i) H(i+1) . . . H(i+ib-1)
00296 *
00297           IF ( NT .LE. NB ) THEN
00298                CALL DLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1,
00299      \$                     A( I, I ), LDA, TAU( I ), WORK(I), LBWORK )
00300           ELSE
00301                CALL DLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1,
00302      \$                     A( I, I ), LDA, TAU( I ),
00303      \$                     WORK(LBWORK*NB+1), NT )
00304           END IF
00305
00306 *
00307 *         Apply H' to A(1:M,M+1:N) from the left
00308 *
00309           DO 40 J = 1, K-NX, NB
00310
00311                IB = MIN( K-J+1, NB )
00312
00313                CALL DLARFB( 'Left', 'Transpose', 'Forward',
00314      \$                     'Columnwise', M-J+1, N-M, IB,
00315      \$                     A( J, J ), LDA, WORK(J), LBWORK,
00316      \$                     A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
00317      \$                     N-M)
00318
00319 40       CONTINUE
00320
00321          IF ( NT.LE.NB ) THEN
00322              CALL DLARFB( 'Left', 'Transpose', 'Forward',
00323      \$                   'Columnwise', M-J+1, N-M, K-J+1,
00324      \$                   A( J, J ), LDA, WORK(J), LBWORK,
00325      \$                   A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
00326      \$                   N-M)
00327          ELSE
00328              CALL DLARFB( 'Left', 'Transpose', 'Forward',
00329      \$                   'Columnwise', M-J+1, N-M, K-J+1,
00330      \$                   A( J, J ), LDA,
00331      \$                   WORK(LBWORK*NB+1),
00332      \$                   NT, A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
00333      \$                   N-M)
00334          END IF
00335
00336       END IF
00337
00338       WORK( 1 ) = IWS
00339       RETURN
00340 *
00341 *     End of DGEQRF
00342 *
00343       END
00344
```