LAPACK: SRC/ssbtrd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SSBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
00002      $                   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          UPLO, VECT
00011       INTEGER            INFO, KD, LDAB, LDQ, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               AB( LDAB, * ), D( * ), E( * ), Q( LDQ, * ),
00015      $                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SSBTRD reduces a real symmetric band matrix A to symmetric
00022 *  tridiagonal form T by an orthogonal similarity transformation:
00023 *  Q**T * A * Q = T.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  VECT    (input) CHARACTER*1
00029 *          = 'N':  do not form Q;
00030 *          = 'V':  form Q;
00031 *          = 'U':  update a matrix X, by forming X*Q.
00032 *
00033 *  UPLO    (input) CHARACTER*1
00034 *          = 'U':  Upper triangle of A is stored;
00035 *          = 'L':  Lower triangle of A is stored.
00036 *
00037 *  N       (input) INTEGER
00038 *          The order of the matrix A.  N >= 0.
00039 *
00040 *  KD      (input) INTEGER
00041 *          The number of superdiagonals of the matrix A if UPLO = 'U',
00042 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00043 *
00044 *  AB      (input/output) REAL array, dimension (LDAB,N)
00045 *          On entry, the upper or lower triangle of the symmetric band
00046 *          matrix A, stored in the first KD+1 rows of the array.  The
00047 *          j-th column of A is stored in the j-th column of the array AB
00048 *          as follows:
00049 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00050 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00051 *          On exit, the diagonal elements of AB are overwritten by the
00052 *          diagonal elements of the tridiagonal matrix T; if KD > 0, the
00053 *          elements on the first superdiagonal (if UPLO = 'U') or the
00054 *          first subdiagonal (if UPLO = 'L') are overwritten by the
00055 *          off-diagonal elements of T; the rest of AB is overwritten by
00056 *          values generated during the reduction.
00057 *
00058 *  LDAB    (input) INTEGER
00059 *          The leading dimension of the array AB.  LDAB >= KD+1.
00060 *
00061 *  D       (output) REAL array, dimension (N)
00062 *          The diagonal elements of the tridiagonal matrix T.
00063 *
00064 *  E       (output) REAL array, dimension (N-1)
00065 *          The off-diagonal elements of the tridiagonal matrix T:
00066 *          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
00067 *
00068 *  Q       (input/output) REAL array, dimension (LDQ,N)
00069 *          On entry, if VECT = 'U', then Q must contain an N-by-N
00070 *          matrix X; if VECT = 'N' or 'V', then Q need not be set.
00071 *
00072 *          On exit:
00073 *          if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
00074 *          if VECT = 'U', Q contains the product X*Q;
00075 *          if VECT = 'N', the array Q is not referenced.
00076 *
00077 *  LDQ     (input) INTEGER
00078 *          The leading dimension of the array Q.
00079 *          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
00080 *
00081 *  WORK    (workspace) REAL array, dimension (N)
00082 *
00083 *  INFO    (output) INTEGER
00084 *          = 0:  successful exit
00085 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00086 *
00087 *  Further Details
00088 *  ===============
00089 *
00090 *  Modified by Linda Kaufman, Bell Labs.
00091 *
00092 *  =====================================================================
00093 *
00094 *     .. Parameters ..
00095       REAL               ZERO, ONE
00096       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00097 *     ..
00098 *     .. Local Scalars ..
00099       LOGICAL            INITQ, UPPER, WANTQ
00100       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
00101      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
00102      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
00103       REAL               TEMP
00104 *     ..
00105 *     .. External Subroutines ..
00106       EXTERNAL           SLAR2V, SLARGV, SLARTG, SLARTV, SLASET, SROT,
00107      $                   XERBLA
00108 *     ..
00109 *     .. Intrinsic Functions ..
00110       INTRINSIC          MAX, MIN
00111 *     ..
00112 *     .. External Functions ..
00113       LOGICAL            LSAME
00114       EXTERNAL           LSAME
00115 *     ..
00116 *     .. Executable Statements ..
00117 *
00118 *     Test the input parameters
00119 *
00120       INITQ = LSAME( VECT, 'V' )
00121       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
00122       UPPER = LSAME( UPLO, 'U' )
00123       KD1 = KD + 1
00124       KDM1 = KD - 1
00125       INCX = LDAB - 1
00126       IQEND = 1
00127 *
00128       INFO = 0
00129       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
00130          INFO = -1
00131       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00132          INFO = -2
00133       ELSE IF( N.LT.0 ) THEN
00134          INFO = -3
00135       ELSE IF( KD.LT.0 ) THEN
00136          INFO = -4
00137       ELSE IF( LDAB.LT.KD1 ) THEN
00138          INFO = -6
00139       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
00140          INFO = -10
00141       END IF
00142       IF( INFO.NE.0 ) THEN
00143          CALL XERBLA( 'SSBTRD', -INFO )
00144          RETURN
00145       END IF
00146 *
00147 *     Quick return if possible
00148 *
00149       IF( N.EQ.0 )
00150      $   RETURN
00151 *
00152 *     Initialize Q to the unit matrix, if needed
00153 *
00154       IF( INITQ )
00155      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00156 *
00157 *     Wherever possible, plane rotations are generated and applied in
00158 *     vector operations of length NR over the index set J1:J2:KD1.
00159 *
00160 *     The cosines and sines of the plane rotations are stored in the
00161 *     arrays D and WORK.
00162 *
00163       INCA = KD1*LDAB
00164       KDN = MIN( N-1, KD )
00165       IF( UPPER ) THEN
00166 *
00167          IF( KD.GT.1 ) THEN
00168 *
00169 *           Reduce to tridiagonal form, working with upper triangle
00170 *
00171             NR = 0
00172             J1 = KDN + 2
00173             J2 = 1
00174 *
00175             DO 90 I = 1, N - 2
00176 *
00177 *              Reduce i-th row of matrix to tridiagonal form
00178 *
00179                DO 80 K = KDN + 1, 2, -1
00180                   J1 = J1 + KDN
00181                   J2 = J2 + KDN
00182 *
00183                   IF( NR.GT.0 ) THEN
00184 *
00185 *                    generate plane rotations to annihilate nonzero
00186 *                    elements which have been created outside the band
00187 *
00188                      CALL SLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
00189      $                            KD1, D( J1 ), KD1 )
00190 *
00191 *                    apply rotations from the right
00192 *
00193 *
00194 *                    Dependent on the the number of diagonals either
00195 *                    SLARTV or SROT is used
00196 *
00197                      IF( NR.GE.2*KD-1 ) THEN
00198                         DO 10 L = 1, KD - 1
00199                            CALL SLARTV( NR, AB( L+1, J1-1 ), INCA,
00200      $                                  AB( L, J1 ), INCA, D( J1 ),
00201      $                                  WORK( J1 ), KD1 )
00202    10                   CONTINUE
00203 *
00204                      ELSE
00205                         JEND = J1 + ( NR-1 )*KD1
00206                         DO 20 JINC = J1, JEND, KD1
00207                            CALL SROT( KDM1, AB( 2, JINC-1 ), 1,
00208      $                                AB( 1, JINC ), 1, D( JINC ),
00209      $                                WORK( JINC ) )
00210    20                   CONTINUE
00211                      END IF
00212                   END IF
00213 *
00214 *
00215                   IF( K.GT.2 ) THEN
00216                      IF( K.LE.N-I+1 ) THEN
00217 *
00218 *                       generate plane rotation to annihilate a(i,i+k-1)
00219 *                       within the band
00220 *
00221                         CALL SLARTG( AB( KD-K+3, I+K-2 ),
00222      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
00223      $                               WORK( I+K-1 ), TEMP )
00224                         AB( KD-K+3, I+K-2 ) = TEMP
00225 *
00226 *                       apply rotation from the right
00227 *
00228                         CALL SROT( K-3, AB( KD-K+4, I+K-2 ), 1,
00229      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
00230      $                             WORK( I+K-1 ) )
00231                      END IF
00232                      NR = NR + 1
00233                      J1 = J1 - KDN - 1
00234                   END IF
00235 *
00236 *                 apply plane rotations from both sides to diagonal
00237 *                 blocks
00238 *
00239                   IF( NR.GT.0 )
00240      $               CALL SLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
00241      $                            AB( KD, J1 ), INCA, D( J1 ),
00242      $                            WORK( J1 ), KD1 )
00243 *
00244 *                 apply plane rotations from the left
00245 *
00246                   IF( NR.GT.0 ) THEN
00247                      IF( 2*KD-1.LT.NR ) THEN
00248 *
00249 *                    Dependent on the the number of diagonals either
00250 *                    SLARTV or SROT is used
00251 *
00252                         DO 30 L = 1, KD - 1
00253                            IF( J2+L.GT.N ) THEN
00254                               NRT = NR - 1
00255                            ELSE
00256                               NRT = NR
00257                            END IF
00258                            IF( NRT.GT.0 )
00259      $                        CALL SLARTV( NRT, AB( KD-L, J1+L ), INCA,
00260      $                                     AB( KD-L+1, J1+L ), INCA,
00261      $                                     D( J1 ), WORK( J1 ), KD1 )
00262    30                   CONTINUE
00263                      ELSE
00264                         J1END = J1 + KD1*( NR-2 )
00265                         IF( J1END.GE.J1 ) THEN
00266                            DO 40 JIN = J1, J1END, KD1
00267                               CALL SROT( KD-1, AB( KD-1, JIN+1 ), INCX,
00268      $                                   AB( KD, JIN+1 ), INCX,
00269      $                                   D( JIN ), WORK( JIN ) )
00270    40                      CONTINUE
00271                         END IF
00272                         LEND = MIN( KDM1, N-J2 )
00273                         LAST = J1END + KD1
00274                         IF( LEND.GT.0 )
00275      $                     CALL SROT( LEND, AB( KD-1, LAST+1 ), INCX,
00276      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
00277      $                                WORK( LAST ) )
00278                      END IF
00279                   END IF
00280 *
00281                   IF( WANTQ ) THEN
00282 *
00283 *                    accumulate product of plane rotations in Q
00284 *
00285                      IF( INITQ ) THEN
00286 *
00287 *                 take advantage of the fact that Q was
00288 *                 initially the Identity matrix
00289 *
00290                         IQEND = MAX( IQEND, J2 )
00291                         I2 = MAX( 0, K-3 )
00292                         IQAEND = 1 + I*KD
00293                         IF( K.EQ.2 )
00294      $                     IQAEND = IQAEND + KD
00295                         IQAEND = MIN( IQAEND, IQEND )
00296                         DO 50 J = J1, J2, KD1
00297                            IBL = I - I2 / KDM1
00298                            I2 = I2 + 1
00299                            IQB = MAX( 1, J-IBL )
00300                            NQ = 1 + IQAEND - IQB
00301                            IQAEND = MIN( IQAEND+KD, IQEND )
00302                            CALL SROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00303      $                                1, D( J ), WORK( J ) )
00304    50                   CONTINUE
00305                      ELSE
00306 *
00307                         DO 60 J = J1, J2, KD1
00308                            CALL SROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00309      $                                D( J ), WORK( J ) )
00310    60                   CONTINUE
00311                      END IF
00312 *
00313                   END IF
00314 *
00315                   IF( J2+KDN.GT.N ) THEN
00316 *
00317 *                    adjust J2 to keep within the bounds of the matrix
00318 *
00319                      NR = NR - 1
00320                      J2 = J2 - KDN - 1
00321                   END IF
00322 *
00323                   DO 70 J = J1, J2, KD1
00324 *
00325 *                    create nonzero element a(j-1,j+kd) outside the band
00326 *                    and store it in WORK
00327 *
00328                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
00329                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
00330    70             CONTINUE
00331    80          CONTINUE
00332    90       CONTINUE
00333          END IF
00334 *
00335          IF( KD.GT.0 ) THEN
00336 *
00337 *           copy off-diagonal elements to E
00338 *
00339             DO 100 I = 1, N - 1
00340                E( I ) = AB( KD, I+1 )
00341   100       CONTINUE
00342          ELSE
00343 *
00344 *           set E to zero if original matrix was diagonal
00345 *
00346             DO 110 I = 1, N - 1
00347                E( I ) = ZERO
00348   110       CONTINUE
00349          END IF
00350 *
00351 *        copy diagonal elements to D
00352 *
00353          DO 120 I = 1, N
00354             D( I ) = AB( KD1, I )
00355   120    CONTINUE
00356 *
00357       ELSE
00358 *
00359          IF( KD.GT.1 ) THEN
00360 *
00361 *           Reduce to tridiagonal form, working with lower triangle
00362 *
00363             NR = 0
00364             J1 = KDN + 2
00365             J2 = 1
00366 *
00367             DO 210 I = 1, N - 2
00368 *
00369 *              Reduce i-th column of matrix to tridiagonal form
00370 *
00371                DO 200 K = KDN + 1, 2, -1
00372                   J1 = J1 + KDN
00373                   J2 = J2 + KDN
00374 *
00375                   IF( NR.GT.0 ) THEN
00376 *
00377 *                    generate plane rotations to annihilate nonzero
00378 *                    elements which have been created outside the band
00379 *
00380                      CALL SLARGV( NR, AB( KD1, J1-KD1 ), INCA,
00381      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
00382 *
00383 *                    apply plane rotations from one side
00384 *
00385 *
00386 *                    Dependent on the the number of diagonals either
00387 *                    SLARTV or SROT is used
00388 *
00389                      IF( NR.GT.2*KD-1 ) THEN
00390                         DO 130 L = 1, KD - 1
00391                            CALL SLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
00392      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
00393      $                                  D( J1 ), WORK( J1 ), KD1 )
00394   130                   CONTINUE
00395                      ELSE
00396                         JEND = J1 + KD1*( NR-1 )
00397                         DO 140 JINC = J1, JEND, KD1
00398                            CALL SROT( KDM1, AB( KD, JINC-KD ), INCX,
00399      $                                AB( KD1, JINC-KD ), INCX,
00400      $                                D( JINC ), WORK( JINC ) )
00401   140                   CONTINUE
00402                      END IF
00403 *
00404                   END IF
00405 *
00406                   IF( K.GT.2 ) THEN
00407                      IF( K.LE.N-I+1 ) THEN
00408 *
00409 *                       generate plane rotation to annihilate a(i+k-1,i)
00410 *                       within the band
00411 *
00412                         CALL SLARTG( AB( K-1, I ), AB( K, I ),
00413      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
00414                         AB( K-1, I ) = TEMP
00415 *
00416 *                       apply rotation from the left
00417 *
00418                         CALL SROT( K-3, AB( K-2, I+1 ), LDAB-1,
00419      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
00420      $                             WORK( I+K-1 ) )
00421                      END IF
00422                      NR = NR + 1
00423                      J1 = J1 - KDN - 1
00424                   END IF
00425 *
00426 *                 apply plane rotations from both sides to diagonal
00427 *                 blocks
00428 *
00429                   IF( NR.GT.0 )
00430      $               CALL SLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
00431      $                            AB( 2, J1-1 ), INCA, D( J1 ),
00432      $                            WORK( J1 ), KD1 )
00433 *
00434 *                 apply plane rotations from the right
00435 *
00436 *
00437 *                    Dependent on the the number of diagonals either
00438 *                    SLARTV or SROT is used
00439 *
00440                   IF( NR.GT.0 ) THEN
00441                      IF( NR.GT.2*KD-1 ) THEN
00442                         DO 150 L = 1, KD - 1
00443                            IF( J2+L.GT.N ) THEN
00444                               NRT = NR - 1
00445                            ELSE
00446                               NRT = NR
00447                            END IF
00448                            IF( NRT.GT.0 )
00449      $                        CALL SLARTV( NRT, AB( L+2, J1-1 ), INCA,
00450      $                                     AB( L+1, J1 ), INCA, D( J1 ),
00451      $                                     WORK( J1 ), KD1 )
00452   150                   CONTINUE
00453                      ELSE
00454                         J1END = J1 + KD1*( NR-2 )
00455                         IF( J1END.GE.J1 ) THEN
00456                            DO 160 J1INC = J1, J1END, KD1
00457                               CALL SROT( KDM1, AB( 3, J1INC-1 ), 1,
00458      $                                   AB( 2, J1INC ), 1, D( J1INC ),
00459      $                                   WORK( J1INC ) )
00460   160                      CONTINUE
00461                         END IF
00462                         LEND = MIN( KDM1, N-J2 )
00463                         LAST = J1END + KD1
00464                         IF( LEND.GT.0 )
00465      $                     CALL SROT( LEND, AB( 3, LAST-1 ), 1,
00466      $                                AB( 2, LAST ), 1, D( LAST ),
00467      $                                WORK( LAST ) )
00468                      END IF
00469                   END IF
00470 *
00471 *
00472 *
00473                   IF( WANTQ ) THEN
00474 *
00475 *                    accumulate product of plane rotations in Q
00476 *
00477                      IF( INITQ ) THEN
00478 *
00479 *                 take advantage of the fact that Q was
00480 *                 initially the Identity matrix
00481 *
00482                         IQEND = MAX( IQEND, J2 )
00483                         I2 = MAX( 0, K-3 )
00484                         IQAEND = 1 + I*KD
00485                         IF( K.EQ.2 )
00486      $                     IQAEND = IQAEND + KD
00487                         IQAEND = MIN( IQAEND, IQEND )
00488                         DO 170 J = J1, J2, KD1
00489                            IBL = I - I2 / KDM1
00490                            I2 = I2 + 1
00491                            IQB = MAX( 1, J-IBL )
00492                            NQ = 1 + IQAEND - IQB
00493                            IQAEND = MIN( IQAEND+KD, IQEND )
00494                            CALL SROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00495      $                                1, D( J ), WORK( J ) )
00496   170                   CONTINUE
00497                      ELSE
00498 *
00499                         DO 180 J = J1, J2, KD1
00500                            CALL SROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00501      $                                D( J ), WORK( J ) )
00502   180                   CONTINUE
00503                      END IF
00504                   END IF
00505 *
00506                   IF( J2+KDN.GT.N ) THEN
00507 *
00508 *                    adjust J2 to keep within the bounds of the matrix
00509 *
00510                      NR = NR - 1
00511                      J2 = J2 - KDN - 1
00512                   END IF
00513 *
00514                   DO 190 J = J1, J2, KD1
00515 *
00516 *                    create nonzero element a(j+kd,j-1) outside the
00517 *                    band and store it in WORK
00518 *
00519                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
00520                      AB( KD1, J ) = D( J )*AB( KD1, J )
00521   190             CONTINUE
00522   200          CONTINUE
00523   210       CONTINUE
00524          END IF
00525 *
00526          IF( KD.GT.0 ) THEN
00527 *
00528 *           copy off-diagonal elements to E
00529 *
00530             DO 220 I = 1, N - 1
00531                E( I ) = AB( 2, I )
00532   220       CONTINUE
00533          ELSE
00534 *
00535 *           set E to zero if original matrix was diagonal
00536 *
00537             DO 230 I = 1, N - 1
00538                E( I ) = ZERO
00539   230       CONTINUE
00540          END IF
00541 *
00542 *        copy diagonal elements to D
00543 *
00544          DO 240 I = 1, N
00545             D( I ) = AB( 1, I )
00546   240    CONTINUE
00547       END IF
00548 *
00549       RETURN
00550 *
00551 *     End of SSBTRD
00552 *
00553       END