File:  [local] / rpl / lapack / lapack / zhbtrd.f
Revision 1.17: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:23 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b ZHBTRD
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZHBTRD + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbtrd.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbtrd.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbtrd.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
   22: *                          WORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          UPLO, VECT
   26: *       INTEGER            INFO, KD, LDAB, LDQ, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   D( * ), E( * )
   30: *       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
   31: *       ..
   32: *
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> ZHBTRD reduces a complex Hermitian band matrix A to real symmetric
   40: *> tridiagonal form T by a unitary similarity transformation:
   41: *> Q**H * A * Q = T.
   42: *> \endverbatim
   43: *
   44: *  Arguments:
   45: *  ==========
   46: *
   47: *> \param[in] VECT
   48: *> \verbatim
   49: *>          VECT is CHARACTER*1
   50: *>          = 'N':  do not form Q;
   51: *>          = 'V':  form Q;
   52: *>          = 'U':  update a matrix X, by forming X*Q.
   53: *> \endverbatim
   54: *>
   55: *> \param[in] UPLO
   56: *> \verbatim
   57: *>          UPLO is CHARACTER*1
   58: *>          = 'U':  Upper triangle of A is stored;
   59: *>          = 'L':  Lower triangle of A is stored.
   60: *> \endverbatim
   61: *>
   62: *> \param[in] N
   63: *> \verbatim
   64: *>          N is INTEGER
   65: *>          The order of the matrix A.  N >= 0.
   66: *> \endverbatim
   67: *>
   68: *> \param[in] KD
   69: *> \verbatim
   70: *>          KD is INTEGER
   71: *>          The number of superdiagonals of the matrix A if UPLO = 'U',
   72: *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
   73: *> \endverbatim
   74: *>
   75: *> \param[in,out] AB
   76: *> \verbatim
   77: *>          AB is COMPLEX*16 array, dimension (LDAB,N)
   78: *>          On entry, the upper or lower triangle of the Hermitian band
   79: *>          matrix A, stored in the first KD+1 rows of the array.  The
   80: *>          j-th column of A is stored in the j-th column of the array AB
   81: *>          as follows:
   82: *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
   83: *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
   84: *>          On exit, the diagonal elements of AB are overwritten by the
   85: *>          diagonal elements of the tridiagonal matrix T; if KD > 0, the
   86: *>          elements on the first superdiagonal (if UPLO = 'U') or the
   87: *>          first subdiagonal (if UPLO = 'L') are overwritten by the
   88: *>          off-diagonal elements of T; the rest of AB is overwritten by
   89: *>          values generated during the reduction.
   90: *> \endverbatim
   91: *>
   92: *> \param[in] LDAB
   93: *> \verbatim
   94: *>          LDAB is INTEGER
   95: *>          The leading dimension of the array AB.  LDAB >= KD+1.
   96: *> \endverbatim
   97: *>
   98: *> \param[out] D
   99: *> \verbatim
  100: *>          D is DOUBLE PRECISION array, dimension (N)
  101: *>          The diagonal elements of the tridiagonal matrix T.
  102: *> \endverbatim
  103: *>
  104: *> \param[out] E
  105: *> \verbatim
  106: *>          E is DOUBLE PRECISION array, dimension (N-1)
  107: *>          The off-diagonal elements of the tridiagonal matrix T:
  108: *>          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
  109: *> \endverbatim
  110: *>
  111: *> \param[in,out] Q
  112: *> \verbatim
  113: *>          Q is COMPLEX*16 array, dimension (LDQ,N)
  114: *>          On entry, if VECT = 'U', then Q must contain an N-by-N
  115: *>          matrix X; if VECT = 'N' or 'V', then Q need not be set.
  116: *>
  117: *>          On exit:
  118: *>          if VECT = 'V', Q contains the N-by-N unitary matrix Q;
  119: *>          if VECT = 'U', Q contains the product X*Q;
  120: *>          if VECT = 'N', the array Q is not referenced.
  121: *> \endverbatim
  122: *>
  123: *> \param[in] LDQ
  124: *> \verbatim
  125: *>          LDQ is INTEGER
  126: *>          The leading dimension of the array Q.
  127: *>          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
  128: *> \endverbatim
  129: *>
  130: *> \param[out] WORK
  131: *> \verbatim
  132: *>          WORK is COMPLEX*16 array, dimension (N)
  133: *> \endverbatim
  134: *>
  135: *> \param[out] INFO
  136: *> \verbatim
  137: *>          INFO is INTEGER
  138: *>          = 0:  successful exit
  139: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  140: *> \endverbatim
  141: *
  142: *  Authors:
  143: *  ========
  144: *
  145: *> \author Univ. of Tennessee
  146: *> \author Univ. of California Berkeley
  147: *> \author Univ. of Colorado Denver
  148: *> \author NAG Ltd.
  149: *
  150: *> \ingroup complex16OTHERcomputational
  151: *
  152: *> \par Further Details:
  153: *  =====================
  154: *>
  155: *> \verbatim
  156: *>
  157: *>  Modified by Linda Kaufman, Bell Labs.
  158: *> \endverbatim
  159: *>
  160: *  =====================================================================
  161:       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
  162:      $                   WORK, INFO )
  163: *
  164: *  -- LAPACK computational routine --
  165: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  166: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  167: *
  168: *     .. Scalar Arguments ..
  169:       CHARACTER          UPLO, VECT
  170:       INTEGER            INFO, KD, LDAB, LDQ, N
  171: *     ..
  172: *     .. Array Arguments ..
  173:       DOUBLE PRECISION   D( * ), E( * )
  174:       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
  175: *     ..
  176: *
  177: *  =====================================================================
  178: *
  179: *     .. Parameters ..
  180:       DOUBLE PRECISION   ZERO
  181:       PARAMETER          ( ZERO = 0.0D+0 )
  182:       COMPLEX*16         CZERO, CONE
  183:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
  184:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  185: *     ..
  186: *     .. Local Scalars ..
  187:       LOGICAL            INITQ, UPPER, WANTQ
  188:       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
  189:      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
  190:      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
  191:       DOUBLE PRECISION   ABST
  192:       COMPLEX*16         T, TEMP
  193: *     ..
  194: *     .. External Subroutines ..
  195:       EXTERNAL           XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
  196:      $                   ZLASET, ZROT, ZSCAL
  197: *     ..
  198: *     .. Intrinsic Functions ..
  199:       INTRINSIC          ABS, DBLE, DCONJG, MAX, MIN
  200: *     ..
  201: *     .. External Functions ..
  202:       LOGICAL            LSAME
  203:       EXTERNAL           LSAME
  204: *     ..
  205: *     .. Executable Statements ..
  206: *
  207: *     Test the input parameters
  208: *
  209:       INITQ = LSAME( VECT, 'V' )
  210:       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
  211:       UPPER = LSAME( UPLO, 'U' )
  212:       KD1 = KD + 1
  213:       KDM1 = KD - 1
  214:       INCX = LDAB - 1
  215:       IQEND = 1
  216: *
  217:       INFO = 0
  218:       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
  219:          INFO = -1
  220:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  221:          INFO = -2
  222:       ELSE IF( N.LT.0 ) THEN
  223:          INFO = -3
  224:       ELSE IF( KD.LT.0 ) THEN
  225:          INFO = -4
  226:       ELSE IF( LDAB.LT.KD1 ) THEN
  227:          INFO = -6
  228:       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
  229:          INFO = -10
  230:       END IF
  231:       IF( INFO.NE.0 ) THEN
  232:          CALL XERBLA( 'ZHBTRD', -INFO )
  233:          RETURN
  234:       END IF
  235: *
  236: *     Quick return if possible
  237: *
  238:       IF( N.EQ.0 )
  239:      $   RETURN
  240: *
  241: *     Initialize Q to the unit matrix, if needed
  242: *
  243:       IF( INITQ )
  244:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
  245: *
  246: *     Wherever possible, plane rotations are generated and applied in
  247: *     vector operations of length NR over the index set J1:J2:KD1.
  248: *
  249: *     The real cosines and complex sines of the plane rotations are
  250: *     stored in the arrays D and WORK.
  251: *
  252:       INCA = KD1*LDAB
  253:       KDN = MIN( N-1, KD )
  254:       IF( UPPER ) THEN
  255: *
  256:          IF( KD.GT.1 ) THEN
  257: *
  258: *           Reduce to complex Hermitian tridiagonal form, working with
  259: *           the upper triangle
  260: *
  261:             NR = 0
  262:             J1 = KDN + 2
  263:             J2 = 1
  264: *
  265:             AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
  266:             DO 90 I = 1, N - 2
  267: *
  268: *              Reduce i-th row of matrix to tridiagonal form
  269: *
  270:                DO 80 K = KDN + 1, 2, -1
  271:                   J1 = J1 + KDN
  272:                   J2 = J2 + KDN
  273: *
  274:                   IF( NR.GT.0 ) THEN
  275: *
  276: *                    generate plane rotations to annihilate nonzero
  277: *                    elements which have been created outside the band
  278: *
  279:                      CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
  280:      $                            KD1, D( J1 ), KD1 )
  281: *
  282: *                    apply rotations from the right
  283: *
  284: *
  285: *                    Dependent on the the number of diagonals either
  286: *                    ZLARTV or ZROT is used
  287: *
  288:                      IF( NR.GE.2*KD-1 ) THEN
  289:                         DO 10 L = 1, KD - 1
  290:                            CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
  291:      $                                  AB( L, J1 ), INCA, D( J1 ),
  292:      $                                  WORK( J1 ), KD1 )
  293:    10                   CONTINUE
  294: *
  295:                      ELSE
  296:                         JEND = J1 + ( NR-1 )*KD1
  297:                         DO 20 JINC = J1, JEND, KD1
  298:                            CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
  299:      $                                AB( 1, JINC ), 1, D( JINC ),
  300:      $                                WORK( JINC ) )
  301:    20                   CONTINUE
  302:                      END IF
  303:                   END IF
  304: *
  305: *
  306:                   IF( K.GT.2 ) THEN
  307:                      IF( K.LE.N-I+1 ) THEN
  308: *
  309: *                       generate plane rotation to annihilate a(i,i+k-1)
  310: *                       within the band
  311: *
  312:                         CALL ZLARTG( AB( KD-K+3, I+K-2 ),
  313:      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
  314:      $                               WORK( I+K-1 ), TEMP )
  315:                         AB( KD-K+3, I+K-2 ) = TEMP
  316: *
  317: *                       apply rotation from the right
  318: *
  319:                         CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
  320:      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
  321:      $                             WORK( I+K-1 ) )
  322:                      END IF
  323:                      NR = NR + 1
  324:                      J1 = J1 - KDN - 1
  325:                   END IF
  326: *
  327: *                 apply plane rotations from both sides to diagonal
  328: *                 blocks
  329: *
  330:                   IF( NR.GT.0 )
  331:      $               CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
  332:      $                            AB( KD, J1 ), INCA, D( J1 ),
  333:      $                            WORK( J1 ), KD1 )
  334: *
  335: *                 apply plane rotations from the left
  336: *
  337:                   IF( NR.GT.0 ) THEN
  338:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
  339:                      IF( 2*KD-1.LT.NR ) THEN
  340: *
  341: *                    Dependent on the the number of diagonals either
  342: *                    ZLARTV or ZROT is used
  343: *
  344:                         DO 30 L = 1, KD - 1
  345:                            IF( J2+L.GT.N ) THEN
  346:                               NRT = NR - 1
  347:                            ELSE
  348:                               NRT = NR
  349:                            END IF
  350:                            IF( NRT.GT.0 )
  351:      $                        CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
  352:      $                                     AB( KD-L+1, J1+L ), INCA,
  353:      $                                     D( J1 ), WORK( J1 ), KD1 )
  354:    30                   CONTINUE
  355:                      ELSE
  356:                         J1END = J1 + KD1*( NR-2 )
  357:                         IF( J1END.GE.J1 ) THEN
  358:                            DO 40 JIN = J1, J1END, KD1
  359:                               CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
  360:      $                                   AB( KD, JIN+1 ), INCX,
  361:      $                                   D( JIN ), WORK( JIN ) )
  362:    40                      CONTINUE
  363:                         END IF
  364:                         LEND = MIN( KDM1, N-J2 )
  365:                         LAST = J1END + KD1
  366:                         IF( LEND.GT.0 )
  367:      $                     CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
  368:      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
  369:      $                                WORK( LAST ) )
  370:                      END IF
  371:                   END IF
  372: *
  373:                   IF( WANTQ ) THEN
  374: *
  375: *                    accumulate product of plane rotations in Q
  376: *
  377:                      IF( INITQ ) THEN
  378: *
  379: *                 take advantage of the fact that Q was
  380: *                 initially the Identity matrix
  381: *
  382:                         IQEND = MAX( IQEND, J2 )
  383:                         I2 = MAX( 0, K-3 )
  384:                         IQAEND = 1 + I*KD
  385:                         IF( K.EQ.2 )
  386:      $                     IQAEND = IQAEND + KD
  387:                         IQAEND = MIN( IQAEND, IQEND )
  388:                         DO 50 J = J1, J2, KD1
  389:                            IBL = I - I2 / KDM1
  390:                            I2 = I2 + 1
  391:                            IQB = MAX( 1, J-IBL )
  392:                            NQ = 1 + IQAEND - IQB
  393:                            IQAEND = MIN( IQAEND+KD, IQEND )
  394:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
  395:      $                                1, D( J ), DCONJG( WORK( J ) ) )
  396:    50                   CONTINUE
  397:                      ELSE
  398: *
  399:                         DO 60 J = J1, J2, KD1
  400:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
  401:      $                                D( J ), DCONJG( WORK( J ) ) )
  402:    60                   CONTINUE
  403:                      END IF
  404: *
  405:                   END IF
  406: *
  407:                   IF( J2+KDN.GT.N ) THEN
  408: *
  409: *                    adjust J2 to keep within the bounds of the matrix
  410: *
  411:                      NR = NR - 1
  412:                      J2 = J2 - KDN - 1
  413:                   END IF
  414: *
  415:                   DO 70 J = J1, J2, KD1
  416: *
  417: *                    create nonzero element a(j-1,j+kd) outside the band
  418: *                    and store it in WORK
  419: *
  420:                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
  421:                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
  422:    70             CONTINUE
  423:    80          CONTINUE
  424:    90       CONTINUE
  425:          END IF
  426: *
  427:          IF( KD.GT.0 ) THEN
  428: *
  429: *           make off-diagonal elements real and copy them to E
  430: *
  431:             DO 100 I = 1, N - 1
  432:                T = AB( KD, I+1 )
  433:                ABST = ABS( T )
  434:                AB( KD, I+1 ) = ABST
  435:                E( I ) = ABST
  436:                IF( ABST.NE.ZERO ) THEN
  437:                   T = T / ABST
  438:                ELSE
  439:                   T = CONE
  440:                END IF
  441:                IF( I.LT.N-1 )
  442:      $            AB( KD, I+2 ) = AB( KD, I+2 )*T
  443:                IF( WANTQ ) THEN
  444:                   CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
  445:                END IF
  446:   100       CONTINUE
  447:          ELSE
  448: *
  449: *           set E to zero if original matrix was diagonal
  450: *
  451:             DO 110 I = 1, N - 1
  452:                E( I ) = ZERO
  453:   110       CONTINUE
  454:          END IF
  455: *
  456: *        copy diagonal elements to D
  457: *
  458:          DO 120 I = 1, N
  459:             D( I ) = DBLE( AB( KD1, I ) )
  460:   120    CONTINUE
  461: *
  462:       ELSE
  463: *
  464:          IF( KD.GT.1 ) THEN
  465: *
  466: *           Reduce to complex Hermitian tridiagonal form, working with
  467: *           the lower triangle
  468: *
  469:             NR = 0
  470:             J1 = KDN + 2
  471:             J2 = 1
  472: *
  473:             AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
  474:             DO 210 I = 1, N - 2
  475: *
  476: *              Reduce i-th column of matrix to tridiagonal form
  477: *
  478:                DO 200 K = KDN + 1, 2, -1
  479:                   J1 = J1 + KDN
  480:                   J2 = J2 + KDN
  481: *
  482:                   IF( NR.GT.0 ) THEN
  483: *
  484: *                    generate plane rotations to annihilate nonzero
  485: *                    elements which have been created outside the band
  486: *
  487:                      CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
  488:      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
  489: *
  490: *                    apply plane rotations from one side
  491: *
  492: *
  493: *                    Dependent on the the number of diagonals either
  494: *                    ZLARTV or ZROT is used
  495: *
  496:                      IF( NR.GT.2*KD-1 ) THEN
  497:                         DO 130 L = 1, KD - 1
  498:                            CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
  499:      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
  500:      $                                  D( J1 ), WORK( J1 ), KD1 )
  501:   130                   CONTINUE
  502:                      ELSE
  503:                         JEND = J1 + KD1*( NR-1 )
  504:                         DO 140 JINC = J1, JEND, KD1
  505:                            CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
  506:      $                                AB( KD1, JINC-KD ), INCX,
  507:      $                                D( JINC ), WORK( JINC ) )
  508:   140                   CONTINUE
  509:                      END IF
  510: *
  511:                   END IF
  512: *
  513:                   IF( K.GT.2 ) THEN
  514:                      IF( K.LE.N-I+1 ) THEN
  515: *
  516: *                       generate plane rotation to annihilate a(i+k-1,i)
  517: *                       within the band
  518: *
  519:                         CALL ZLARTG( AB( K-1, I ), AB( K, I ),
  520:      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
  521:                         AB( K-1, I ) = TEMP
  522: *
  523: *                       apply rotation from the left
  524: *
  525:                         CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
  526:      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
  527:      $                             WORK( I+K-1 ) )
  528:                      END IF
  529:                      NR = NR + 1
  530:                      J1 = J1 - KDN - 1
  531:                   END IF
  532: *
  533: *                 apply plane rotations from both sides to diagonal
  534: *                 blocks
  535: *
  536:                   IF( NR.GT.0 )
  537:      $               CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
  538:      $                            AB( 2, J1-1 ), INCA, D( J1 ),
  539:      $                            WORK( J1 ), KD1 )
  540: *
  541: *                 apply plane rotations from the right
  542: *
  543: *
  544: *                    Dependent on the the number of diagonals either
  545: *                    ZLARTV or ZROT is used
  546: *
  547:                   IF( NR.GT.0 ) THEN
  548:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
  549:                      IF( NR.GT.2*KD-1 ) THEN
  550:                         DO 150 L = 1, KD - 1
  551:                            IF( J2+L.GT.N ) THEN
  552:                               NRT = NR - 1
  553:                            ELSE
  554:                               NRT = NR
  555:                            END IF
  556:                            IF( NRT.GT.0 )
  557:      $                        CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
  558:      $                                     AB( L+1, J1 ), INCA, D( J1 ),
  559:      $                                     WORK( J1 ), KD1 )
  560:   150                   CONTINUE
  561:                      ELSE
  562:                         J1END = J1 + KD1*( NR-2 )
  563:                         IF( J1END.GE.J1 ) THEN
  564:                            DO 160 J1INC = J1, J1END, KD1
  565:                               CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
  566:      $                                   AB( 2, J1INC ), 1, D( J1INC ),
  567:      $                                   WORK( J1INC ) )
  568:   160                      CONTINUE
  569:                         END IF
  570:                         LEND = MIN( KDM1, N-J2 )
  571:                         LAST = J1END + KD1
  572:                         IF( LEND.GT.0 )
  573:      $                     CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
  574:      $                                AB( 2, LAST ), 1, D( LAST ),
  575:      $                                WORK( LAST ) )
  576:                      END IF
  577:                   END IF
  578: *
  579: *
  580: *
  581:                   IF( WANTQ ) THEN
  582: *
  583: *                    accumulate product of plane rotations in Q
  584: *
  585:                      IF( INITQ ) THEN
  586: *
  587: *                 take advantage of the fact that Q was
  588: *                 initially the Identity matrix
  589: *
  590:                         IQEND = MAX( IQEND, J2 )
  591:                         I2 = MAX( 0, K-3 )
  592:                         IQAEND = 1 + I*KD
  593:                         IF( K.EQ.2 )
  594:      $                     IQAEND = IQAEND + KD
  595:                         IQAEND = MIN( IQAEND, IQEND )
  596:                         DO 170 J = J1, J2, KD1
  597:                            IBL = I - I2 / KDM1
  598:                            I2 = I2 + 1
  599:                            IQB = MAX( 1, J-IBL )
  600:                            NQ = 1 + IQAEND - IQB
  601:                            IQAEND = MIN( IQAEND+KD, IQEND )
  602:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
  603:      $                                1, D( J ), WORK( J ) )
  604:   170                   CONTINUE
  605:                      ELSE
  606: *
  607:                         DO 180 J = J1, J2, KD1
  608:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
  609:      $                                D( J ), WORK( J ) )
  610:   180                   CONTINUE
  611:                      END IF
  612:                   END IF
  613: *
  614:                   IF( J2+KDN.GT.N ) THEN
  615: *
  616: *                    adjust J2 to keep within the bounds of the matrix
  617: *
  618:                      NR = NR - 1
  619:                      J2 = J2 - KDN - 1
  620:                   END IF
  621: *
  622:                   DO 190 J = J1, J2, KD1
  623: *
  624: *                    create nonzero element a(j+kd,j-1) outside the
  625: *                    band and store it in WORK
  626: *
  627:                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
  628:                      AB( KD1, J ) = D( J )*AB( KD1, J )
  629:   190             CONTINUE
  630:   200          CONTINUE
  631:   210       CONTINUE
  632:          END IF
  633: *
  634:          IF( KD.GT.0 ) THEN
  635: *
  636: *           make off-diagonal elements real and copy them to E
  637: *
  638:             DO 220 I = 1, N - 1
  639:                T = AB( 2, I )
  640:                ABST = ABS( T )
  641:                AB( 2, I ) = ABST
  642:                E( I ) = ABST
  643:                IF( ABST.NE.ZERO ) THEN
  644:                   T = T / ABST
  645:                ELSE
  646:                   T = CONE
  647:                END IF
  648:                IF( I.LT.N-1 )
  649:      $            AB( 2, I+1 ) = AB( 2, I+1 )*T
  650:                IF( WANTQ ) THEN
  651:                   CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
  652:                END IF
  653:   220       CONTINUE
  654:          ELSE
  655: *
  656: *           set E to zero if original matrix was diagonal
  657: *
  658:             DO 230 I = 1, N - 1
  659:                E( I ) = ZERO
  660:   230       CONTINUE
  661:          END IF
  662: *
  663: *        copy diagonal elements to D
  664: *
  665:          DO 240 I = 1, N
  666:             D( I ) = DBLE( AB( 1, I ) )
  667:   240    CONTINUE
  668:       END IF
  669: *
  670:       RETURN
  671: *
  672: *     End of ZHBTRD
  673: *
  674:       END

CVSweb interface <joel.bertrand@systella.fr>