 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLABRD( M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, 00002 $ LDY ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER LDA, LDX, LDY, M, N, NB 00011 * .. 00012 * .. Array Arguments .. 00013 REAL D( * ), E( * ) 00014 COMPLEX A( LDA, * ), TAUP( * ), TAUQ( * ), X( LDX, * ), 00015 $ Y( LDY, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CLABRD reduces the first NB rows and columns of a complex general 00022 * m by n matrix A to upper or lower real bidiagonal form by a unitary 00023 * transformation Q**H * A * P, and returns the matrices X and Y which 00024 * are needed to apply the transformation to the unreduced part of A. 00025 * 00026 * If m >= n, A is reduced to upper bidiagonal form; if m < n, to lower 00027 * bidiagonal form. 00028 * 00029 * This is an auxiliary routine called by CGEBRD 00030 * 00031 * Arguments 00032 * ========= 00033 * 00034 * M (input) INTEGER 00035 * The number of rows in the matrix A. 00036 * 00037 * N (input) INTEGER 00038 * The number of columns in the matrix A. 00039 * 00040 * NB (input) INTEGER 00041 * The number of leading rows and columns of A to be reduced. 00042 * 00043 * A (input/output) COMPLEX array, dimension (LDA,N) 00044 * On entry, the m by n general matrix to be reduced. 00045 * On exit, the first NB rows and columns of the matrix are 00046 * overwritten; the rest of the array is unchanged. 00047 * If m >= n, elements on and below the diagonal in the first NB 00048 * columns, with the array TAUQ, represent the unitary 00049 * matrix Q as a product of elementary reflectors; and 00050 * elements above the diagonal in the first NB rows, with the 00051 * array TAUP, represent the unitary matrix P as a product 00052 * of elementary reflectors. 00053 * If m < n, elements below the diagonal in the first NB 00054 * columns, with the array TAUQ, represent the unitary 00055 * matrix Q as a product of elementary reflectors, and 00056 * elements on and above the diagonal in the first NB rows, 00057 * with the array TAUP, represent the unitary matrix P as 00058 * a product of elementary reflectors. 00059 * See Further Details. 00060 * 00061 * LDA (input) INTEGER 00062 * The leading dimension of the array A. LDA >= max(1,M). 00063 * 00064 * D (output) REAL array, dimension (NB) 00065 * The diagonal elements of the first NB rows and columns of 00066 * the reduced matrix. D(i) = A(i,i). 00067 * 00068 * E (output) REAL array, dimension (NB) 00069 * The off-diagonal elements of the first NB rows and columns of 00070 * the reduced matrix. 00071 * 00072 * TAUQ (output) COMPLEX array dimension (NB) 00073 * The scalar factors of the elementary reflectors which 00074 * represent the unitary matrix Q. See Further Details. 00075 * 00076 * TAUP (output) COMPLEX array, dimension (NB) 00077 * The scalar factors of the elementary reflectors which 00078 * represent the unitary matrix P. See Further Details. 00079 * 00080 * X (output) COMPLEX array, dimension (LDX,NB) 00081 * The m-by-nb matrix X required to update the unreduced part 00082 * of A. 00083 * 00084 * LDX (input) INTEGER 00085 * The leading dimension of the array X. LDX >= max(1,M). 00086 * 00087 * Y (output) COMPLEX array, dimension (LDY,NB) 00088 * The n-by-nb matrix Y required to update the unreduced part 00089 * of A. 00090 * 00091 * LDY (input) INTEGER 00092 * The leading dimension of the array Y. LDY >= max(1,N). 00093 * 00094 * Further Details 00095 * =============== 00096 * 00097 * The matrices Q and P are represented as products of elementary 00098 * reflectors: 00099 * 00100 * Q = H(1) H(2) . . . H(nb) and P = G(1) G(2) . . . G(nb) 00101 * 00102 * Each H(i) and G(i) has the form: 00103 * 00104 * H(i) = I - tauq * v * v**H and G(i) = I - taup * u * u**H 00105 * 00106 * where tauq and taup are complex scalars, and v and u are complex 00107 * vectors. 00108 * 00109 * If m >= n, v(1:i-1) = 0, v(i) = 1, and v(i:m) is stored on exit in 00110 * A(i:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+1:n) is stored on exit in 00111 * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). 00112 * 00113 * If m < n, v(1:i) = 0, v(i+1) = 1, and v(i+1:m) is stored on exit in 00114 * A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i:n) is stored on exit in 00115 * A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i). 00116 * 00117 * The elements of the vectors v and u together form the m-by-nb matrix 00118 * V and the nb-by-n matrix U**H which are needed, with X and Y, to apply 00119 * the transformation to the unreduced part of the matrix, using a block 00120 * update of the form: A := A - V*Y**H - X*U**H. 00121 * 00122 * The contents of A on exit are illustrated by the following examples 00123 * with nb = 2: 00124 * 00125 * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n): 00126 * 00127 * ( 1 1 u1 u1 u1 ) ( 1 u1 u1 u1 u1 u1 ) 00128 * ( v1 1 1 u2 u2 ) ( 1 1 u2 u2 u2 u2 ) 00129 * ( v1 v2 a a a ) ( v1 1 a a a a ) 00130 * ( v1 v2 a a a ) ( v1 v2 a a a a ) 00131 * ( v1 v2 a a a ) ( v1 v2 a a a a ) 00132 * ( v1 v2 a a a ) 00133 * 00134 * where a denotes an element of the original matrix which is unchanged, 00135 * vi denotes an element of the vector defining H(i), and ui an element 00136 * of the vector defining G(i). 00137 * 00138 * ===================================================================== 00139 * 00140 * .. Parameters .. 00141 COMPLEX ZERO, ONE 00142 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), 00143 $ ONE = ( 1.0E+0, 0.0E+0 ) ) 00144 * .. 00145 * .. Local Scalars .. 00146 INTEGER I 00147 COMPLEX ALPHA 00148 * .. 00149 * .. External Subroutines .. 00150 EXTERNAL CGEMV, CLACGV, CLARFG, CSCAL 00151 * .. 00152 * .. Intrinsic Functions .. 00153 INTRINSIC MIN 00154 * .. 00155 * .. Executable Statements .. 00156 * 00157 * Quick return if possible 00158 * 00159 IF( M.LE.0 .OR. N.LE.0 ) 00160 $ RETURN 00161 * 00162 IF( M.GE.N ) THEN 00163 * 00164 * Reduce to upper bidiagonal form 00165 * 00166 DO 10 I = 1, NB 00167 * 00168 * Update A(i:m,i) 00169 * 00170 CALL CLACGV( I-1, Y( I, 1 ), LDY ) 00171 CALL CGEMV( 'No transpose', M-I+1, I-1, -ONE, A( I, 1 ), 00172 $ LDA, Y( I, 1 ), LDY, ONE, A( I, I ), 1 ) 00173 CALL CLACGV( I-1, Y( I, 1 ), LDY ) 00174 CALL CGEMV( 'No transpose', M-I+1, I-1, -ONE, X( I, 1 ), 00175 $ LDX, A( 1, I ), 1, ONE, A( I, I ), 1 ) 00176 * 00177 * Generate reflection Q(i) to annihilate A(i+1:m,i) 00178 * 00179 ALPHA = A( I, I ) 00180 CALL CLARFG( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1, 00181 $ TAUQ( I ) ) 00182 D( I ) = ALPHA 00183 IF( I.LT.N ) THEN 00184 A( I, I ) = ONE 00185 * 00186 * Compute Y(i+1:n,i) 00187 * 00188 CALL CGEMV( 'Conjugate transpose', M-I+1, N-I, ONE, 00189 $ A( I, I+1 ), LDA, A( I, I ), 1, ZERO, 00190 $ Y( I+1, I ), 1 ) 00191 CALL CGEMV( 'Conjugate transpose', M-I+1, I-1, ONE, 00192 $ A( I, 1 ), LDA, A( I, I ), 1, ZERO, 00193 $ Y( 1, I ), 1 ) 00194 CALL CGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ), 00195 $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 ) 00196 CALL CGEMV( 'Conjugate transpose', M-I+1, I-1, ONE, 00197 $ X( I, 1 ), LDX, A( I, I ), 1, ZERO, 00198 $ Y( 1, I ), 1 ) 00199 CALL CGEMV( 'Conjugate transpose', I-1, N-I, -ONE, 00200 $ A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE, 00201 $ Y( I+1, I ), 1 ) 00202 CALL CSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 ) 00203 * 00204 * Update A(i,i+1:n) 00205 * 00206 CALL CLACGV( N-I, A( I, I+1 ), LDA ) 00207 CALL CLACGV( I, A( I, 1 ), LDA ) 00208 CALL CGEMV( 'No transpose', N-I, I, -ONE, Y( I+1, 1 ), 00209 $ LDY, A( I, 1 ), LDA, ONE, A( I, I+1 ), LDA ) 00210 CALL CLACGV( I, A( I, 1 ), LDA ) 00211 CALL CLACGV( I-1, X( I, 1 ), LDX ) 00212 CALL CGEMV( 'Conjugate transpose', I-1, N-I, -ONE, 00213 $ A( 1, I+1 ), LDA, X( I, 1 ), LDX, ONE, 00214 $ A( I, I+1 ), LDA ) 00215 CALL CLACGV( I-1, X( I, 1 ), LDX ) 00216 * 00217 * Generate reflection P(i) to annihilate A(i,i+2:n) 00218 * 00219 ALPHA = A( I, I+1 ) 00220 CALL CLARFG( N-I, ALPHA, A( I, MIN( I+2, N ) ), 00221 $ LDA, TAUP( I ) ) 00222 E( I ) = ALPHA 00223 A( I, I+1 ) = ONE 00224 * 00225 * Compute X(i+1:m,i) 00226 * 00227 CALL CGEMV( 'No transpose', M-I, N-I, ONE, A( I+1, I+1 ), 00228 $ LDA, A( I, I+1 ), LDA, ZERO, X( I+1, I ), 1 ) 00229 CALL CGEMV( 'Conjugate transpose', N-I, I, ONE, 00230 $ Y( I+1, 1 ), LDY, A( I, I+1 ), LDA, ZERO, 00231 $ X( 1, I ), 1 ) 00232 CALL CGEMV( 'No transpose', M-I, I, -ONE, A( I+1, 1 ), 00233 $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 00234 CALL CGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ), 00235 $ LDA, A( I, I+1 ), LDA, ZERO, X( 1, I ), 1 ) 00236 CALL CGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ), 00237 $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 00238 CALL CSCAL( M-I, TAUP( I ), X( I+1, I ), 1 ) 00239 CALL CLACGV( N-I, A( I, I+1 ), LDA ) 00240 END IF 00241 10 CONTINUE 00242 ELSE 00243 * 00244 * Reduce to lower bidiagonal form 00245 * 00246 DO 20 I = 1, NB 00247 * 00248 * Update A(i,i:n) 00249 * 00250 CALL CLACGV( N-I+1, A( I, I ), LDA ) 00251 CALL CLACGV( I-1, A( I, 1 ), LDA ) 00252 CALL CGEMV( 'No transpose', N-I+1, I-1, -ONE, Y( I, 1 ), 00253 $ LDY, A( I, 1 ), LDA, ONE, A( I, I ), LDA ) 00254 CALL CLACGV( I-1, A( I, 1 ), LDA ) 00255 CALL CLACGV( I-1, X( I, 1 ), LDX ) 00256 CALL CGEMV( 'Conjugate transpose', I-1, N-I+1, -ONE, 00257 $ A( 1, I ), LDA, X( I, 1 ), LDX, ONE, A( I, I ), 00258 $ LDA ) 00259 CALL CLACGV( I-1, X( I, 1 ), LDX ) 00260 * 00261 * Generate reflection P(i) to annihilate A(i,i+1:n) 00262 * 00263 ALPHA = A( I, I ) 00264 CALL CLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA, 00265 $ TAUP( I ) ) 00266 D( I ) = ALPHA 00267 IF( I.LT.M ) THEN 00268 A( I, I ) = ONE 00269 * 00270 * Compute X(i+1:m,i) 00271 * 00272 CALL CGEMV( 'No transpose', M-I, N-I+1, ONE, A( I+1, I ), 00273 $ LDA, A( I, I ), LDA, ZERO, X( I+1, I ), 1 ) 00274 CALL CGEMV( 'Conjugate transpose', N-I+1, I-1, ONE, 00275 $ Y( I, 1 ), LDY, A( I, I ), LDA, ZERO, 00276 $ X( 1, I ), 1 ) 00277 CALL CGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ), 00278 $ LDA, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 00279 CALL CGEMV( 'No transpose', I-1, N-I+1, ONE, A( 1, I ), 00280 $ LDA, A( I, I ), LDA, ZERO, X( 1, I ), 1 ) 00281 CALL CGEMV( 'No transpose', M-I, I-1, -ONE, X( I+1, 1 ), 00282 $ LDX, X( 1, I ), 1, ONE, X( I+1, I ), 1 ) 00283 CALL CSCAL( M-I, TAUP( I ), X( I+1, I ), 1 ) 00284 CALL CLACGV( N-I+1, A( I, I ), LDA ) 00285 * 00286 * Update A(i+1:m,i) 00287 * 00288 CALL CLACGV( I-1, Y( I, 1 ), LDY ) 00289 CALL CGEMV( 'No transpose', M-I, I-1, -ONE, A( I+1, 1 ), 00290 $ LDA, Y( I, 1 ), LDY, ONE, A( I+1, I ), 1 ) 00291 CALL CLACGV( I-1, Y( I, 1 ), LDY ) 00292 CALL CGEMV( 'No transpose', M-I, I, -ONE, X( I+1, 1 ), 00293 $ LDX, A( 1, I ), 1, ONE, A( I+1, I ), 1 ) 00294 * 00295 * Generate reflection Q(i) to annihilate A(i+2:m,i) 00296 * 00297 ALPHA = A( I+1, I ) 00298 CALL CLARFG( M-I, ALPHA, A( MIN( I+2, M ), I ), 1, 00299 $ TAUQ( I ) ) 00300 E( I ) = ALPHA 00301 A( I+1, I ) = ONE 00302 * 00303 * Compute Y(i+1:n,i) 00304 * 00305 CALL CGEMV( 'Conjugate transpose', M-I, N-I, ONE, 00306 $ A( I+1, I+1 ), LDA, A( I+1, I ), 1, ZERO, 00307 $ Y( I+1, I ), 1 ) 00308 CALL CGEMV( 'Conjugate transpose', M-I, I-1, ONE, 00309 $ A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO, 00310 $ Y( 1, I ), 1 ) 00311 CALL CGEMV( 'No transpose', N-I, I-1, -ONE, Y( I+1, 1 ), 00312 $ LDY, Y( 1, I ), 1, ONE, Y( I+1, I ), 1 ) 00313 CALL CGEMV( 'Conjugate transpose', M-I, I, ONE, 00314 $ X( I+1, 1 ), LDX, A( I+1, I ), 1, ZERO, 00315 $ Y( 1, I ), 1 ) 00316 CALL CGEMV( 'Conjugate transpose', I, N-I, -ONE, 00317 $ A( 1, I+1 ), LDA, Y( 1, I ), 1, ONE, 00318 $ Y( I+1, I ), 1 ) 00319 CALL CSCAL( N-I, TAUQ( I ), Y( I+1, I ), 1 ) 00320 ELSE 00321 CALL CLACGV( N-I+1, A( I, I ), LDA ) 00322 END IF 00323 20 CONTINUE 00324 END IF 00325 RETURN 00326 * 00327 * End of CLABRD 00328 * 00329 END
1.7.3 
