SRC/sgeqrt2.f

   1 *> \brief \b SGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGEQRT2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqrt2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqrt2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqrt2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGEQRT2( M, N, A, LDA, T, LDT, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER   INFO, LDA, LDT, M, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       REAL   A( LDA, * ), T( LDT, * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> SGEQRT2 computes a QR 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 >= N.
  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 REAL array, dimension (LDA,N)
  58 *>          On entry, the real M-by-N matrix A.  On exit, the elements on and
  59 *>          above the diagonal contain the N-by-N upper triangular matrix R; the
  60 *>          elements below the diagonal are the columns of V.  See below for
  61 *>          further details.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] LDA
  65 *> \verbatim
  66 *>          LDA is INTEGER
  67 *>          The leading dimension of the array A.  LDA >= max(1,M).
  68 *> \endverbatim
  69 *>
  70 *> \param[out] T
  71 *> \verbatim
  72 *>          T is REAL array, dimension (LDT,N)
  73 *>          The N-by-N upper triangular factor of the block reflector.
  74 *>          The elements on and above the diagonal contain the block
  75 *>          reflector T; the elements below the diagonal are not used.
  76 *>          See below for further details.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] LDT
  80 *> \verbatim
  81 *>          LDT is INTEGER
  82 *>          The leading dimension of the array T.  LDT >= max(1,N).
  83 *> \endverbatim
  84 *>
  85 *> \param[out] INFO
  86 *> \verbatim
  87 *>          INFO is INTEGER
  88 *>          = 0: successful exit
  89 *>          < 0: if INFO = -i, the i-th argument had an illegal value
  90 *> \endverbatim
  91 *
  92 *  Authors:
  93 *  ========
  94 *
  95 *> \author Univ. of Tennessee
  96 *> \author Univ. of California Berkeley
  97 *> \author Univ. of Colorado Denver
  98 *> \author NAG Ltd.
  99 *
 100 *> \date September 2012
 101 *
 102 *> \ingroup realGEcomputational
 103 *
 104 *> \par Further Details:
 105 *  =====================
 106 *>
 107 *> \verbatim
 108 *>
 109 *>  The matrix V stores the elementary reflectors H(i) in the i-th column
 110 *>  below the diagonal. For example, if M=5 and N=3, the matrix V is
 111 *>
 112 *>               V = (  1       )
 113 *>                   ( v1  1    )
 114 *>                   ( v1 v2  1 )
 115 *>                   ( v1 v2 v3 )
 116 *>                   ( v1 v2 v3 )
 117 *>
 118 *>  where the vi's represent the vectors which define H(i), which are returned
 119 *>  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
 120 *>  block reflector H is then given by
 121 *>
 122 *>               H = I - V * T * V**T
 123 *>
 124 *>  where V**T is the transpose of V.
 125 *> \endverbatim
 126 *>
 127 *  =====================================================================
 128       SUBROUTINE SGEQRT2( M, N, A, LDA, T, LDT, INFO )
 129 *
 130 *  -- LAPACK computational routine (version 3.4.2) --
 131 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 132 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 133 *     September 2012
 134 *
 135 *     .. Scalar Arguments ..
 136       INTEGER   INFO, LDA, LDT, M, N
 137 *     ..
 138 *     .. Array Arguments ..
 139       REAL   A( LDA, * ), T( LDT, * )
 140 *     ..
 141 *
 142 *  =====================================================================
 143 *
 144 *     .. Parameters ..
 145       REAL  ONE, ZERO
 146       PARAMETER( ONE = 1.0, ZERO = 0.0 )
 147 *     ..
 148 *     .. Local Scalars ..
 149       INTEGER   I, K
 150       REAL   AII, ALPHA
 151 *     ..
 152 *     .. External Subroutines ..
 153       EXTERNAL  SLARFG, SGEMV, SGER, STRMV, XERBLA
 154 *     ..
 155 *     .. Executable Statements ..
 156 *
 157 *     Test the input arguments
 158 *
 159       INFO = 0
 160       IF( M.LT.0 ) THEN
 161          INFO = -1
 162       ELSE IF( N.LT.0 ) THEN
 163          INFO = -2
 164       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 165          INFO = -4
 166       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
 167          INFO = -6
 168       END IF
 169       IF( INFO.NE.0 ) THEN
 170          CALL XERBLA( 'SGEQRT2', -INFO )
 171          RETURN
 172       END IF
 173 *
 174       K = MIN( M, N )
 175 *
 176       DO I = 1, K
 177 *
 178 *        Generate elem. refl. H(i) to annihilate A(i+1:m,i), tau(I) -> T(I,1)
 179 *
 180          CALL SLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
 181      $                T( I, 1 ) )
 182          IF( I.LT.N ) THEN
 183 *
 184 *           Apply H(i) to A(I:M,I+1:N) from the left
 185 *
 186             AII = A( I, I )
 187             A( I, I ) = ONE
 188 *
 189 *           W(1:N-I) := A(I:M,I+1:N)^H * A(I:M,I) [W = T(:,N)]
 190 *
 191             CALL SGEMV( 'T',M-I+1, N-I, ONE, A( I, I+1 ), LDA,
 192      $                  A( I, I ), 1, ZERO, T( 1, N ), 1 )
 193 *
 194 *           A(I:M,I+1:N) = A(I:m,I+1:N) + alpha*A(I:M,I)*W(1:N-1)^H
 195 *
 196             ALPHA = -(T( I, 1 ))
 197             CALL SGER( M-I+1, N-I, ALPHA, A( I, I ), 1,
 198      $           T( 1, N ), 1, A( I, I+1 ), LDA )
 199             A( I, I ) = AII
 200          END IF
 201       END DO
 202 *
 203       DO I = 2, N
 204          AII = A( I, I )
 205          A( I, I ) = ONE
 206 *
 207 *        T(1:I-1,I) := alpha * A(I:M,1:I-1)**T * A(I:M,I)
 208 *
 209          ALPHA = -T( I, 1 )
 210          CALL SGEMV( 'T', M-I+1, I-1, ALPHA, A( I, 1 ), LDA,
 211      $               A( I, I ), 1, ZERO, T( 1, I ), 1 )
 212          A( I, I ) = AII
 213 *
 214 *        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
 215 *
 216          CALL STRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
 217 *
 218 *           T(I,I) = tau(I)
 219 *
 220             T( I, I ) = T( I, 1 )
 221             T( I, 1) = ZERO
 222       END DO
 223
 224 *
 225 *     End of SGEQRT2
 226 *
 227       END