SRC/slabrd.f

   1 *> \brief \b SLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLABRD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slabrd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slabrd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slabrd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
  22 *                          LDY )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            LDA, LDX, LDY, M, N, NB
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               A( LDA, * ), D( * ), E( * ), TAUP( * ),
  29 *      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLABRD reduces the first NB rows and columns of a real general
  39 *> m by n matrix A to upper or lower bidiagonal form by an orthogonal
  40 *> transformation Q**T * A * P, and returns the matrices X and Y which
  41 *> are needed to apply the transformation to the unreduced part of A.
  42 *>
  43 *> If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
  44 *> bidiagonal form.
  45 *>
  46 *> This is an auxiliary routine called by SGEBRD
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] M
  53 *> \verbatim
  54 *>          M is INTEGER
  55 *>          The number of rows in the matrix A.
  56 *> \endverbatim
  57 *>
  58 *> \param[in] N
  59 *> \verbatim
  60 *>          N is INTEGER
  61 *>          The number of columns in the matrix A.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] NB
  65 *> \verbatim
  66 *>          NB is INTEGER
  67 *>          The number of leading rows and columns of A to be reduced.
  68 *> \endverbatim
  69 *>
  70 *> \param[in,out] A
  71 *> \verbatim
  72 *>          A is REAL array, dimension (LDA,N)
  73 *>          On entry, the m by n general matrix to be reduced.
  74 *>          On exit, the first NB rows and columns of the matrix are
  75 *>          overwritten; the rest of the array is unchanged.
  76 *>          If m >= n, elements on and below the diagonal in the first NB
  77 *>            columns, with the array TAUQ, represent the orthogonal
  78 *>            matrix Q as a product of elementary reflectors; and
  79 *>            elements above the diagonal in the first NB rows, with the
  80 *>            array TAUP, represent the orthogonal matrix P as a product
  81 *>            of elementary reflectors.
  82 *>          If m < n, elements below the diagonal in the first NB
  83 *>            columns, with the array TAUQ, represent the orthogonal
  84 *>            matrix Q as a product of elementary reflectors, and
  85 *>            elements on and above the diagonal in the first NB rows,
  86 *>            with the array TAUP, represent the orthogonal matrix P as
  87 *>            a product of elementary reflectors.
  88 *>          See Further Details.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] LDA
  92 *> \verbatim
  93 *>          LDA is INTEGER
  94 *>          The leading dimension of the array A.  LDA >= max(1,M).
  95 *> \endverbatim
  96 *>
  97 *> \param[out] D
  98 *> \verbatim
  99 *>          D is REAL array, dimension (NB)
 100 *>          The diagonal elements of the first NB rows and columns of
 101 *>          the reduced matrix.  D(i) = A(i,i).
 102 *> \endverbatim
 103 *>
 104 *> \param[out] E
 105 *> \verbatim
 106 *>          E is REAL array, dimension (NB)
 107 *>          The off-diagonal elements of the first NB rows and columns of
 108 *>          the reduced matrix.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] TAUQ
 112 *> \verbatim
 113 *>          TAUQ is REAL array dimension (NB)
 114 *>          The scalar factors of the elementary reflectors which
 115 *>          represent the orthogonal matrix Q. See Further Details.
 116 *> \endverbatim
 117 *>
 118 *> \param[out] TAUP
 119 *> \verbatim
 120 *>          TAUP is REAL array, dimension (NB)
 121 *>          The scalar factors of the elementary reflectors which
 122 *>          represent the orthogonal matrix P. See Further Details.
 123 *> \endverbatim
 124 *>
 125 *> \param[out] X
 126 *> \verbatim
 127 *>          X is REAL array, dimension (LDX,NB)
 128 *>          The m-by-nb matrix X required to update the unreduced part
 129 *>          of A.
 130 *> \endverbatim
 131 *>
 132 *> \param[in] LDX
 133 *> \verbatim
 134 *>          LDX is INTEGER
 135 *>          The leading dimension of the array X. LDX >= max(1,M).
 136 *> \endverbatim
 137 *>
 138 *> \param[out] Y
 139 *> \verbatim
 140 *>          Y is REAL array, dimension (LDY,NB)
 141 *>          The n-by-nb matrix Y required to update the unreduced part
 142 *>          of A.
 143 *> \endverbatim
 144 *>
 145 *> \param[in] LDY
 146 *> \verbatim
 147 *>          LDY is INTEGER
 148 *>          The leading dimension of the array Y. LDY >= max(1,N).
 149 *> \endverbatim
 150 *
 151 *  Authors:
 152 *  ========
 153 *
 154 *> \author Univ. of Tennessee
 155 *> \author Univ. of California Berkeley
 156 *> \author Univ. of Colorado Denver
 157 *> \author NAG Ltd.
 158 *
 159 *> \date September 2012
 160 *
 161 *> \ingroup realOTHERauxiliary
 162 *
 163 *> \par Further Details:
 164 *  =====================
 165 *>
 166 *> \verbatim
 167 *>
 168 *>  The matrices Q and P are represented as products of elementary
 169 *>  reflectors:
 170 *>
 171 *>     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)
 172 *>
 173 *>  Each H(i) and G(i) has the form:
 174 *>
 175 *>     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T
 176 *>
 177 *>  where tauq and taup are real scalars, and v and u are real vectors.
 178 *>
 179 *>  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
 180 *>  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
 181 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
 182 *>
 183 *>  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
 184 *>  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
 185 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
 186 *>
 187 *>  The elements of the vectors v and u together form the m-by-nb matrix
 188 *>  V and the nb-by-n matrix U**T which are needed, with X and Y, to apply
 189 *>  the transformation to the unreduced part of the matrix, using a block
 190 *>  update of the form:  A := A - V*Y**T - X*U**T.
 191 *>
 192 *>  The contents of A on exit are illustrated by the following examples
 193 *>  with nb = 2:
 194 *>
 195 *>  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
 196 *>
 197 *>    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
 198 *>    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
 199 *>    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
 200 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
 201 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
 202 *>    (  v1  v2  a   a   a  )
 203 *>
 204 *>  where a denotes an element of the original matrix which is unchanged,
 205 *>  vi denotes an element of the vector defining H(i), and ui an element
 206 *>  of the vector defining G(i).
 207 *> \endverbatim
 208 *>
 209 *  =====================================================================
 210       SUBROUTINE SLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
 211      $                   LDY )
 212 *
 213 *  -- LAPACK auxiliary routine (version 3.4.2) --
 214 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 215 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 216 *     September 2012
 217 *
 218 *     .. Scalar Arguments ..
 219       INTEGER            LDA, LDX, LDY, M, N, NB
 220 *     ..
 221 *     .. Array Arguments ..
 222       REAL               A( LDA, * ), D( * ), E( * ), TAUP( * ),
 223      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
 224 *     ..
 225 *
 226 *  =====================================================================
 227 *
 228 *     .. Parameters ..
 229       REAL               ZERO, ONE
 230       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
 231 *     ..
 232 *     .. Local Scalars ..
 233       INTEGER            I
 234 *     ..
 235 *     .. External Subroutines ..
 236       EXTERNAL           SGEMV, SLARFG, SSCAL
 237 *     ..
 238 *     .. Intrinsic Functions ..
 239       INTRINSIC          MIN
 240 *     ..
 241 *     .. Executable Statements ..
 242 *
 243 *     Quick return if possible
 244 *
 245       IF( M.LE.0 .OR. N.LE.0 )
 246      $   RETURN
 247 *
 248       IF( M.GE.N ) THEN
 249 *
 250 *        Reduce to upper bidiagonal form
 251 *
 252          DO 10 I = 1, NB
 253 *
 254 *           Update A(i:m,i)
 255 *
 256             CALL SGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
 257      $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
 258             CALL SGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
 259      $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
 260 *
 261 *           Generate reflection Q(i) to annihilate A(i+1:m,i)
 262 *
 263             CALL SLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
 264      $                   TAUQ( I ) )
 265             D( I ) = A( I, I )
 266             IF( I.LT.N ) THEN
 267                A( I, I ) = ONE
 268 *
 269 *              Compute Y(i+1:n,i)
 270 *
 271                CALL SGEMV( 'Transpose', M-I+1, N-I, ONE, A( I, I+1 ),
 272      $                     LDA, A( I, I ), 1, ZERO, Y( I+1, I ), 1 )
 273                CALL SGEMV( 'Transpose', M-I+1, I-1, ONE, A( I, 1 ), LDA,
 274      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
 275                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
 276      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
 277                CALL SGEMV( 'Transpose', M-I+1, I-1, ONE, X( I, 1 ), LDX,
 278      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
 279                CALL SGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
 280      $                     LDA, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
 281                CALL SSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
 282 *
 283 *              Update A(i,i+1:n)
 284 *
 285                CALL SGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
 286      $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
 287                CALL SGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
 288      $                     LDA, X( I, 1 ), LDX, ONE, A( I, I+1 ), LDA )
 289 *
 290 *              Generate reflection P(i) to annihilate A(i,i+2:n)
 291 *
 292                CALL SLARFG( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ),
 293      $                      LDA, TAUP( I ) )
 294                E( I ) = A( I, I+1 )
 295                A( I, I+1 ) = ONE
 296 *
 297 *              Compute X(i+1:m,i)
 298 *
 299                CALL SGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
 300      $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
 301                CALL SGEMV( 'Transpose', N-I, I, ONE, Y( I+1, 1 ), LDY,
 302      $                     A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
 303                CALL SGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
 304      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
 305                CALL SGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
 306      $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
 307                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
 308      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
 309                CALL SSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
 310             END IF
 311    10    CONTINUE
 312       ELSE
 313 *
 314 *        Reduce to lower bidiagonal form
 315 *
 316          DO 20 I = 1, NB
 317 *
 318 *           Update A(i,i:n)
 319 *
 320             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
 321      $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
 322             CALL SGEMV( 'Transpose', I-1, N-I+1, -ONE, A( 1, I ), LDA,
 323      $                  X( I, 1 ), LDX, ONE, A( I, I ), LDA )
 324 *
 325 *           Generate reflection P(i) to annihilate A(i,i+1:n)
 326 *
 327             CALL SLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA,
 328      $                   TAUP( I ) )
 329             D( I ) = A( I, I )
 330             IF( I.LT.M ) THEN
 331                A( I, I ) = ONE
 332 *
 333 *              Compute X(i+1:m,i)
 334 *
 335                CALL SGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
 336      $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
 337                CALL SGEMV( 'Transpose', N-I+1, I-1, ONE, Y( I, 1 ), LDY,
 338      $                     A( I, I ), LDA, ZERO, X( 1, I ), 1 )
 339                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
 340      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
 341                CALL SGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
 342      $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
 343                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
 344      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
 345                CALL SSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
 346 *
 347 *              Update A(i+1:m,i)
 348 *
 349                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
 350      $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
 351                CALL SGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
 352      $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
 353 *
 354 *              Generate reflection Q(i) to annihilate A(i+2:m,i)
 355 *
 356                CALL SLARFG( M-I, A( I+1, I ), A( MIN( I+2, M ), I ), 1,
 357      $                      TAUQ( I ) )
 358                E( I ) = A( I+1, I )
 359                A( I+1, I ) = ONE
 360 *
 361 *              Compute Y(i+1:n,i)
 362 *
 363                CALL SGEMV( 'Transpose', M-I, N-I, ONE, A( I+1, I+1 ),
 364      $                     LDA, A( I+1, I ), 1, ZERO, Y( I+1, I ), 1 )
 365                CALL SGEMV( 'Transpose', M-I, I-1, ONE, A( I+1, 1 ), LDA,
 366      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
 367                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
 368      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
 369                CALL SGEMV( 'Transpose', M-I, I, ONE, X( I+1, 1 ), LDX,
 370      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
 371                CALL SGEMV( 'Transpose', I, N-I, -ONE, A( 1, I+1 ), LDA,
 372      $                     Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
 373                CALL SSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
 374             END IF
 375    20    CONTINUE
 376       END IF
 377       RETURN
 378 *
 379 *     End of SLABRD
 380 *
 381       END