![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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