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