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