SRC/cgeqrfp.f

   1 *> \brief \b CGEQRFP
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CGEQRFP + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeqrfp.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeqrfp.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeqrfp.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, LDA, LWORK, M, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> CGEQRFP computes a QR factorization of a complex M-by-N matrix A:
  37 *> A = Q * R. The diagonal entries of R are real and nonnegative.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] M
  44 *> \verbatim
  45 *>          M is INTEGER
  46 *>          The number of rows of the matrix A.  M >= 0.
  47 *> \endverbatim
  48 *>
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>          The number of columns of the matrix A.  N >= 0.
  53 *> \endverbatim
  54 *>
  55 *> \param[in,out] A
  56 *> \verbatim
  57 *>          A is COMPLEX array, dimension (LDA,N)
  58 *>          On entry, the M-by-N matrix A.
  59 *>          On exit, the elements on and above the diagonal of the array
  60 *>          contain the min(M,N)-by-N upper trapezoidal matrix R (R is
  61 *>          upper triangular if m >= n). The diagonal entries of R
  62 *>          are real and nonnegative; the elements below the diagonal,
  63 *>          with the array TAU, represent the unitary matrix Q as a
  64 *>          product of min(m,n) elementary reflectors (see Further
  65 *>          Details).
  66 *> \endverbatim
  67 *>
  68 *> \param[in] LDA
  69 *> \verbatim
  70 *>          LDA is INTEGER
  71 *>          The leading dimension of the array A.  LDA >= max(1,M).
  72 *> \endverbatim
  73 *>
  74 *> \param[out] TAU
  75 *> \verbatim
  76 *>          TAU is COMPLEX array, dimension (min(M,N))
  77 *>          The scalar factors of the elementary reflectors (see Further
  78 *>          Details).
  79 *> \endverbatim
  80 *>
  81 *> \param[out] WORK
  82 *> \verbatim
  83 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
  84 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] LWORK
  88 *> \verbatim
  89 *>          LWORK is INTEGER
  90 *>          The dimension of the array WORK.  LWORK >= max(1,N).
  91 *>          For optimum performance LWORK >= N*NB, where NB is
  92 *>          the optimal blocksize.
  93 *>
  94 *>          If LWORK = -1, then a workspace query is assumed; the routine
  95 *>          only calculates the optimal size of the WORK array, returns
  96 *>          this value as the first entry of the WORK array, and no error
  97 *>          message related to LWORK is issued by XERBLA.
  98 *> \endverbatim
  99 *>
 100 *> \param[out] INFO
 101 *> \verbatim
 102 *>          INFO is INTEGER
 103 *>          = 0:  successful exit
 104 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 105 *> \endverbatim
 106 *
 107 *  Authors:
 108 *  ========
 109 *
 110 *> \author Univ. of Tennessee
 111 *> \author Univ. of California Berkeley
 112 *> \author Univ. of Colorado Denver
 113 *> \author NAG Ltd.
 114 *
 115 *> \date November 2015
 116 *
 117 *> \ingroup complexGEcomputational
 118 *
 119 *> \par Further Details:
 120 *  =====================
 121 *>
 122 *> \verbatim
 123 *>
 124 *>  The matrix Q is represented as a product of elementary reflectors
 125 *>
 126 *>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
 127 *>
 128 *>  Each H(i) has the form
 129 *>
 130 *>     H(i) = I - tau * v * v**H
 131 *>
 132 *>  where tau is a complex scalar, and v is a complex vector with
 133 *>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
 134 *>  and tau in TAU(i).
 135 *>
 136 *> See Lapack Working Note 203 for details
 137 *> \endverbatim
 138 *>
 139 *  =====================================================================
 140       SUBROUTINE CGEQRFP( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 141 *
 142 *  -- LAPACK computational routine (version 3.6.0) --
 143 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 144 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 145 *     November 2015
 146 *
 147 *     .. Scalar Arguments ..
 148       INTEGER            INFO, LDA, LWORK, M, N
 149 *     ..
 150 *     .. Array Arguments ..
 151       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
 152 *     ..
 153 *
 154 *  =====================================================================
 155 *
 156 *     .. Local Scalars ..
 157       LOGICAL            LQUERY
 158       INTEGER            I, IB, IINFO, IWS, K, LDWORK, LWKOPT, NB,
 159      $                   NBMIN, NX
 160 *     ..
 161 *     .. External Subroutines ..
 162       EXTERNAL           CGEQR2P, CLARFB, CLARFT, XERBLA
 163 *     ..
 164 *     .. Intrinsic Functions ..
 165       INTRINSIC          MAX, MIN
 166 *     ..
 167 *     .. External Functions ..
 168       INTEGER            ILAENV
 169       EXTERNAL           ILAENV
 170 *     ..
 171 *     .. Executable Statements ..
 172 *
 173 *     Test the input arguments
 174 *
 175       INFO = 0
 176       NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
 177       LWKOPT = N*NB
 178       WORK( 1 ) = LWKOPT
 179       LQUERY = ( LWORK.EQ.-1 )
 180       IF( M.LT.0 ) THEN
 181          INFO = -1
 182       ELSE IF( N.LT.0 ) THEN
 183          INFO = -2
 184       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 185          INFO = -4
 186       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 187          INFO = -7
 188       END IF
 189       IF( INFO.NE.0 ) THEN
 190          CALL XERBLA( 'CGEQRFP', -INFO )
 191          RETURN
 192       ELSE IF( LQUERY ) THEN
 193          RETURN
 194       END IF
 195 *
 196 *     Quick return if possible
 197 *
 198       K = MIN( M, N )
 199       IF( K.EQ.0 ) THEN
 200          WORK( 1 ) = 1
 201          RETURN
 202       END IF
 203 *
 204       NBMIN = 2
 205       NX = 0
 206       IWS = N
 207       IF( NB.GT.1 .AND. NB.LT.K ) THEN
 208 *
 209 *        Determine when to cross over from blocked to unblocked code.
 210 *
 211          NX = MAX( 0, ILAENV( 3, 'CGEQRF', ' ', M, N, -1, -1 ) )
 212          IF( NX.LT.K ) THEN
 213 *
 214 *           Determine if workspace is large enough for blocked code.
 215 *
 216             LDWORK = N
 217             IWS = LDWORK*NB
 218             IF( LWORK.LT.IWS ) THEN
 219 *
 220 *              Not enough workspace to use optimal NB:  reduce NB and
 221 *              determine the minimum value of NB.
 222 *
 223                NB = LWORK / LDWORK
 224                NBMIN = MAX( 2, ILAENV( 2, 'CGEQRF', ' ', M, N, -1,
 225      $                 -1 ) )
 226             END IF
 227          END IF
 228       END IF
 229 *
 230       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
 231 *
 232 *        Use blocked code initially
 233 *
 234          DO 10 I = 1, K - NX, NB
 235             IB = MIN( K-I+1, NB )
 236 *
 237 *           Compute the QR factorization of the current block
 238 *           A(i:m,i:i+ib-1)
 239 *
 240             CALL CGEQR2P( M-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
 241      $                   IINFO )
 242             IF( I+IB.LE.N ) THEN
 243 *
 244 *              Form the triangular factor of the block reflector
 245 *              H = H(i) H(i+1) . . . H(i+ib-1)
 246 *
 247                CALL CLARFT( 'Forward', 'Columnwise', M-I+1, IB,
 248      $                      A( I, I ), LDA, TAU( I ), WORK, LDWORK )
 249 *
 250 *              Apply H**H to A(i:m,i+ib:n) from the left
 251 *
 252                CALL CLARFB( 'Left', 'Conjugate transpose', 'Forward',
 253      $                      'Columnwise', M-I+1, N-I-IB+1, IB,
 254      $                      A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
 255      $                      LDA, WORK( IB+1 ), LDWORK )
 256             END IF
 257    10    CONTINUE
 258       ELSE
 259          I = 1
 260       END IF
 261 *
 262 *     Use unblocked code to factor the last or only block.
 263 *
 264       IF( I.LE.K )
 265      $   CALL CGEQR2P( M-I+1, N-I+1, A( I, I ), LDA, TAU( I ), WORK,
 266      $                IINFO )
 267 *
 268       WORK( 1 ) = IWS
 269       RETURN
 270 *
 271 *     End of CGEQRFP
 272 *
 273       END