LAPACK: SRC/chbtrd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CHBTRD( 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               D( * ), E( * )
00015       COMPLEX            AB( LDAB, * ), Q( LDQ, * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CHBTRD reduces a complex Hermitian band matrix A to real symmetric
00022 *  tridiagonal form T by a unitary similarity transformation:
00023 *  Q**H * 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) COMPLEX array, dimension (LDAB,N)
00045 *          On entry, the upper or lower triangle of the Hermitian 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) COMPLEX 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 unitary 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) COMPLEX 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
00096       PARAMETER          ( ZERO = 0.0E+0 )
00097       COMPLEX            CZERO, CONE
00098       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00099      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00100 *     ..
00101 *     .. Local Scalars ..
00102       LOGICAL            INITQ, UPPER, WANTQ
00103       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
00104      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
00105      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
00106       REAL               ABST
00107       COMPLEX            T, TEMP
00108 *     ..
00109 *     .. External Subroutines ..
00110       EXTERNAL           CLACGV, CLAR2V, CLARGV, CLARTG, CLARTV, CLASET,
00111      $                   CROT, CSCAL, XERBLA
00112 *     ..
00113 *     .. Intrinsic Functions ..
00114       INTRINSIC          ABS, CONJG, MAX, MIN, REAL
00115 *     ..
00116 *     .. External Functions ..
00117       LOGICAL            LSAME
00118       EXTERNAL           LSAME
00119 *     ..
00120 *     .. Executable Statements ..
00121 *
00122 *     Test the input parameters
00123 *
00124       INITQ = LSAME( VECT, 'V' )
00125       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
00126       UPPER = LSAME( UPLO, 'U' )
00127       KD1 = KD + 1
00128       KDM1 = KD - 1
00129       INCX = LDAB - 1
00130       IQEND = 1
00131 *
00132       INFO = 0
00133       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
00134          INFO = -1
00135       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00136          INFO = -2
00137       ELSE IF( N.LT.0 ) THEN
00138          INFO = -3
00139       ELSE IF( KD.LT.0 ) THEN
00140          INFO = -4
00141       ELSE IF( LDAB.LT.KD1 ) THEN
00142          INFO = -6
00143       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
00144          INFO = -10
00145       END IF
00146       IF( INFO.NE.0 ) THEN
00147          CALL XERBLA( 'CHBTRD', -INFO )
00148          RETURN
00149       END IF
00150 *
00151 *     Quick return if possible
00152 *
00153       IF( N.EQ.0 )
00154      $   RETURN
00155 *
00156 *     Initialize Q to the unit matrix, if needed
00157 *
00158       IF( INITQ )
00159      $   CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
00160 *
00161 *     Wherever possible, plane rotations are generated and applied in
00162 *     vector operations of length NR over the index set J1:J2:KD1.
00163 *
00164 *     The real cosines and complex sines of the plane rotations are
00165 *     stored in the arrays D and WORK.
00166 *
00167       INCA = KD1*LDAB
00168       KDN = MIN( N-1, KD )
00169       IF( UPPER ) THEN
00170 *
00171          IF( KD.GT.1 ) THEN
00172 *
00173 *           Reduce to complex Hermitian tridiagonal form, working with
00174 *           the upper triangle
00175 *
00176             NR = 0
00177             J1 = KDN + 2
00178             J2 = 1
00179 *
00180             AB( KD1, 1 ) = REAL( AB( KD1, 1 ) )
00181             DO 90 I = 1, N - 2
00182 *
00183 *              Reduce i-th row of matrix to tridiagonal form
00184 *
00185                DO 80 K = KDN + 1, 2, -1
00186                   J1 = J1 + KDN
00187                   J2 = J2 + KDN
00188 *
00189                   IF( NR.GT.0 ) THEN
00190 *
00191 *                    generate plane rotations to annihilate nonzero
00192 *                    elements which have been created outside the band
00193 *
00194                      CALL CLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
00195      $                            KD1, D( J1 ), KD1 )
00196 *
00197 *                    apply rotations from the right
00198 *
00199 *
00200 *                    Dependent on the the number of diagonals either
00201 *                    CLARTV or CROT is used
00202 *
00203                      IF( NR.GE.2*KD-1 ) THEN
00204                         DO 10 L = 1, KD - 1
00205                            CALL CLARTV( NR, AB( L+1, J1-1 ), INCA,
00206      $                                  AB( L, J1 ), INCA, D( J1 ),
00207      $                                  WORK( J1 ), KD1 )
00208    10                   CONTINUE
00209 *
00210                      ELSE
00211                         JEND = J1 + ( NR-1 )*KD1
00212                         DO 20 JINC = J1, JEND, KD1
00213                            CALL CROT( KDM1, AB( 2, JINC-1 ), 1,
00214      $                                AB( 1, JINC ), 1, D( JINC ),
00215      $                                WORK( JINC ) )
00216    20                   CONTINUE
00217                      END IF
00218                   END IF
00219 *
00220 *
00221                   IF( K.GT.2 ) THEN
00222                      IF( K.LE.N-I+1 ) THEN
00223 *
00224 *                       generate plane rotation to annihilate a(i,i+k-1)
00225 *                       within the band
00226 *
00227                         CALL CLARTG( AB( KD-K+3, I+K-2 ),
00228      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
00229      $                               WORK( I+K-1 ), TEMP )
00230                         AB( KD-K+3, I+K-2 ) = TEMP
00231 *
00232 *                       apply rotation from the right
00233 *
00234                         CALL CROT( K-3, AB( KD-K+4, I+K-2 ), 1,
00235      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
00236      $                             WORK( I+K-1 ) )
00237                      END IF
00238                      NR = NR + 1
00239                      J1 = J1 - KDN - 1
00240                   END IF
00241 *
00242 *                 apply plane rotations from both sides to diagonal
00243 *                 blocks
00244 *
00245                   IF( NR.GT.0 )
00246      $               CALL CLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
00247      $                            AB( KD, J1 ), INCA, D( J1 ),
00248      $                            WORK( J1 ), KD1 )
00249 *
00250 *                 apply plane rotations from the left
00251 *
00252                   IF( NR.GT.0 ) THEN
00253                      CALL CLACGV( NR, WORK( J1 ), KD1 )
00254                      IF( 2*KD-1.LT.NR ) THEN
00255 *
00256 *                    Dependent on the the number of diagonals either
00257 *                    CLARTV or CROT is used
00258 *
00259                         DO 30 L = 1, KD - 1
00260                            IF( J2+L.GT.N ) THEN
00261                               NRT = NR - 1
00262                            ELSE
00263                               NRT = NR
00264                            END IF
00265                            IF( NRT.GT.0 )
00266      $                        CALL CLARTV( NRT, AB( KD-L, J1+L ), INCA,
00267      $                                     AB( KD-L+1, J1+L ), INCA,
00268      $                                     D( J1 ), WORK( J1 ), KD1 )
00269    30                   CONTINUE
00270                      ELSE
00271                         J1END = J1 + KD1*( NR-2 )
00272                         IF( J1END.GE.J1 ) THEN
00273                            DO 40 JIN = J1, J1END, KD1
00274                               CALL CROT( KD-1, AB( KD-1, JIN+1 ), INCX,
00275      $                                   AB( KD, JIN+1 ), INCX,
00276      $                                   D( JIN ), WORK( JIN ) )
00277    40                      CONTINUE
00278                         END IF
00279                         LEND = MIN( KDM1, N-J2 )
00280                         LAST = J1END + KD1
00281                         IF( LEND.GT.0 )
00282      $                     CALL CROT( LEND, AB( KD-1, LAST+1 ), INCX,
00283      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
00284      $                                WORK( LAST ) )
00285                      END IF
00286                   END IF
00287 *
00288                   IF( WANTQ ) THEN
00289 *
00290 *                    accumulate product of plane rotations in Q
00291 *
00292                      IF( INITQ ) THEN
00293 *
00294 *                 take advantage of the fact that Q was
00295 *                 initially the Identity matrix
00296 *
00297                         IQEND = MAX( IQEND, J2 )
00298                         I2 = MAX( 0, K-3 )
00299                         IQAEND = 1 + I*KD
00300                         IF( K.EQ.2 )
00301      $                     IQAEND = IQAEND + KD
00302                         IQAEND = MIN( IQAEND, IQEND )
00303                         DO 50 J = J1, J2, KD1
00304                            IBL = I - I2 / KDM1
00305                            I2 = I2 + 1
00306                            IQB = MAX( 1, J-IBL )
00307                            NQ = 1 + IQAEND - IQB
00308                            IQAEND = MIN( IQAEND+KD, IQEND )
00309                            CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00310      $                                1, D( J ), CONJG( WORK( J ) ) )
00311    50                   CONTINUE
00312                      ELSE
00313 *
00314                         DO 60 J = J1, J2, KD1
00315                            CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00316      $                                D( J ), CONJG( WORK( J ) ) )
00317    60                   CONTINUE
00318                      END IF
00319 *
00320                   END IF
00321 *
00322                   IF( J2+KDN.GT.N ) THEN
00323 *
00324 *                    adjust J2 to keep within the bounds of the matrix
00325 *
00326                      NR = NR - 1
00327                      J2 = J2 - KDN - 1
00328                   END IF
00329 *
00330                   DO 70 J = J1, J2, KD1
00331 *
00332 *                    create nonzero element a(j-1,j+kd) outside the band
00333 *                    and store it in WORK
00334 *
00335                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
00336                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
00337    70             CONTINUE
00338    80          CONTINUE
00339    90       CONTINUE
00340          END IF
00341 *
00342          IF( KD.GT.0 ) THEN
00343 *
00344 *           make off-diagonal elements real and copy them to E
00345 *
00346             DO 100 I = 1, N - 1
00347                T = AB( KD, I+1 )
00348                ABST = ABS( T )
00349                AB( KD, I+1 ) = ABST
00350                E( I ) = ABST
00351                IF( ABST.NE.ZERO ) THEN
00352                   T = T / ABST
00353                ELSE
00354                   T = CONE
00355                END IF
00356                IF( I.LT.N-1 )
00357      $            AB( KD, I+2 ) = AB( KD, I+2 )*T
00358                IF( WANTQ ) THEN
00359                   CALL CSCAL( N, CONJG( T ), Q( 1, I+1 ), 1 )
00360                END IF
00361   100       CONTINUE
00362          ELSE
00363 *
00364 *           set E to zero if original matrix was diagonal
00365 *
00366             DO 110 I = 1, N - 1
00367                E( I ) = ZERO
00368   110       CONTINUE
00369          END IF
00370 *
00371 *        copy diagonal elements to D
00372 *
00373          DO 120 I = 1, N
00374             D( I ) = AB( KD1, I )
00375   120    CONTINUE
00376 *
00377       ELSE
00378 *
00379          IF( KD.GT.1 ) THEN
00380 *
00381 *           Reduce to complex Hermitian tridiagonal form, working with
00382 *           the lower triangle
00383 *
00384             NR = 0
00385             J1 = KDN + 2
00386             J2 = 1
00387 *
00388             AB( 1, 1 ) = REAL( AB( 1, 1 ) )
00389             DO 210 I = 1, N - 2
00390 *
00391 *              Reduce i-th column of matrix to tridiagonal form
00392 *
00393                DO 200 K = KDN + 1, 2, -1
00394                   J1 = J1 + KDN
00395                   J2 = J2 + KDN
00396 *
00397                   IF( NR.GT.0 ) THEN
00398 *
00399 *                    generate plane rotations to annihilate nonzero
00400 *                    elements which have been created outside the band
00401 *
00402                      CALL CLARGV( NR, AB( KD1, J1-KD1 ), INCA,
00403      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
00404 *
00405 *                    apply plane rotations from one side
00406 *
00407 *
00408 *                    Dependent on the the number of diagonals either
00409 *                    CLARTV or CROT is used
00410 *
00411                      IF( NR.GT.2*KD-1 ) THEN
00412                         DO 130 L = 1, KD - 1
00413                            CALL CLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
00414      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
00415      $                                  D( J1 ), WORK( J1 ), KD1 )
00416   130                   CONTINUE
00417                      ELSE
00418                         JEND = J1 + KD1*( NR-1 )
00419                         DO 140 JINC = J1, JEND, KD1
00420                            CALL CROT( KDM1, AB( KD, JINC-KD ), INCX,
00421      $                                AB( KD1, JINC-KD ), INCX,
00422      $                                D( JINC ), WORK( JINC ) )
00423   140                   CONTINUE
00424                      END IF
00425 *
00426                   END IF
00427 *
00428                   IF( K.GT.2 ) THEN
00429                      IF( K.LE.N-I+1 ) THEN
00430 *
00431 *                       generate plane rotation to annihilate a(i+k-1,i)
00432 *                       within the band
00433 *
00434                         CALL CLARTG( AB( K-1, I ), AB( K, I ),
00435      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
00436                         AB( K-1, I ) = TEMP
00437 *
00438 *                       apply rotation from the left
00439 *
00440                         CALL CROT( K-3, AB( K-2, I+1 ), LDAB-1,
00441      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
00442      $                             WORK( I+K-1 ) )
00443                      END IF
00444                      NR = NR + 1
00445                      J1 = J1 - KDN - 1
00446                   END IF
00447 *
00448 *                 apply plane rotations from both sides to diagonal
00449 *                 blocks
00450 *
00451                   IF( NR.GT.0 )
00452      $               CALL CLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
00453      $                            AB( 2, J1-1 ), INCA, D( J1 ),
00454      $                            WORK( J1 ), KD1 )
00455 *
00456 *                 apply plane rotations from the right
00457 *
00458 *
00459 *                    Dependent on the the number of diagonals either
00460 *                    CLARTV or CROT is used
00461 *
00462                   IF( NR.GT.0 ) THEN
00463                      CALL CLACGV( NR, WORK( J1 ), KD1 )
00464                      IF( NR.GT.2*KD-1 ) THEN
00465                         DO 150 L = 1, KD - 1
00466                            IF( J2+L.GT.N ) THEN
00467                               NRT = NR - 1
00468                            ELSE
00469                               NRT = NR
00470                            END IF
00471                            IF( NRT.GT.0 )
00472      $                        CALL CLARTV( NRT, AB( L+2, J1-1 ), INCA,
00473      $                                     AB( L+1, J1 ), INCA, D( J1 ),
00474      $                                     WORK( J1 ), KD1 )
00475   150                   CONTINUE
00476                      ELSE
00477                         J1END = J1 + KD1*( NR-2 )
00478                         IF( J1END.GE.J1 ) THEN
00479                            DO 160 J1INC = J1, J1END, KD1
00480                               CALL CROT( KDM1, AB( 3, J1INC-1 ), 1,
00481      $                                   AB( 2, J1INC ), 1, D( J1INC ),
00482      $                                   WORK( J1INC ) )
00483   160                      CONTINUE
00484                         END IF
00485                         LEND = MIN( KDM1, N-J2 )
00486                         LAST = J1END + KD1
00487                         IF( LEND.GT.0 )
00488      $                     CALL CROT( LEND, AB( 3, LAST-1 ), 1,
00489      $                                AB( 2, LAST ), 1, D( LAST ),
00490      $                                WORK( LAST ) )
00491                      END IF
00492                   END IF
00493 *
00494 *
00495 *
00496                   IF( WANTQ ) THEN
00497 *
00498 *                    accumulate product of plane rotations in Q
00499 *
00500                      IF( INITQ ) THEN
00501 *
00502 *                 take advantage of the fact that Q was
00503 *                 initially the Identity matrix
00504 *
00505                         IQEND = MAX( IQEND, J2 )
00506                         I2 = MAX( 0, K-3 )
00507                         IQAEND = 1 + I*KD
00508                         IF( K.EQ.2 )
00509      $                     IQAEND = IQAEND + KD
00510                         IQAEND = MIN( IQAEND, IQEND )
00511                         DO 170 J = J1, J2, KD1
00512                            IBL = I - I2 / KDM1
00513                            I2 = I2 + 1
00514                            IQB = MAX( 1, J-IBL )
00515                            NQ = 1 + IQAEND - IQB
00516                            IQAEND = MIN( IQAEND+KD, IQEND )
00517                            CALL CROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00518      $                                1, D( J ), WORK( J ) )
00519   170                   CONTINUE
00520                      ELSE
00521 *
00522                         DO 180 J = J1, J2, KD1
00523                            CALL CROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00524      $                                D( J ), WORK( J ) )
00525   180                   CONTINUE
00526                      END IF
00527                   END IF
00528 *
00529                   IF( J2+KDN.GT.N ) THEN
00530 *
00531 *                    adjust J2 to keep within the bounds of the matrix
00532 *
00533                      NR = NR - 1
00534                      J2 = J2 - KDN - 1
00535                   END IF
00536 *
00537                   DO 190 J = J1, J2, KD1
00538 *
00539 *                    create nonzero element a(j+kd,j-1) outside the
00540 *                    band and store it in WORK
00541 *
00542                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
00543                      AB( KD1, J ) = D( J )*AB( KD1, J )
00544   190             CONTINUE
00545   200          CONTINUE
00546   210       CONTINUE
00547          END IF
00548 *
00549          IF( KD.GT.0 ) THEN
00550 *
00551 *           make off-diagonal elements real and copy them to E
00552 *
00553             DO 220 I = 1, N - 1
00554                T = AB( 2, I )
00555                ABST = ABS( T )
00556                AB( 2, I ) = ABST
00557                E( I ) = ABST
00558                IF( ABST.NE.ZERO ) THEN
00559                   T = T / ABST
00560                ELSE
00561                   T = CONE
00562                END IF
00563                IF( I.LT.N-1 )
00564      $            AB( 2, I+1 ) = AB( 2, I+1 )*T
00565                IF( WANTQ ) THEN
00566                   CALL CSCAL( N, T, Q( 1, I+1 ), 1 )
00567                END IF
00568   220       CONTINUE
00569          ELSE
00570 *
00571 *           set E to zero if original matrix was diagonal
00572 *
00573             DO 230 I = 1, N - 1
00574                E( I ) = ZERO
00575   230       CONTINUE
00576          END IF
00577 *
00578 *        copy diagonal elements to D
00579 *
00580          DO 240 I = 1, N
00581             D( I ) = AB( 1, I )
00582   240    CONTINUE
00583       END IF
00584 *
00585       RETURN
00586 *
00587 *     End of CHBTRD
00588 *
00589       END