SRC/VARIANTS/qr/LL/cgeqrf.f

   1 C> \brief \b CGEQRF VARIANT: left-looking Level 3 BLAS version of the algorithm.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       SUBROUTINE CGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       INTEGER            INFO, LDA, LWORK, M, N
  15 *       ..
  16 *       .. Array Arguments ..
  17 *       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
  18 *       ..
  19 *
  20 *  Purpose
  21 *  =======
  22 *
  23 C>\details \b Purpose:
  24 C>\verbatim
  25 C>
  26 C> CGEQRF computes a QR factorization of a real M-by-N matrix A:
  27 C> A = Q * R.
  28 C>
  29 C> This is the left-looking Level 3 BLAS version of the algorithm.
  30 C>
  31 C>\endverbatim
  32 *
  33 *  Arguments:
  34 *  ==========
  35 *
  36 C> \param[in] M
  37 C> \verbatim
  38 C>          M is INTEGER
  39 C>          The number of rows of the matrix A.  M >= 0.
  40 C> \endverbatim
  41 C>
  42 C> \param[in] N
  43 C> \verbatim
  44 C>          N is INTEGER
  45 C>          The number of columns of the matrix A.  N >= 0.
  46 C> \endverbatim
  47 C>
  48 C> \param[in,out] A
  49 C> \verbatim
  50 C>          A is COMPLEX array, dimension (LDA,N)
  51 C>          On entry, the M-by-N matrix A.
  52 C>          On exit, the elements on and above the diagonal of the array
  53 C>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
  54 C>          upper triangular if m >= n); the elements below the diagonal,
  55 C>          with the array TAU, represent the orthogonal matrix Q as a
  56 C>          product of min(m,n) elementary reflectors (see Further
  57 C>          Details).
  58 C> \endverbatim
  59 C>
  60 C> \param[in] LDA
  61 C> \verbatim
  62 C>          LDA is INTEGER
  63 C>          The leading dimension of the array A.  LDA >= max(1,M).
  64 C> \endverbatim
  65 C>
  66 C> \param[out] TAU
  67 C> \verbatim
  68 C>          TAU is COMPLEX array, dimension (min(M,N))
  69 C>          The scalar factors of the elementary reflectors (see Further
  70 C>          Details).
  71 C> \endverbatim
  72 C>
  73 C> \param[out] WORK
  74 C> \verbatim
  75 C>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
  76 C>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  77 C> \endverbatim
  78 C>
  79 C> \param[in] LWORK
  80 C> \verbatim
  81 C>          LWORK is INTEGER
  82 C> \endverbatim
  83 C> \verbatim
  84 C>          The dimension of the array WORK. The dimension can be divided into three parts.
  85 C> \endverbatim
  86 C> \verbatim
  87 C>          1) The part for the triangular factor T. If the very last T is not bigger
  88 C>             than any of the rest, then this part is NB x ceiling(K/NB), otherwise,
  89 C>             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T
  90 C> \endverbatim
  91 C> \verbatim
  92 C>          2) The part for the very last T when T is bigger than any of the rest T.
  93 C>             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB,
  94 C>             where K = min(M,N), NX is calculated by
  95 C>                   NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) )
  96 C> \endverbatim
  97 C> \verbatim
  98 C>          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB)
  99 C> \endverbatim
 100 C> \verbatim
 101 C>          So LWORK = part1 + part2 + part3
 102 C> \endverbatim
 103 C> \verbatim
 104 C>          If LWORK = -1, then a workspace query is assumed; the routine
 105 C>          only calculates the optimal size of the WORK array, returns
 106 C>          this value as the first entry of the WORK array, and no error
 107 C>          message related to LWORK is issued by XERBLA.
 108 C> \endverbatim
 109 C>
 110 C> \param[out] INFO
 111 C> \verbatim
 112 C>          INFO is INTEGER
 113 C>          = 0:  successful exit
 114 C>          < 0:  if INFO = -i, the i-th argument had an illegal value
 115 C> \endverbatim
 116 C>
 117 *
 118 *  Authors:
 119 *  ========
 120 *
 121 C> \author Univ. of Tennessee
 122 C> \author Univ. of California Berkeley
 123 C> \author Univ. of Colorado Denver
 124 C> \author NAG Ltd.
 125 *
 126 C> \date November 2011
 127 *
 128 C> \ingroup variantsGEcomputational
 129 *
 130 *  Further Details
 131 *  ===============
 132 C>\details \b Further \b Details
 133 C> \verbatim
 134 C>
 135 C>  The matrix Q is represented as a product of elementary reflectors
 136 C>
 137 C>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
 138 C>
 139 C>  Each H(i) has the form
 140 C>
 141 C>     H(i) = I - tau * v * v'
 142 C>
 143 C>  where tau is a real scalar, and v is a real vector with
 144 C>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
 145 C>  and tau in TAU(i).
 146 C>
 147 C> \endverbatim
 148 C>
 149 *  =====================================================================
 150       SUBROUTINE CGEQRF ( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 151 *
 152 *  -- LAPACK computational routine (version 3.1) --
 153 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 154 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 155 *     November 2011
 156 *
 157 *     .. Scalar Arguments ..
 158       INTEGER            INFO, LDA, LWORK, M, N
 159 *     ..
 160 *     .. Array Arguments ..
 161       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
 162 *     ..
 163 *
 164 *  =====================================================================
 165 *
 166 *     .. Local Scalars ..
 167       LOGICAL            LQUERY
 168       INTEGER            I, IB, IINFO, IWS, J, K, LWKOPT, NB,
 169      $                   NBMIN, NX, LBWORK, NT, LLWORK
 170 *     ..
 171 *     .. External Subroutines ..
 172       EXTERNAL           CGEQR2, CLARFB, CLARFT, XERBLA
 173 *     ..
 174 *     .. Intrinsic Functions ..
 175       INTRINSIC          MAX, MIN
 176 *     ..
 177 *     .. External Functions ..
 178       INTEGER            ILAENV
 179       REAL               SCEIL
 180       EXTERNAL           ILAENV, SCEIL
 181 *     ..
 182 *     .. Executable Statements ..
 183
 184       INFO = 0
 185       NBMIN = 2
 186       NX = 0
 187       IWS = N
 188       K = MIN( M, N )
 189       NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
 190
 191       IF( NB.GT.1 .AND. NB.LT.K ) THEN
 192 *
 193 *        Determine when to cross over from blocked to unblocked code.
 194 *
 195          NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) )
 196       END IF
 197 *
 198 *     Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.:
 199 *
 200 *            NB=3     2NB=6       K=10
 201 *            |        |           |
 202 *      1--2--3--4--5--6--7--8--9--10
 203 *                  |     \________/
 204 *               K-NX=5      NT=4
 205 *
 206 *     So here 4 x 4 is the last T stored in the workspace
 207 *
 208       NT = K-SCEIL(REAL(K-NX)/REAL(NB))*NB
 209
 210 *
 211 *     optimal workspace = space for dlarfb + space for normal T's + space for the last T
 212 *
 213       LLWORK = MAX (MAX((N-M)*K, (N-M)*NB), MAX(K*NB, NB*NB))
 214       LLWORK = SCEIL(REAL(LLWORK)/REAL(NB))
 215
 216       IF ( NT.GT.NB ) THEN
 217
 218           LBWORK = K-NT
 219 *
 220 *         Optimal workspace for dlarfb = MAX(1,N)*NT
 221 *
 222           LWKOPT = (LBWORK+LLWORK)*NB
 223           WORK( 1 ) = (LWKOPT+NT*NT)
 224
 225       ELSE
 226
 227           LBWORK = SCEIL(REAL(K)/REAL(NB))*NB
 228           LWKOPT = (LBWORK+LLWORK-NB)*NB
 229           WORK( 1 ) = LWKOPT
 230
 231       END IF
 232
 233 *
 234 *     Test the input arguments
 235 *
 236       LQUERY = ( LWORK.EQ.-1 )
 237       IF( M.LT.0 ) THEN
 238          INFO = -1
 239       ELSE IF( N.LT.0 ) THEN
 240          INFO = -2
 241       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 242          INFO = -4
 243       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 244          INFO = -7
 245       END IF
 246       IF( INFO.NE.0 ) THEN
 247          CALL XERBLA( 'CGEQRF', -INFO )
 248          RETURN
 249       ELSE IF( LQUERY ) THEN
 250          RETURN
 251       END IF
 252 *
 253 *     Quick return if possible
 254 *
 255       IF( K.EQ.0 ) THEN
 256          WORK( 1 ) = 1
 257          RETURN
 258       END IF
 259 *
 260       IF( NB.GT.1 .AND. NB.LT.K ) THEN
 261
 262          IF( NX.LT.K ) THEN
 263 *
 264 *           Determine if workspace is large enough for blocked code.
 265 *
 266             IF ( NT.LE.NB ) THEN
 267                 IWS = (LBWORK+LLWORK-NB)*NB
 268             ELSE
 269                 IWS = (LBWORK+LLWORK)*NB+NT*NT
 270             END IF
 271
 272             IF( LWORK.LT.IWS ) THEN
 273 *
 274 *              Not enough workspace to use optimal NB:  reduce NB and
 275 *              determine the minimum value of NB.
 276 *
 277                IF ( NT.LE.NB ) THEN
 278                     NB = LWORK / (LLWORK+(LBWORK-NB))
 279                ELSE
 280                     NB = (LWORK-NT*NT)/(LBWORK+LLWORK)
 281                END IF
 282
 283                NBMIN = MAX( 2, ILAENV( 2, 'CGEQRF', ' ', M, N, -1,
 284      $                 -1 ) )
 285             END IF
 286          END IF
 287       END IF
 288 *
 289       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
 290 *
 291 *        Use blocked code initially
 292 *
 293          DO 10 I = 1, K - NX, NB
 294             IB = MIN( K-I+1, NB )
 295 *
 296 *           Update the current column using old T's
 297 *
 298             DO 20 J = 1, I - NB, NB
 299 *
 300 *              Apply H' to A(J:M,I:I+IB-1) from the left
 301 *
 302                CALL CLARFB( 'Left', 'Transpose', 'Forward',
 303      $                      'Columnwise', M-J+1, IB, NB,
 304      $                      A( J, J ), LDA, WORK(J), LBWORK,
 305      $                      A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1),
 306      $                      IB)
 307
 308 20          CONTINUE
 309 *
 310 *           Compute the QR factorization of the current block
 311 *           A(I:M,I:I+IB-1)
 312 *
 313             CALL CGEQR2( M-I+1, IB, A( I, I ), LDA, TAU( I ),
 314      $                        WORK(LBWORK*NB+NT*NT+1), IINFO )
 315
 316             IF( I+IB.LE.N ) THEN
 317 *
 318 *              Form the triangular factor of the block reflector
 319 *              H = H(i) H(i+1) . . . H(i+ib-1)
 320 *
 321                CALL CLARFT( 'Forward', 'Columnwise', M-I+1, IB,
 322      $                      A( I, I ), LDA, TAU( I ),
 323      $                      WORK(I), LBWORK )
 324 *
 325             END IF
 326    10    CONTINUE
 327       ELSE
 328          I = 1
 329       END IF
 330 *
 331 *     Use unblocked code to factor the last or only block.
 332 *
 333       IF( I.LE.K ) THEN
 334
 335          IF ( I .NE. 1 )   THEN
 336
 337              DO 30 J = 1, I - NB, NB
 338 *
 339 *                Apply H' to A(J:M,I:K) from the left
 340 *
 341                  CALL CLARFB( 'Left', 'Transpose', 'Forward',
 342      $                       'Columnwise', M-J+1, K-I+1, NB,
 343      $                       A( J, J ), LDA, WORK(J), LBWORK,
 344      $                       A( J, I ), LDA, WORK(LBWORK*NB+NT*NT+1),
 345      $                       K-I+1)
 346 30           CONTINUE
 347
 348              CALL CGEQR2( M-I+1, K-I+1, A( I, I ), LDA, TAU( I ),
 349      $                   WORK(LBWORK*NB+NT*NT+1),IINFO )
 350
 351          ELSE
 352 *
 353 *        Use unblocked code to factor the last or only block.
 354 *
 355          CALL CGEQR2( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ),
 356      $               WORK,IINFO )
 357
 358          END IF
 359       END IF
 360
 361
 362 *
 363 *     Apply update to the column M+1:N when N > M
 364 *
 365       IF ( M.LT.N .AND. I.NE.1) THEN
 366 *
 367 *         Form the last triangular factor of the block reflector
 368 *         H = H(i) H(i+1) . . . H(i+ib-1)
 369 *
 370           IF ( NT .LE. NB ) THEN
 371                CALL CLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1,
 372      $                     A( I, I ), LDA, TAU( I ), WORK(I), LBWORK )
 373           ELSE
 374                CALL CLARFT( 'Forward', 'Columnwise', M-I+1, K-I+1,
 375      $                     A( I, I ), LDA, TAU( I ),
 376      $                     WORK(LBWORK*NB+1), NT )
 377           END IF
 378
 379 *
 380 *         Apply H' to A(1:M,M+1:N) from the left
 381 *
 382           DO 40 J = 1, K-NX, NB
 383
 384                IB = MIN( K-J+1, NB )
 385
 386                CALL CLARFB( 'Left', 'Transpose', 'Forward',
 387      $                     'Columnwise', M-J+1, N-M, IB,
 388      $                     A( J, J ), LDA, WORK(J), LBWORK,
 389      $                     A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
 390      $                     N-M)
 391
 392 40       CONTINUE
 393
 394          IF ( NT.LE.NB ) THEN
 395              CALL CLARFB( 'Left', 'Transpose', 'Forward',
 396      $                   'Columnwise', M-J+1, N-M, K-J+1,
 397      $                   A( J, J ), LDA, WORK(J), LBWORK,
 398      $                   A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
 399      $                   N-M)
 400          ELSE
 401              CALL CLARFB( 'Left', 'Transpose', 'Forward',
 402      $                   'Columnwise', M-J+1, N-M, K-J+1,
 403      $                   A( J, J ), LDA,
 404      $                   WORK(LBWORK*NB+1),
 405      $                   NT, A( J, M+1 ), LDA, WORK(LBWORK*NB+NT*NT+1),
 406      $                   N-M)
 407          END IF
 408
 409       END IF
 410
 411       WORK( 1 ) = IWS
 412       RETURN
 413 *
 414 *     End of CGEQRF
 415 *
 416       END