LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slabrd.f
Go to the documentation of this file.
00001 *> \brief \b SLABRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLABRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slabrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slabrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slabrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
00022 *                          LDY )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            LDA, LDX, LDY, M, N, NB
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               A( LDA, * ), D( * ), E( * ), TAUP( * ),
00029 *      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SLABRD reduces the first NB rows and columns of a real general
00039 *> m by n matrix A to upper or lower bidiagonal form by an orthogonal
00040 *> transformation Q**T * A * P, and returns the matrices X and Y which
00041 *> are needed to apply the transformation to the unreduced part of A.
00042 *>
00043 *> If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower
00044 *> bidiagonal form.
00045 *>
00046 *> This is an auxiliary routine called by SGEBRD
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] M
00053 *> \verbatim
00054 *>          M is INTEGER
00055 *>          The number of rows in the matrix A.
00056 *> \endverbatim
00057 *>
00058 *> \param[in] N
00059 *> \verbatim
00060 *>          N is INTEGER
00061 *>          The number of columns in the matrix A.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] NB
00065 *> \verbatim
00066 *>          NB is INTEGER
00067 *>          The number of leading rows and columns of A to be reduced.
00068 *> \endverbatim
00069 *>
00070 *> \param[in,out] A
00071 *> \verbatim
00072 *>          A is REAL array, dimension (LDA,N)
00073 *>          On entry, the m by n general matrix to be reduced.
00074 *>          On exit, the first NB rows and columns of the matrix are
00075 *>          overwritten; the rest of the array is unchanged.
00076 *>          If m >= n, elements on and below the diagonal in the first NB
00077 *>            columns, with the array TAUQ, represent the orthogonal
00078 *>            matrix Q as a product of elementary reflectors; and
00079 *>            elements above the diagonal in the first NB rows, with the
00080 *>            array TAUP, represent the orthogonal matrix P as a product
00081 *>            of elementary reflectors.
00082 *>          If m < n, elements below the diagonal in the first NB
00083 *>            columns, with the array TAUQ, represent the orthogonal
00084 *>            matrix Q as a product of elementary reflectors, and
00085 *>            elements on and above the diagonal in the first NB rows,
00086 *>            with the array TAUP, represent the orthogonal matrix P as
00087 *>            a product of elementary reflectors.
00088 *>          See Further Details.
00089 *> \endverbatim
00090 *>
00091 *> \param[in] LDA
00092 *> \verbatim
00093 *>          LDA is INTEGER
00094 *>          The leading dimension of the array A.  LDA >= max(1,M).
00095 *> \endverbatim
00096 *>
00097 *> \param[out] D
00098 *> \verbatim
00099 *>          D is REAL array, dimension (NB)
00100 *>          The diagonal elements of the first NB rows and columns of
00101 *>          the reduced matrix.  D(i) = A(i,i).
00102 *> \endverbatim
00103 *>
00104 *> \param[out] E
00105 *> \verbatim
00106 *>          E is REAL array, dimension (NB)
00107 *>          The off-diagonal elements of the first NB rows and columns of
00108 *>          the reduced matrix.
00109 *> \endverbatim
00110 *>
00111 *> \param[out] TAUQ
00112 *> \verbatim
00113 *>          TAUQ is REAL array dimension (NB)
00114 *>          The scalar factors of the elementary reflectors which
00115 *>          represent the orthogonal matrix Q. See Further Details.
00116 *> \endverbatim
00117 *>
00118 *> \param[out] TAUP
00119 *> \verbatim
00120 *>          TAUP is REAL array, dimension (NB)
00121 *>          The scalar factors of the elementary reflectors which
00122 *>          represent the orthogonal matrix P. See Further Details.
00123 *> \endverbatim
00124 *>
00125 *> \param[out] X
00126 *> \verbatim
00127 *>          X is REAL array, dimension (LDX,NB)
00128 *>          The m-by-nb matrix X required to update the unreduced part
00129 *>          of A.
00130 *> \endverbatim
00131 *>
00132 *> \param[in] LDX
00133 *> \verbatim
00134 *>          LDX is INTEGER
00135 *>          The leading dimension of the array X. LDX >= max(1,M).
00136 *> \endverbatim
00137 *>
00138 *> \param[out] Y
00139 *> \verbatim
00140 *>          Y is REAL array, dimension (LDY,NB)
00141 *>          The n-by-nb matrix Y required to update the unreduced part
00142 *>          of A.
00143 *> \endverbatim
00144 *>
00145 *> \param[in] LDY
00146 *> \verbatim
00147 *>          LDY is INTEGER
00148 *>          The leading dimension of the array Y. LDY >= max(1,N).
00149 *> \endverbatim
00150 *
00151 *  Authors:
00152 *  ========
00153 *
00154 *> \author Univ. of Tennessee 
00155 *> \author Univ. of California Berkeley 
00156 *> \author Univ. of Colorado Denver 
00157 *> \author NAG Ltd. 
00158 *
00159 *> \date November 2011
00160 *
00161 *> \ingroup realOTHERauxiliary
00162 *
00163 *> \par Further Details:
00164 *  =====================
00165 *>
00166 *> \verbatim
00167 *>
00168 *>  The matrices Q and P are represented as products of elementary
00169 *>  reflectors:
00170 *>
00171 *>     Q = H(1) H(2) . . . H(nb)  and  P = G(1) G(2) . . . G(nb)
00172 *>
00173 *>  Each H(i) and G(i) has the form:
00174 *>
00175 *>     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T
00176 *>
00177 *>  where tauq and taup are real scalars, and v and u are real vectors.
00178 *>
00179 *>  If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in
00180 *>  A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in
00181 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00182 *>
00183 *>  If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in
00184 *>  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in
00185 *>  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
00186 *>
00187 *>  The elements of the vectors v and u together form the m-by-nb matrix
00188 *>  V and the nb-by-n matrix U**T which are needed, with X and Y, to apply
00189 *>  the transformation to the unreduced part of the matrix, using a block
00190 *>  update of the form:  A := A - V*Y**T - X*U**T.
00191 *>
00192 *>  The contents of A on exit are illustrated by the following examples
00193 *>  with nb = 2:
00194 *>
00195 *>  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
00196 *>
00197 *>    (  1   1   u1  u1  u1 )           (  1   u1  u1  u1  u1  u1 )
00198 *>    (  v1  1   1   u2  u2 )           (  1   1   u2  u2  u2  u2 )
00199 *>    (  v1  v2  a   a   a  )           (  v1  1   a   a   a   a  )
00200 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
00201 *>    (  v1  v2  a   a   a  )           (  v1  v2  a   a   a   a  )
00202 *>    (  v1  v2  a   a   a  )
00203 *>
00204 *>  where a denotes an element of the original matrix which is unchanged,
00205 *>  vi denotes an element of the vector defining H(i), and ui an element
00206 *>  of the vector defining G(i).
00207 *> \endverbatim
00208 *>
00209 *  =====================================================================
00210       SUBROUTINE SLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y,
00211      $                   LDY )
00212 *
00213 *  -- LAPACK auxiliary routine (version 3.4.0) --
00214 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00215 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00216 *     November 2011
00217 *
00218 *     .. Scalar Arguments ..
00219       INTEGER            LDA, LDX, LDY, M, N, NB
00220 *     ..
00221 *     .. Array Arguments ..
00222       REAL               A( LDA, * ), D( * ), E( * ), TAUP( * ),
00223      $                   TAUQ( * ), X( LDX, * ), Y( LDY, * )
00224 *     ..
00225 *
00226 *  =====================================================================
00227 *
00228 *     .. Parameters ..
00229       REAL               ZERO, ONE
00230       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00231 *     ..
00232 *     .. Local Scalars ..
00233       INTEGER            I
00234 *     ..
00235 *     .. External Subroutines ..
00236       EXTERNAL           SGEMV, SLARFG, SSCAL
00237 *     ..
00238 *     .. Intrinsic Functions ..
00239       INTRINSIC          MIN
00240 *     ..
00241 *     .. Executable Statements ..
00242 *
00243 *     Quick return if possible
00244 *
00245       IF( M.LE.0 .OR. N.LE.0 )
00246      $   RETURN
00247 *
00248       IF( M.GE.N ) THEN
00249 *
00250 *        Reduce to upper bidiagonal form
00251 *
00252          DO 10 I = 1, NB
00253 *
00254 *           Update A(i:m,i)
00255 *
00256             CALL SGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ),
00257      $                  LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 )
00258             CALL SGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ),
00259      $                  LDX, A( 1, I ), 1, ONE, A( I, I ), 1 )
00260 *
00261 *           Generate reflection Q(i) to annihilate A(i+1:m,i)
00262 *
00263             CALL SLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
00264      $                   TAUQ( I ) )
00265             D( I ) = A( I, I )
00266             IF( I.LT.N ) THEN
00267                A( I, I ) = ONE
00268 *
00269 *              Compute Y(i+1:n,i)
00270 *
00271                CALL SGEMV( 'Transpose', M-I+1, N-I, ONE, A( I, I+1 ),
00272      $                     LDA, A( I, I ), 1, ZERO, Y( I+1, I ), 1 )
00273                CALL SGEMV( 'Transpose', M-I+1, I-1, ONE, A( I, 1 ), LDA,
00274      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
00275                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
00276      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00277                CALL SGEMV( 'Transpose', M-I+1, I-1, ONE, X( I, 1 ), LDX,
00278      $                     A( I, I ), 1, ZERO, Y( 1, I ), 1 )
00279                CALL SGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
00280      $                     LDA, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00281                CALL SSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
00282 *
00283 *              Update A(i,i+1:n)
00284 *
00285                CALL SGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ),
00286      $                     LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA )
00287                CALL SGEMV( 'Transpose', I-1, N-I, -ONE, A( 1, I+1 ),
00288      $                     LDA, X( I, 1 ), LDX, ONE, A( I, I+1 ), LDA )
00289 *
00290 *              Generate reflection P(i) to annihilate A(i,i+2:n)
00291 *
00292                CALL SLARFG( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ),
00293      $                      LDA, TAUP( I ) )
00294                E( I ) = A( I, I+1 )
00295                A( I, I+1 ) = ONE
00296 *
00297 *              Compute X(i+1:m,i)
00298 *
00299                CALL SGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ),
00300      $                     LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 )
00301                CALL SGEMV( 'Transpose', N-I, I, ONE, Y( I+1, 1 ), LDY,
00302      $                     A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
00303                CALL SGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ),
00304      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00305                CALL SGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
00306      $                     LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 )
00307                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
00308      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00309                CALL SSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
00310             END IF
00311    10    CONTINUE
00312       ELSE
00313 *
00314 *        Reduce to lower bidiagonal form
00315 *
00316          DO 20 I = 1, NB
00317 *
00318 *           Update A(i,i:n)
00319 *
00320             CALL SGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ),
00321      $                  LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA )
00322             CALL SGEMV( 'Transpose', I-1, N-I+1, -ONE, A( 1, I ), LDA,
00323      $                  X( I, 1 ), LDX, ONE, A( I, I ), LDA )
00324 *
00325 *           Generate reflection P(i) to annihilate A(i,i+1:n)
00326 *
00327             CALL SLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA,
00328      $                   TAUP( I ) )
00329             D( I ) = A( I, I )
00330             IF( I.LT.M ) THEN
00331                A( I, I ) = ONE
00332 *
00333 *              Compute X(i+1:m,i)
00334 *
00335                CALL SGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ),
00336      $                     LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 )
00337                CALL SGEMV( 'Transpose', N-I+1, I-1, ONE, Y( I, 1 ), LDY,
00338      $                     A( I, I ), LDA, ZERO, X( 1, I ), 1 )
00339                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
00340      $                     LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00341                CALL SGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ),
00342      $                     LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 )
00343                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ),
00344      $                     LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 )
00345                CALL SSCAL( M-I, TAUP( I ), X( I+1, I ), 1 )
00346 *
00347 *              Update A(i+1:m,i)
00348 *
00349                CALL SGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ),
00350      $                     LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 )
00351                CALL SGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ),
00352      $                     LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 )
00353 *
00354 *              Generate reflection Q(i) to annihilate A(i+2:m,i)
00355 *
00356                CALL SLARFG( M-I, A( I+1, I ), A( MIN( I+2, M ), I ), 1,
00357      $                      TAUQ( I ) )
00358                E( I ) = A( I+1, I )
00359                A( I+1, I ) = ONE
00360 *
00361 *              Compute Y(i+1:n,i)
00362 *
00363                CALL SGEMV( 'Transpose', M-I, N-I, ONE, A( I+1, I+1 ),
00364      $                     LDA, A( I+1, I ), 1, ZERO, Y( I+1, I ), 1 )
00365                CALL SGEMV( 'Transpose', M-I, I-1, ONE, A( I+1, 1 ), LDA,
00366      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
00367                CALL SGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ),
00368      $                     LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00369                CALL SGEMV( 'Transpose', M-I, I, ONE, X( I+1, 1 ), LDX,
00370      $                     A( I+1, I ), 1, ZERO, Y( 1, I ), 1 )
00371                CALL SGEMV( 'Transpose', I, N-I, -ONE, A( 1, I+1 ), LDA,
00372      $                     Y( 1, I ), 1, ONE, Y( I+1, I ), 1 )
00373                CALL SSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 )
00374             END IF
00375    20    CONTINUE
00376       END IF
00377       RETURN
00378 *
00379 *     End of SLABRD
00380 *
00381       END
 All Files Functions