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