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