|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q, 00002 $ LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER VECT 00011 INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 00015 COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ), 00016 $ Q( LDQ, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZGBBRD reduces a complex general m-by-n band matrix A to real upper 00023 * bidiagonal form B by a unitary transformation: Q' * A * P = B. 00024 * 00025 * The routine computes B, and optionally forms Q or P', or computes 00026 * Q'*C for a given matrix C. 00027 * 00028 * Arguments 00029 * ========= 00030 * 00031 * VECT (input) CHARACTER*1 00032 * Specifies whether or not the matrices Q and P' are to be 00033 * formed. 00034 * = 'N': do not form Q or P'; 00035 * = 'Q': form Q only; 00036 * = 'P': form P' only; 00037 * = 'B': form both. 00038 * 00039 * M (input) INTEGER 00040 * The number of rows of the matrix A. M >= 0. 00041 * 00042 * N (input) INTEGER 00043 * The number of columns of the matrix A. N >= 0. 00044 * 00045 * NCC (input) INTEGER 00046 * The number of columns of the matrix C. NCC >= 0. 00047 * 00048 * KL (input) INTEGER 00049 * The number of subdiagonals of the matrix A. KL >= 0. 00050 * 00051 * KU (input) INTEGER 00052 * The number of superdiagonals of the matrix A. KU >= 0. 00053 * 00054 * AB (input/output) COMPLEX*16 array, dimension (LDAB,N) 00055 * On entry, the m-by-n band matrix A, stored in rows 1 to 00056 * KL+KU+1. The j-th column of A is stored in the j-th column of 00057 * the array AB as follows: 00058 * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl). 00059 * On exit, A is overwritten by values generated during the 00060 * reduction. 00061 * 00062 * LDAB (input) INTEGER 00063 * The leading dimension of the array A. LDAB >= KL+KU+1. 00064 * 00065 * D (output) DOUBLE PRECISION array, dimension (min(M,N)) 00066 * The diagonal elements of the bidiagonal matrix B. 00067 * 00068 * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1) 00069 * The superdiagonal elements of the bidiagonal matrix B. 00070 * 00071 * Q (output) COMPLEX*16 array, dimension (LDQ,M) 00072 * If VECT = 'Q' or 'B', the m-by-m unitary matrix Q. 00073 * If VECT = 'N' or 'P', the array Q is not referenced. 00074 * 00075 * LDQ (input) INTEGER 00076 * The leading dimension of the array Q. 00077 * LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise. 00078 * 00079 * PT (output) COMPLEX*16 array, dimension (LDPT,N) 00080 * If VECT = 'P' or 'B', the n-by-n unitary matrix P'. 00081 * If VECT = 'N' or 'Q', the array PT is not referenced. 00082 * 00083 * LDPT (input) INTEGER 00084 * The leading dimension of the array PT. 00085 * LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise. 00086 * 00087 * C (input/output) COMPLEX*16 array, dimension (LDC,NCC) 00088 * On entry, an m-by-ncc matrix C. 00089 * On exit, C is overwritten by Q'*C. 00090 * C is not referenced if NCC = 0. 00091 * 00092 * LDC (input) INTEGER 00093 * The leading dimension of the array C. 00094 * LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0. 00095 * 00096 * WORK (workspace) COMPLEX*16 array, dimension (max(M,N)) 00097 * 00098 * RWORK (workspace) DOUBLE PRECISION array, dimension (max(M,N)) 00099 * 00100 * INFO (output) INTEGER 00101 * = 0: successful exit. 00102 * < 0: if INFO = -i, the i-th argument had an illegal value. 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 DOUBLE PRECISION ZERO 00108 PARAMETER ( ZERO = 0.0D+0 ) 00109 COMPLEX*16 CZERO, CONE 00110 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00111 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00112 * .. 00113 * .. Local Scalars .. 00114 LOGICAL WANTB, WANTC, WANTPT, WANTQ 00115 INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1, 00116 $ KUN, L, MINMN, ML, ML0, MU, MU0, NR, NRT 00117 DOUBLE PRECISION ABST, RC 00118 COMPLEX*16 RA, RB, RS, T 00119 * .. 00120 * .. External Subroutines .. 00121 EXTERNAL XERBLA, ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT, 00122 $ ZSCAL 00123 * .. 00124 * .. Intrinsic Functions .. 00125 INTRINSIC ABS, DCONJG, MAX, MIN 00126 * .. 00127 * .. External Functions .. 00128 LOGICAL LSAME 00129 EXTERNAL LSAME 00130 * .. 00131 * .. Executable Statements .. 00132 * 00133 * Test the input parameters 00134 * 00135 WANTB = LSAME( VECT, 'B' ) 00136 WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB 00137 WANTPT = LSAME( VECT, 'P' ) .OR. WANTB 00138 WANTC = NCC.GT.0 00139 KLU1 = KL + KU + 1 00140 INFO = 0 00141 IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) ) 00142 $ THEN 00143 INFO = -1 00144 ELSE IF( M.LT.0 ) THEN 00145 INFO = -2 00146 ELSE IF( N.LT.0 ) THEN 00147 INFO = -3 00148 ELSE IF( NCC.LT.0 ) THEN 00149 INFO = -4 00150 ELSE IF( KL.LT.0 ) THEN 00151 INFO = -5 00152 ELSE IF( KU.LT.0 ) THEN 00153 INFO = -6 00154 ELSE IF( LDAB.LT.KLU1 ) THEN 00155 INFO = -8 00156 ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN 00157 INFO = -12 00158 ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN 00159 INFO = -14 00160 ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN 00161 INFO = -16 00162 END IF 00163 IF( INFO.NE.0 ) THEN 00164 CALL XERBLA( 'ZGBBRD', -INFO ) 00165 RETURN 00166 END IF 00167 * 00168 * Initialize Q and P' to the unit matrix, if needed 00169 * 00170 IF( WANTQ ) 00171 $ CALL ZLASET( 'Full', M, M, CZERO, CONE, Q, LDQ ) 00172 IF( WANTPT ) 00173 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, PT, LDPT ) 00174 * 00175 * Quick return if possible. 00176 * 00177 IF( M.EQ.0 .OR. N.EQ.0 ) 00178 $ RETURN 00179 * 00180 MINMN = MIN( M, N ) 00181 * 00182 IF( KL+KU.GT.1 ) THEN 00183 * 00184 * Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce 00185 * first to lower bidiagonal form and then transform to upper 00186 * bidiagonal 00187 * 00188 IF( KU.GT.0 ) THEN 00189 ML0 = 1 00190 MU0 = 2 00191 ELSE 00192 ML0 = 2 00193 MU0 = 1 00194 END IF 00195 * 00196 * Wherever possible, plane rotations are generated and applied in 00197 * vector operations of length NR over the index set J1:J2:KLU1. 00198 * 00199 * The complex sines of the plane rotations are stored in WORK, 00200 * and the real cosines in RWORK. 00201 * 00202 KLM = MIN( M-1, KL ) 00203 KUN = MIN( N-1, KU ) 00204 KB = KLM + KUN 00205 KB1 = KB + 1 00206 INCA = KB1*LDAB 00207 NR = 0 00208 J1 = KLM + 2 00209 J2 = 1 - KUN 00210 * 00211 DO 90 I = 1, MINMN 00212 * 00213 * Reduce i-th column and i-th row of matrix to bidiagonal form 00214 * 00215 ML = KLM + 1 00216 MU = KUN + 1 00217 DO 80 KK = 1, KB 00218 J1 = J1 + KB 00219 J2 = J2 + KB 00220 * 00221 * generate plane rotations to annihilate nonzero elements 00222 * which have been created below the band 00223 * 00224 IF( NR.GT.0 ) 00225 $ CALL ZLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA, 00226 $ WORK( J1 ), KB1, RWORK( J1 ), KB1 ) 00227 * 00228 * apply plane rotations from the left 00229 * 00230 DO 10 L = 1, KB 00231 IF( J2-KLM+L-1.GT.N ) THEN 00232 NRT = NR - 1 00233 ELSE 00234 NRT = NR 00235 END IF 00236 IF( NRT.GT.0 ) 00237 $ CALL ZLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA, 00238 $ AB( KLU1-L+1, J1-KLM+L-1 ), INCA, 00239 $ RWORK( J1 ), WORK( J1 ), KB1 ) 00240 10 CONTINUE 00241 * 00242 IF( ML.GT.ML0 ) THEN 00243 IF( ML.LE.M-I+1 ) THEN 00244 * 00245 * generate plane rotation to annihilate a(i+ml-1,i) 00246 * within the band, and apply rotation from the left 00247 * 00248 CALL ZLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ), 00249 $ RWORK( I+ML-1 ), WORK( I+ML-1 ), RA ) 00250 AB( KU+ML-1, I ) = RA 00251 IF( I.LT.N ) 00252 $ CALL ZROT( MIN( KU+ML-2, N-I ), 00253 $ AB( KU+ML-2, I+1 ), LDAB-1, 00254 $ AB( KU+ML-1, I+1 ), LDAB-1, 00255 $ RWORK( I+ML-1 ), WORK( I+ML-1 ) ) 00256 END IF 00257 NR = NR + 1 00258 J1 = J1 - KB1 00259 END IF 00260 * 00261 IF( WANTQ ) THEN 00262 * 00263 * accumulate product of plane rotations in Q 00264 * 00265 DO 20 J = J1, J2, KB1 00266 CALL ZROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1, 00267 $ RWORK( J ), DCONJG( WORK( J ) ) ) 00268 20 CONTINUE 00269 END IF 00270 * 00271 IF( WANTC ) THEN 00272 * 00273 * apply plane rotations to C 00274 * 00275 DO 30 J = J1, J2, KB1 00276 CALL ZROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC, 00277 $ RWORK( J ), WORK( J ) ) 00278 30 CONTINUE 00279 END IF 00280 * 00281 IF( J2+KUN.GT.N ) THEN 00282 * 00283 * adjust J2 to keep within the bounds of the matrix 00284 * 00285 NR = NR - 1 00286 J2 = J2 - KB1 00287 END IF 00288 * 00289 DO 40 J = J1, J2, KB1 00290 * 00291 * create nonzero element a(j-1,j+ku) above the band 00292 * and store it in WORK(n+1:2*n) 00293 * 00294 WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN ) 00295 AB( 1, J+KUN ) = RWORK( J )*AB( 1, J+KUN ) 00296 40 CONTINUE 00297 * 00298 * generate plane rotations to annihilate nonzero elements 00299 * which have been generated above the band 00300 * 00301 IF( NR.GT.0 ) 00302 $ CALL ZLARGV( NR, AB( 1, J1+KUN-1 ), INCA, 00303 $ WORK( J1+KUN ), KB1, RWORK( J1+KUN ), 00304 $ KB1 ) 00305 * 00306 * apply plane rotations from the right 00307 * 00308 DO 50 L = 1, KB 00309 IF( J2+L-1.GT.M ) THEN 00310 NRT = NR - 1 00311 ELSE 00312 NRT = NR 00313 END IF 00314 IF( NRT.GT.0 ) 00315 $ CALL ZLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA, 00316 $ AB( L, J1+KUN ), INCA, 00317 $ RWORK( J1+KUN ), WORK( J1+KUN ), KB1 ) 00318 50 CONTINUE 00319 * 00320 IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN 00321 IF( MU.LE.N-I+1 ) THEN 00322 * 00323 * generate plane rotation to annihilate a(i,i+mu-1) 00324 * within the band, and apply rotation from the right 00325 * 00326 CALL ZLARTG( AB( KU-MU+3, I+MU-2 ), 00327 $ AB( KU-MU+2, I+MU-1 ), 00328 $ RWORK( I+MU-1 ), WORK( I+MU-1 ), RA ) 00329 AB( KU-MU+3, I+MU-2 ) = RA 00330 CALL ZROT( MIN( KL+MU-2, M-I ), 00331 $ AB( KU-MU+4, I+MU-2 ), 1, 00332 $ AB( KU-MU+3, I+MU-1 ), 1, 00333 $ RWORK( I+MU-1 ), WORK( I+MU-1 ) ) 00334 END IF 00335 NR = NR + 1 00336 J1 = J1 - KB1 00337 END IF 00338 * 00339 IF( WANTPT ) THEN 00340 * 00341 * accumulate product of plane rotations in P' 00342 * 00343 DO 60 J = J1, J2, KB1 00344 CALL ZROT( N, PT( J+KUN-1, 1 ), LDPT, 00345 $ PT( J+KUN, 1 ), LDPT, RWORK( J+KUN ), 00346 $ DCONJG( WORK( J+KUN ) ) ) 00347 60 CONTINUE 00348 END IF 00349 * 00350 IF( J2+KB.GT.M ) THEN 00351 * 00352 * adjust J2 to keep within the bounds of the matrix 00353 * 00354 NR = NR - 1 00355 J2 = J2 - KB1 00356 END IF 00357 * 00358 DO 70 J = J1, J2, KB1 00359 * 00360 * create nonzero element a(j+kl+ku,j+ku-1) below the 00361 * band and store it in WORK(1:n) 00362 * 00363 WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN ) 00364 AB( KLU1, J+KUN ) = RWORK( J+KUN )*AB( KLU1, J+KUN ) 00365 70 CONTINUE 00366 * 00367 IF( ML.GT.ML0 ) THEN 00368 ML = ML - 1 00369 ELSE 00370 MU = MU - 1 00371 END IF 00372 80 CONTINUE 00373 90 CONTINUE 00374 END IF 00375 * 00376 IF( KU.EQ.0 .AND. KL.GT.0 ) THEN 00377 * 00378 * A has been reduced to complex lower bidiagonal form 00379 * 00380 * Transform lower bidiagonal form to upper bidiagonal by applying 00381 * plane rotations from the left, overwriting superdiagonal 00382 * elements on subdiagonal elements 00383 * 00384 DO 100 I = 1, MIN( M-1, N ) 00385 CALL ZLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA ) 00386 AB( 1, I ) = RA 00387 IF( I.LT.N ) THEN 00388 AB( 2, I ) = RS*AB( 1, I+1 ) 00389 AB( 1, I+1 ) = RC*AB( 1, I+1 ) 00390 END IF 00391 IF( WANTQ ) 00392 $ CALL ZROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC, 00393 $ DCONJG( RS ) ) 00394 IF( WANTC ) 00395 $ CALL ZROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC, 00396 $ RS ) 00397 100 CONTINUE 00398 ELSE 00399 * 00400 * A has been reduced to complex upper bidiagonal form or is 00401 * diagonal 00402 * 00403 IF( KU.GT.0 .AND. M.LT.N ) THEN 00404 * 00405 * Annihilate a(m,m+1) by applying plane rotations from the 00406 * right 00407 * 00408 RB = AB( KU, M+1 ) 00409 DO 110 I = M, 1, -1 00410 CALL ZLARTG( AB( KU+1, I ), RB, RC, RS, RA ) 00411 AB( KU+1, I ) = RA 00412 IF( I.GT.1 ) THEN 00413 RB = -DCONJG( RS )*AB( KU, I ) 00414 AB( KU, I ) = RC*AB( KU, I ) 00415 END IF 00416 IF( WANTPT ) 00417 $ CALL ZROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT, 00418 $ RC, DCONJG( RS ) ) 00419 110 CONTINUE 00420 END IF 00421 END IF 00422 * 00423 * Make diagonal and superdiagonal elements real, storing them in D 00424 * and E 00425 * 00426 T = AB( KU+1, 1 ) 00427 DO 120 I = 1, MINMN 00428 ABST = ABS( T ) 00429 D( I ) = ABST 00430 IF( ABST.NE.ZERO ) THEN 00431 T = T / ABST 00432 ELSE 00433 T = CONE 00434 END IF 00435 IF( WANTQ ) 00436 $ CALL ZSCAL( M, T, Q( 1, I ), 1 ) 00437 IF( WANTC ) 00438 $ CALL ZSCAL( NCC, DCONJG( T ), C( I, 1 ), LDC ) 00439 IF( I.LT.MINMN ) THEN 00440 IF( KU.EQ.0 .AND. KL.EQ.0 ) THEN 00441 E( I ) = ZERO 00442 T = AB( 1, I+1 ) 00443 ELSE 00444 IF( KU.EQ.0 ) THEN 00445 T = AB( 2, I )*DCONJG( T ) 00446 ELSE 00447 T = AB( KU, I+1 )*DCONJG( T ) 00448 END IF 00449 ABST = ABS( T ) 00450 E( I ) = ABST 00451 IF( ABST.NE.ZERO ) THEN 00452 T = T / ABST 00453 ELSE 00454 T = CONE 00455 END IF 00456 IF( WANTPT ) 00457 $ CALL ZSCAL( N, T, PT( I+1, 1 ), LDPT ) 00458 T = AB( KU+1, I+1 )*DCONJG( T ) 00459 END IF 00460 END IF 00461 120 CONTINUE 00462 RETURN 00463 * 00464 * End of ZGBBRD 00465 * 00466 END
1.7.3 
