SRC/sgelqt.f

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