SRC/cgelqt.f

   1 *
   2 *  Definition:
   3 *  ===========
   4 *
   5 *       SUBROUTINE CGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
   6 *
   7 *       .. Scalar Arguments ..
   8 *       INTEGER INFO, LDA, LDT, M, N, MB
   9 *       ..
  10 *       .. Array Arguments ..
  11 *       COMPLEX A( LDA, * ), T( LDT, * ), WORK( * )
  12 *       ..
  13 *
  14 *
  15 *> \par Purpose:
  16 *  =============
  17 *>
  18 *> \verbatim
  19 *>
  20 *> CGELQT computes a blocked LQ factorization of a complex M-by-N matrix A
  21 *> using the compact WY representation of Q.
  22 *> \endverbatim
  23 *
  24 *  Arguments:
  25 *  ==========
  26 *
  27 *> \param[in] M
  28 *> \verbatim
  29 *>          M is INTEGER
  30 *>          The number of rows of the matrix A.  M >= 0.
  31 *> \endverbatim
  32 *>
  33 *> \param[in] N
  34 *> \verbatim
  35 *>          N is INTEGER
  36 *>          The number of columns of the matrix A.  N >= 0.
  37 *> \endverbatim
  38 *>
  39 *> \param[in] MB
  40 *> \verbatim
  41 *>          MB is INTEGER
  42 *>          The block size to be used in the blocked QR.  MIN(M,N) >= MB >= 1.
  43 *> \endverbatim
  44 *>
  45 *> \param[in,out] A
  46 *> \verbatim
  47 *>          A is COMPLEX array, dimension (LDA,N)
  48 *>          On entry, the M-by-N matrix A.
  49 *>          On exit, the elements on and below the diagonal of the array
  50 *>          contain the M-by-MIN(M,N) lower trapezoidal matrix L (L is
  51 *>          lower triangular if M <= N); the elements above the diagonal
  52 *>          are the rows of V.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] LDA
  56 *> \verbatim
  57 *>          LDA is INTEGER
  58 *>          The leading dimension of the array A.  LDA >= max(1,M).
  59 *> \endverbatim
  60 *>
  61 *> \param[out] T
  62 *> \verbatim
  63 *>          T is COMPLEX array, dimension (LDT,MIN(M,N))
  64 *>          The upper triangular block reflectors stored in compact form
  65 *>          as a sequence of upper triangular blocks.  See below
  66 *>          for further details.
  67 *> \endverbatim
  68 *>
  69 *> \param[in] LDT
  70 *> \verbatim
  71 *>          LDT is INTEGER
  72 *>          The leading dimension of the array T.  LDT >= MB.
  73 *> \endverbatim
  74 *>
  75 *> \param[out] WORK
  76 *> \verbatim
  77 *>          WORK is COMPLEX array, dimension (MB*N)
  78 *> \endverbatim
  79 *>
  80 *> \param[out] INFO
  81 *> \verbatim
  82 *>          INFO is INTEGER
  83 *>          = 0:  successful exit
  84 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  85 *> \endverbatim
  86 *
  87 *  Authors:
  88 *  ========
  89 *
  90 *> \author Univ. of Tennessee
  91 *> \author Univ. of California Berkeley
  92 *> \author Univ. of Colorado Denver
  93 *> \author NAG Ltd.
  94 *
  95 *> \date November 2013
  96 *
  97 *> \ingroup doubleGEcomputational
  98 *
  99 *> \par Further Details:
 100 *  =====================
 101 *>
 102 *> \verbatim
 103 *>
 104 *>  The matrix V stores the elementary reflectors H(i) in the i-th column
 105 *>  below the diagonal. For example, if M=5 and N=3, the matrix V is
 106 *>
 107 *>               V = (  1  v1 v1 v1 v1 )
 108 *>                   (     1  v2 v2 v2 )
 109 *>                   (         1 v3 v3 )
 110 *>
 111 *>
 112 *>  where the vi's represent the vectors which define H(i), which are returned
 113 *>  in the matrix A.  The 1's along the diagonal of V are not stored in A.
 114 *>  Let K=MIN(M,N).  The number of blocks is B = ceiling(K/NB), where each
 115 *>  block is of order NB except for the last block, which is of order
 116 *>  IB = K - (B-1)*NB.  For each of the B blocks, a upper triangular block
 117 *>  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB
 118 *>  for the last block) T's are stored in the NB-by-N matrix T as
 119 *>
 120 *>               T = (T1 T2 ... TB).
 121 *> \endverbatim
 122 *>
 123 *  =====================================================================
 124       SUBROUTINE CGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
 125 *
 126 *  -- LAPACK computational routine (version 3.5.0) --
 127 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 128 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 129 *     November 2013
 130 *
 131 *     .. Scalar Arguments ..
 132       INTEGER INFO, LDA, LDT, M, N, MB
 133 *     ..
 134 *     .. Array Arguments ..
 135       COMPLEX A( LDA, * ), T( LDT, * ), WORK( * )
 136 *     ..
 137 *
 138 * =====================================================================
 139 *
 140 *     ..
 141 *     .. Local Scalars ..
 142       INTEGER    I, IB, IINFO, K
 143 *     ..
 144 *     .. External Subroutines ..
 145       EXTERNAL   CGELQT3, CLARFB, XERBLA
 146 *     ..
 147 *     .. Executable Statements ..
 148 *
 149 *     Test the input arguments
 150 *
 151       INFO = 0
 152       IF( M.LT.0 ) THEN
 153          INFO = -1
 154       ELSE IF( N.LT.0 ) THEN
 155          INFO = -2
 156       ELSE IF( MB.LT.1 .OR. (MB.GT.MIN(M,N) .AND. MIN(M,N).GT.0 ))THEN
 157          INFO = -3
 158       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 159          INFO = -5
 160       ELSE IF( LDT.LT.MB ) THEN
 161          INFO = -7
 162       END IF
 163       IF( INFO.NE.0 ) THEN
 164          CALL XERBLA( 'CGELQT', -INFO )
 165          RETURN
 166       END IF
 167 *
 168 *     Quick return if possible
 169 *
 170       K = MIN( M, N )
 171       IF( K.EQ.0 ) RETURN
 172 *
 173 *     Blocked loop of length K
 174 *
 175       DO I = 1, K,  MB
 176          IB = MIN( K-I+1, MB )
 177 *
 178 *     Compute the LQ factorization of the current block A(I:M,I:I+IB-1)
 179 *
 180          CALL CGELQT3( IB, N-I+1, A(I,I), LDA, T(1,I), LDT, IINFO )
 181          IF( I+IB.LE.M ) THEN
 182 *
 183 *     Update by applying H**T to A(I:M,I+IB:N) from the right
 184 *
 185          CALL CLARFB( 'R', 'N', 'F', 'R', M-I-IB+1, N-I+1, IB,
 186      $                   A( I, I ), LDA, T( 1, I ), LDT,
 187      $                   A( I+IB, I ), LDA, WORK , M-I-IB+1 )
 188          END IF
 189       END DO
 190       RETURN
 191 *
 192 *     End of CGELQT
 193 *
 194       END