SRC/dgelqt.f

   1 *> \brief \b DGELQT
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGEQRT + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelqt.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelqt.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelqt.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER INFO, LDA, LDT, M, N, MB
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DGELQT computes a blocked LQ factorization of a real M-by-N matrix A
  37 *> using the compact WY representation of Q.
  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] MB
  56 *> \verbatim
  57 *>          MB is INTEGER
  58 *>          The block size to be used in the blocked QR.  MIN(M,N) >= MB >= 1.
  59 *> \endverbatim
  60 *>
  61 *> \param[in,out] A
  62 *> \verbatim
  63 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  64 *>          On entry, the M-by-N matrix A.
  65 *>          On exit, the elements on and below the diagonal of the array
  66 *>          contain the M-by-MIN(M,N) lower trapezoidal matrix L (L is
  67 *>          lower triangular if M <= N); the elements above the diagonal
  68 *>          are the rows of V.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] LDA
  72 *> \verbatim
  73 *>          LDA is INTEGER
  74 *>          The leading dimension of the array A.  LDA >= max(1,M).
  75 *> \endverbatim
  76 *>
  77 *> \param[out] T
  78 *> \verbatim
  79 *>          T is DOUBLE PRECISION array, dimension (LDT,MIN(M,N))
  80 *>          The upper triangular block reflectors stored in compact form
  81 *>          as a sequence of upper triangular blocks.  See below
  82 *>          for further details.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] LDT
  86 *> \verbatim
  87 *>          LDT is INTEGER
  88 *>          The leading dimension of the array T.  LDT >= MB.
  89 *> \endverbatim
  90 *>
  91 *> \param[out] WORK
  92 *> \verbatim
  93 *>          WORK is DOUBLE PRECISION array, dimension (MB*N)
  94 *> \endverbatim
  95 *>
  96 *> \param[out] INFO
  97 *> \verbatim
  98 *>          INFO is INTEGER
  99 *>          = 0:  successful exit
 100 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 101 *> \endverbatim
 102 *
 103 *  Authors:
 104 *  ========
 105 *
 106 *> \author Univ. of Tennessee
 107 *> \author Univ. of California Berkeley
 108 *> \author Univ. of Colorado Denver
 109 *> \author NAG Ltd.
 110 *
 111 *> \date November 2013
 112 *
 113 *> \ingroup doubleGEcomputational
 114 *
 115 *> \par Further Details:
 116 *  =====================
 117 *>
 118 *> \verbatim
 119 *>
 120 *>  The matrix V stores the elementary reflectors H(i) in the i-th column
 121 *>  below the diagonal. For example, if M=5 and N=3, the matrix V is
 122 *>
 123 *>               V = (  1  v1 v1 v1 v1 )
 124 *>                   (     1  v2 v2 v2 )
 125 *>                   (         1 v3 v3 )
 126 *>
 127 *>
 128 *>  where the vi's represent the vectors which define H(i), which are returned
 129 *>  in the matrix A.  The 1's along the diagonal of V are not stored in A.
 130 *>  Let K=MIN(M,N).  The number of blocks is B = ceiling(K/NB), where each
 131 *>  block is of order NB except for the last block, which is of order
 132 *>  IB = K - (B-1)*NB.  For each of the B blocks, a upper triangular block
 133 *>  reflector factor is computed: T1, T2, ..., TB.  The NB-by-NB (and IB-by-IB
 134 *>  for the last block) T's are stored in the NB-by-N matrix T as
 135 *>
 136 *>               T = (T1 T2 ... TB).
 137 *> \endverbatim
 138 *>
 139 *  =====================================================================
 140       SUBROUTINE DGELQT( M, N, MB, A, LDA, T, LDT, WORK, INFO )
 141 *
 142 *  -- LAPACK computational routine (version 3.5.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 2013
 146 *
 147 *     .. Scalar Arguments ..
 148       INTEGER INFO, LDA, LDT, M, N, MB
 149 *     ..
 150 *     .. Array Arguments ..
 151       DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
 152 *     ..
 153 *
 154 * =====================================================================
 155 *
 156 *     ..
 157 *     .. Local Scalars ..
 158       INTEGER    I, IB, IINFO, K
 159 *     ..
 160 *     .. External Subroutines ..
 161       EXTERNAL   DGEQRT2, DGEQRT3, DLARFB, XERBLA
 162 *     ..
 163 *     .. Executable Statements ..
 164 *
 165 *     Test the input arguments
 166 *
 167       INFO = 0
 168       IF( M.LT.0 ) THEN
 169          INFO = -1
 170       ELSE IF( N.LT.0 ) THEN
 171          INFO = -2
 172       ELSE IF( MB.LT.1 .OR. ( MB.GT.MIN(M,N) .AND. MIN(M,N).GT.0 ) )THEN
 173          INFO = -3
 174       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 175          INFO = -5
 176       ELSE IF( LDT.LT.MB ) THEN
 177          INFO = -7
 178       END IF
 179       IF( INFO.NE.0 ) THEN
 180          CALL XERBLA( 'DGELQT', -INFO )
 181          RETURN
 182       END IF
 183 *
 184 *     Quick return if possible
 185 *
 186       K = MIN( M, N )
 187       IF( K.EQ.0 ) RETURN
 188 *
 189 *     Blocked loop of length K
 190 *
 191       DO I = 1, K,  MB
 192          IB = MIN( K-I+1, MB )
 193 *
 194 *     Compute the LQ factorization of the current block A(I:M,I:I+IB-1)
 195 *
 196          CALL DGELQT3( IB, N-I+1, A(I,I), LDA, T(1,I), LDT, IINFO )
 197          IF( I+IB.LE.M ) THEN
 198 *
 199 *     Update by applying H**T to A(I:M,I+IB:N) from the right
 200 *
 201          CALL DLARFB( 'R', 'N', 'F', 'R', M-I-IB+1, N-I+1, IB,
 202      $                   A( I, I ), LDA, T( 1, I ), LDT,
 203      $                   A( I+IB, I ), LDA, WORK , M-I-IB+1 )
 204          END IF
 205       END DO
 206       RETURN
 207 *
 208 *     End of DGELQT
 209 *
 210       END