LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhbtrd.f
Go to the documentation of this file.
00001 *> \brief \b ZHBTRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHBTRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbtrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbtrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbtrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
00022 *                          WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO, VECT
00026 *       INTEGER            INFO, KD, LDAB, LDQ, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   D( * ), E( * )
00030 *       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> ZHBTRD reduces a complex Hermitian band matrix A to real symmetric
00040 *> tridiagonal form T by a unitary similarity transformation:
00041 *> Q**H * A * Q = T.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] VECT
00048 *> \verbatim
00049 *>          VECT is CHARACTER*1
00050 *>          = 'N':  do not form Q;
00051 *>          = 'V':  form Q;
00052 *>          = 'U':  update a matrix X, by forming X*Q.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] UPLO
00056 *> \verbatim
00057 *>          UPLO is CHARACTER*1
00058 *>          = 'U':  Upper triangle of A is stored;
00059 *>          = 'L':  Lower triangle of A is stored.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] N
00063 *> \verbatim
00064 *>          N is INTEGER
00065 *>          The order of the matrix A.  N >= 0.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] KD
00069 *> \verbatim
00070 *>          KD is INTEGER
00071 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00072 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00073 *> \endverbatim
00074 *>
00075 *> \param[in,out] AB
00076 *> \verbatim
00077 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
00078 *>          On entry, the upper or lower triangle of the Hermitian band
00079 *>          matrix A, stored in the first KD+1 rows of the array.  The
00080 *>          j-th column of A is stored in the j-th column of the array AB
00081 *>          as follows:
00082 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00083 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00084 *>          On exit, the diagonal elements of AB are overwritten by the
00085 *>          diagonal elements of the tridiagonal matrix T; if KD > 0, the
00086 *>          elements on the first superdiagonal (if UPLO = 'U') or the
00087 *>          first subdiagonal (if UPLO = 'L') are overwritten by the
00088 *>          off-diagonal elements of T; the rest of AB is overwritten by
00089 *>          values generated during the reduction.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDAB
00093 *> \verbatim
00094 *>          LDAB is INTEGER
00095 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00096 *> \endverbatim
00097 *>
00098 *> \param[out] D
00099 *> \verbatim
00100 *>          D is DOUBLE PRECISION array, dimension (N)
00101 *>          The diagonal elements of the tridiagonal matrix T.
00102 *> \endverbatim
00103 *>
00104 *> \param[out] E
00105 *> \verbatim
00106 *>          E is DOUBLE PRECISION array, dimension (N-1)
00107 *>          The off-diagonal elements of the tridiagonal matrix T:
00108 *>          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
00109 *> \endverbatim
00110 *>
00111 *> \param[in,out] Q
00112 *> \verbatim
00113 *>          Q is COMPLEX*16 array, dimension (LDQ,N)
00114 *>          On entry, if VECT = 'U', then Q must contain an N-by-N
00115 *>          matrix X; if VECT = 'N' or 'V', then Q need not be set.
00116 *>
00117 *>          On exit:
00118 *>          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
00119 *>          if VECT = 'U', Q contains the product X*Q;
00120 *>          if VECT = 'N', the array Q is not referenced.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDQ
00124 *> \verbatim
00125 *>          LDQ is INTEGER
00126 *>          The leading dimension of the array Q.
00127 *>          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] WORK
00131 *> \verbatim
00132 *>          WORK is COMPLEX*16 array, dimension (N)
00133 *> \endverbatim
00134 *>
00135 *> \param[out] INFO
00136 *> \verbatim
00137 *>          INFO is INTEGER
00138 *>          = 0:  successful exit
00139 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00140 *> \endverbatim
00141 *
00142 *  Authors:
00143 *  ========
00144 *
00145 *> \author Univ. of Tennessee 
00146 *> \author Univ. of California Berkeley 
00147 *> \author Univ. of Colorado Denver 
00148 *> \author NAG Ltd. 
00149 *
00150 *> \date November 2011
00151 *
00152 *> \ingroup complex16OTHERcomputational
00153 *
00154 *> \par Further Details:
00155 *  =====================
00156 *>
00157 *> \verbatim
00158 *>
00159 *>  Modified by Linda Kaufman, Bell Labs.
00160 *> \endverbatim
00161 *>
00162 *  =====================================================================
00163       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
00164      $                   WORK, INFO )
00165 *
00166 *  -- LAPACK computational routine (version 3.4.0) --
00167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00169 *     November 2011
00170 *
00171 *     .. Scalar Arguments ..
00172       CHARACTER          UPLO, VECT
00173       INTEGER            INFO, KD, LDAB, LDQ, N
00174 *     ..
00175 *     .. Array Arguments ..
00176       DOUBLE PRECISION   D( * ), E( * )
00177       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
00178 *     ..
00179 *
00180 *  =====================================================================
00181 *
00182 *     .. Parameters ..
00183       DOUBLE PRECISION   ZERO
00184       PARAMETER          ( ZERO = 0.0D+0 )
00185       COMPLEX*16         CZERO, CONE
00186       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00187      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00188 *     ..
00189 *     .. Local Scalars ..
00190       LOGICAL            INITQ, UPPER, WANTQ
00191       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
00192      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
00193      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
00194       DOUBLE PRECISION   ABST
00195       COMPLEX*16         T, TEMP
00196 *     ..
00197 *     .. External Subroutines ..
00198       EXTERNAL           XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
00199      $                   ZLASET, ZROT, ZSCAL
00200 *     ..
00201 *     .. Intrinsic Functions ..
00202       INTRINSIC          ABS, DBLE, DCONJG, MAX, MIN
00203 *     ..
00204 *     .. External Functions ..
00205       LOGICAL            LSAME
00206       EXTERNAL           LSAME
00207 *     ..
00208 *     .. Executable Statements ..
00209 *
00210 *     Test the input parameters
00211 *
00212       INITQ = LSAME( VECT, 'V' )
00213       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
00214       UPPER = LSAME( UPLO, 'U' )
00215       KD1 = KD + 1
00216       KDM1 = KD - 1
00217       INCX = LDAB - 1
00218       IQEND = 1
00219 *
00220       INFO = 0
00221       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
00222          INFO = -1
00223       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00224          INFO = -2
00225       ELSE IF( N.LT.0 ) THEN
00226          INFO = -3
00227       ELSE IF( KD.LT.0 ) THEN
00228          INFO = -4
00229       ELSE IF( LDAB.LT.KD1 ) THEN
00230          INFO = -6
00231       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
00232          INFO = -10
00233       END IF
00234       IF( INFO.NE.0 ) THEN
00235          CALL XERBLA( 'ZHBTRD', -INFO )
00236          RETURN
00237       END IF
00238 *
00239 *     Quick return if possible
00240 *
00241       IF( N.EQ.0 )
00242      $   RETURN
00243 *
00244 *     Initialize Q to the unit matrix, if needed
00245 *
00246       IF( INITQ )
00247      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
00248 *
00249 *     Wherever possible, plane rotations are generated and applied in
00250 *     vector operations of length NR over the index set J1:J2:KD1.
00251 *
00252 *     The real cosines and complex sines of the plane rotations are
00253 *     stored in the arrays D and WORK.
00254 *
00255       INCA = KD1*LDAB
00256       KDN = MIN( N-1, KD )
00257       IF( UPPER ) THEN
00258 *
00259          IF( KD.GT.1 ) THEN
00260 *
00261 *           Reduce to complex Hermitian tridiagonal form, working with
00262 *           the upper triangle
00263 *
00264             NR = 0
00265             J1 = KDN + 2
00266             J2 = 1
00267 *
00268             AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
00269             DO 90 I = 1, N - 2
00270 *
00271 *              Reduce i-th row of matrix to tridiagonal form
00272 *
00273                DO 80 K = KDN + 1, 2, -1
00274                   J1 = J1 + KDN
00275                   J2 = J2 + KDN
00276 *
00277                   IF( NR.GT.0 ) THEN
00278 *
00279 *                    generate plane rotations to annihilate nonzero
00280 *                    elements which have been created outside the band
00281 *
00282                      CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
00283      $                            KD1, D( J1 ), KD1 )
00284 *
00285 *                    apply rotations from the right
00286 *
00287 *
00288 *                    Dependent on the the number of diagonals either
00289 *                    ZLARTV or ZROT is used
00290 *
00291                      IF( NR.GE.2*KD-1 ) THEN
00292                         DO 10 L = 1, KD - 1
00293                            CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
00294      $                                  AB( L, J1 ), INCA, D( J1 ),
00295      $                                  WORK( J1 ), KD1 )
00296    10                   CONTINUE
00297 *
00298                      ELSE
00299                         JEND = J1 + ( NR-1 )*KD1
00300                         DO 20 JINC = J1, JEND, KD1
00301                            CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
00302      $                                AB( 1, JINC ), 1, D( JINC ),
00303      $                                WORK( JINC ) )
00304    20                   CONTINUE
00305                      END IF
00306                   END IF
00307 *
00308 *
00309                   IF( K.GT.2 ) THEN
00310                      IF( K.LE.N-I+1 ) THEN
00311 *
00312 *                       generate plane rotation to annihilate a(i,i+k-1)
00313 *                       within the band
00314 *
00315                         CALL ZLARTG( AB( KD-K+3, I+K-2 ),
00316      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
00317      $                               WORK( I+K-1 ), TEMP )
00318                         AB( KD-K+3, I+K-2 ) = TEMP
00319 *
00320 *                       apply rotation from the right
00321 *
00322                         CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
00323      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
00324      $                             WORK( I+K-1 ) )
00325                      END IF
00326                      NR = NR + 1
00327                      J1 = J1 - KDN - 1
00328                   END IF
00329 *
00330 *                 apply plane rotations from both sides to diagonal
00331 *                 blocks
00332 *
00333                   IF( NR.GT.0 )
00334      $               CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
00335      $                            AB( KD, J1 ), INCA, D( J1 ),
00336      $                            WORK( J1 ), KD1 )
00337 *
00338 *                 apply plane rotations from the left
00339 *
00340                   IF( NR.GT.0 ) THEN
00341                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
00342                      IF( 2*KD-1.LT.NR ) THEN
00343 *
00344 *                    Dependent on the the number of diagonals either
00345 *                    ZLARTV or ZROT is used
00346 *
00347                         DO 30 L = 1, KD - 1
00348                            IF( J2+L.GT.N ) THEN
00349                               NRT = NR - 1
00350                            ELSE
00351                               NRT = NR
00352                            END IF
00353                            IF( NRT.GT.0 )
00354      $                        CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
00355      $                                     AB( KD-L+1, J1+L ), INCA,
00356      $                                     D( J1 ), WORK( J1 ), KD1 )
00357    30                   CONTINUE
00358                      ELSE
00359                         J1END = J1 + KD1*( NR-2 )
00360                         IF( J1END.GE.J1 ) THEN
00361                            DO 40 JIN = J1, J1END, KD1
00362                               CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
00363      $                                   AB( KD, JIN+1 ), INCX,
00364      $                                   D( JIN ), WORK( JIN ) )
00365    40                      CONTINUE
00366                         END IF
00367                         LEND = MIN( KDM1, N-J2 )
00368                         LAST = J1END + KD1
00369                         IF( LEND.GT.0 )
00370      $                     CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
00371      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
00372      $                                WORK( LAST ) )
00373                      END IF
00374                   END IF
00375 *
00376                   IF( WANTQ ) THEN
00377 *
00378 *                    accumulate product of plane rotations in Q
00379 *
00380                      IF( INITQ ) THEN
00381 *
00382 *                 take advantage of the fact that Q was
00383 *                 initially the Identity matrix
00384 *
00385                         IQEND = MAX( IQEND, J2 )
00386                         I2 = MAX( 0, K-3 )
00387                         IQAEND = 1 + I*KD
00388                         IF( K.EQ.2 )
00389      $                     IQAEND = IQAEND + KD
00390                         IQAEND = MIN( IQAEND, IQEND )
00391                         DO 50 J = J1, J2, KD1
00392                            IBL = I - I2 / KDM1
00393                            I2 = I2 + 1
00394                            IQB = MAX( 1, J-IBL )
00395                            NQ = 1 + IQAEND - IQB
00396                            IQAEND = MIN( IQAEND+KD, IQEND )
00397                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00398      $                                1, D( J ), DCONJG( WORK( J ) ) )
00399    50                   CONTINUE
00400                      ELSE
00401 *
00402                         DO 60 J = J1, J2, KD1
00403                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00404      $                                D( J ), DCONJG( WORK( J ) ) )
00405    60                   CONTINUE
00406                      END IF
00407 *
00408                   END IF
00409 *
00410                   IF( J2+KDN.GT.N ) THEN
00411 *
00412 *                    adjust J2 to keep within the bounds of the matrix
00413 *
00414                      NR = NR - 1
00415                      J2 = J2 - KDN - 1
00416                   END IF
00417 *
00418                   DO 70 J = J1, J2, KD1
00419 *
00420 *                    create nonzero element a(j-1,j+kd) outside the band
00421 *                    and store it in WORK
00422 *
00423                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
00424                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
00425    70             CONTINUE
00426    80          CONTINUE
00427    90       CONTINUE
00428          END IF
00429 *
00430          IF( KD.GT.0 ) THEN
00431 *
00432 *           make off-diagonal elements real and copy them to E
00433 *
00434             DO 100 I = 1, N - 1
00435                T = AB( KD, I+1 )
00436                ABST = ABS( T )
00437                AB( KD, I+1 ) = ABST
00438                E( I ) = ABST
00439                IF( ABST.NE.ZERO ) THEN
00440                   T = T / ABST
00441                ELSE
00442                   T = CONE
00443                END IF
00444                IF( I.LT.N-1 )
00445      $            AB( KD, I+2 ) = AB( KD, I+2 )*T
00446                IF( WANTQ ) THEN
00447                   CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
00448                END IF
00449   100       CONTINUE
00450          ELSE
00451 *
00452 *           set E to zero if original matrix was diagonal
00453 *
00454             DO 110 I = 1, N - 1
00455                E( I ) = ZERO
00456   110       CONTINUE
00457          END IF
00458 *
00459 *        copy diagonal elements to D
00460 *
00461          DO 120 I = 1, N
00462             D( I ) = AB( KD1, I )
00463   120    CONTINUE
00464 *
00465       ELSE
00466 *
00467          IF( KD.GT.1 ) THEN
00468 *
00469 *           Reduce to complex Hermitian tridiagonal form, working with
00470 *           the lower triangle
00471 *
00472             NR = 0
00473             J1 = KDN + 2
00474             J2 = 1
00475 *
00476             AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
00477             DO 210 I = 1, N - 2
00478 *
00479 *              Reduce i-th column of matrix to tridiagonal form
00480 *
00481                DO 200 K = KDN + 1, 2, -1
00482                   J1 = J1 + KDN
00483                   J2 = J2 + KDN
00484 *
00485                   IF( NR.GT.0 ) THEN
00486 *
00487 *                    generate plane rotations to annihilate nonzero
00488 *                    elements which have been created outside the band
00489 *
00490                      CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
00491      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
00492 *
00493 *                    apply plane rotations from one side
00494 *
00495 *
00496 *                    Dependent on the the number of diagonals either
00497 *                    ZLARTV or ZROT is used
00498 *
00499                      IF( NR.GT.2*KD-1 ) THEN
00500                         DO 130 L = 1, KD - 1
00501                            CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
00502      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
00503      $                                  D( J1 ), WORK( J1 ), KD1 )
00504   130                   CONTINUE
00505                      ELSE
00506                         JEND = J1 + KD1*( NR-1 )
00507                         DO 140 JINC = J1, JEND, KD1
00508                            CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
00509      $                                AB( KD1, JINC-KD ), INCX,
00510      $                                D( JINC ), WORK( JINC ) )
00511   140                   CONTINUE
00512                      END IF
00513 *
00514                   END IF
00515 *
00516                   IF( K.GT.2 ) THEN
00517                      IF( K.LE.N-I+1 ) THEN
00518 *
00519 *                       generate plane rotation to annihilate a(i+k-1,i)
00520 *                       within the band
00521 *
00522                         CALL ZLARTG( AB( K-1, I ), AB( K, I ),
00523      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
00524                         AB( K-1, I ) = TEMP
00525 *
00526 *                       apply rotation from the left
00527 *
00528                         CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
00529      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
00530      $                             WORK( I+K-1 ) )
00531                      END IF
00532                      NR = NR + 1
00533                      J1 = J1 - KDN - 1
00534                   END IF
00535 *
00536 *                 apply plane rotations from both sides to diagonal
00537 *                 blocks
00538 *
00539                   IF( NR.GT.0 )
00540      $               CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
00541      $                            AB( 2, J1-1 ), INCA, D( J1 ),
00542      $                            WORK( J1 ), KD1 )
00543 *
00544 *                 apply plane rotations from the right
00545 *
00546 *
00547 *                    Dependent on the the number of diagonals either
00548 *                    ZLARTV or ZROT is used
00549 *
00550                   IF( NR.GT.0 ) THEN
00551                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
00552                      IF( NR.GT.2*KD-1 ) THEN
00553                         DO 150 L = 1, KD - 1
00554                            IF( J2+L.GT.N ) THEN
00555                               NRT = NR - 1
00556                            ELSE
00557                               NRT = NR
00558                            END IF
00559                            IF( NRT.GT.0 )
00560      $                        CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
00561      $                                     AB( L+1, J1 ), INCA, D( J1 ),
00562      $                                     WORK( J1 ), KD1 )
00563   150                   CONTINUE
00564                      ELSE
00565                         J1END = J1 + KD1*( NR-2 )
00566                         IF( J1END.GE.J1 ) THEN
00567                            DO 160 J1INC = J1, J1END, KD1
00568                               CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
00569      $                                   AB( 2, J1INC ), 1, D( J1INC ),
00570      $                                   WORK( J1INC ) )
00571   160                      CONTINUE
00572                         END IF
00573                         LEND = MIN( KDM1, N-J2 )
00574                         LAST = J1END + KD1
00575                         IF( LEND.GT.0 )
00576      $                     CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
00577      $                                AB( 2, LAST ), 1, D( LAST ),
00578      $                                WORK( LAST ) )
00579                      END IF
00580                   END IF
00581 *
00582 *
00583 *
00584                   IF( WANTQ ) THEN
00585 *
00586 *                    accumulate product of plane rotations in Q
00587 *
00588                      IF( INITQ ) THEN
00589 *
00590 *                 take advantage of the fact that Q was
00591 *                 initially the Identity matrix
00592 *
00593                         IQEND = MAX( IQEND, J2 )
00594                         I2 = MAX( 0, K-3 )
00595                         IQAEND = 1 + I*KD
00596                         IF( K.EQ.2 )
00597      $                     IQAEND = IQAEND + KD
00598                         IQAEND = MIN( IQAEND, IQEND )
00599                         DO 170 J = J1, J2, KD1
00600                            IBL = I - I2 / KDM1
00601                            I2 = I2 + 1
00602                            IQB = MAX( 1, J-IBL )
00603                            NQ = 1 + IQAEND - IQB
00604                            IQAEND = MIN( IQAEND+KD, IQEND )
00605                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
00606      $                                1, D( J ), WORK( J ) )
00607   170                   CONTINUE
00608                      ELSE
00609 *
00610                         DO 180 J = J1, J2, KD1
00611                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
00612      $                                D( J ), WORK( J ) )
00613   180                   CONTINUE
00614                      END IF
00615                   END IF
00616 *
00617                   IF( J2+KDN.GT.N ) THEN
00618 *
00619 *                    adjust J2 to keep within the bounds of the matrix
00620 *
00621                      NR = NR - 1
00622                      J2 = J2 - KDN - 1
00623                   END IF
00624 *
00625                   DO 190 J = J1, J2, KD1
00626 *
00627 *                    create nonzero element a(j+kd,j-1) outside the
00628 *                    band and store it in WORK
00629 *
00630                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
00631                      AB( KD1, J ) = D( J )*AB( KD1, J )
00632   190             CONTINUE
00633   200          CONTINUE
00634   210       CONTINUE
00635          END IF
00636 *
00637          IF( KD.GT.0 ) THEN
00638 *
00639 *           make off-diagonal elements real and copy them to E
00640 *
00641             DO 220 I = 1, N - 1
00642                T = AB( 2, I )
00643                ABST = ABS( T )
00644                AB( 2, I ) = ABST
00645                E( I ) = ABST
00646                IF( ABST.NE.ZERO ) THEN
00647                   T = T / ABST
00648                ELSE
00649                   T = CONE
00650                END IF
00651                IF( I.LT.N-1 )
00652      $            AB( 2, I+1 ) = AB( 2, I+1 )*T
00653                IF( WANTQ ) THEN
00654                   CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
00655                END IF
00656   220       CONTINUE
00657          ELSE
00658 *
00659 *           set E to zero if original matrix was diagonal
00660 *
00661             DO 230 I = 1, N - 1
00662                E( I ) = ZERO
00663   230       CONTINUE
00664          END IF
00665 *
00666 *        copy diagonal elements to D
00667 *
00668          DO 240 I = 1, N
00669             D( I ) = AB( 1, I )
00670   240    CONTINUE
00671       END IF
00672 *
00673       RETURN
00674 *
00675 *     End of ZHBTRD
00676 *
00677       END
 All Files Functions