File:  [local] / rpl / lapack / lapack / zhbtrd.f
Revision 1.11: download - view: text, annotated - select for diffs - revision graph
Fri Dec 14 14:22:47 2012 UTC (11 years, 5 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_16, rpl-4_1_15, rpl-4_1_14, rpl-4_1_13, rpl-4_1_12, rpl-4_1_11, HEAD
Mise à jour de lapack.

    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: *> \date November 2011
  151: *
  152: *> \ingroup complex16OTHERcomputational
  153: *
  154: *> \par Further Details:
  155: *  =====================
  156: *>
  157: *> \verbatim
  158: *>
  159: *>  Modified by Linda Kaufman, Bell Labs.
  160: *> \endverbatim
  161: *>
  162: *  =====================================================================
  163:       SUBROUTINE ZHBTRD( VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ,
  164:      $                   WORK, INFO )
  165: *
  166: *  -- LAPACK computational routine (version 3.4.0) --
  167: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  168: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  169: *     November 2011
  170: *
  171: *     .. Scalar Arguments ..
  172:       CHARACTER          UPLO, VECT
  173:       INTEGER            INFO, KD, LDAB, LDQ, N
  174: *     ..
  175: *     .. Array Arguments ..
  176:       DOUBLE PRECISION   D( * ), E( * )
  177:       COMPLEX*16         AB( LDAB, * ), Q( LDQ, * ), WORK( * )
  178: *     ..
  179: *
  180: *  =====================================================================
  181: *
  182: *     .. Parameters ..
  183:       DOUBLE PRECISION   ZERO
  184:       PARAMETER          ( ZERO = 0.0D+0 )
  185:       COMPLEX*16         CZERO, CONE
  186:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
  187:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  188: *     ..
  189: *     .. Local Scalars ..
  190:       LOGICAL            INITQ, UPPER, WANTQ
  191:       INTEGER            I, I2, IBL, INCA, INCX, IQAEND, IQB, IQEND, J,
  192:      $                   J1, J1END, J1INC, J2, JEND, JIN, JINC, K, KD1,
  193:      $                   KDM1, KDN, L, LAST, LEND, NQ, NR, NRT
  194:       DOUBLE PRECISION   ABST
  195:       COMPLEX*16         T, TEMP
  196: *     ..
  197: *     .. External Subroutines ..
  198:       EXTERNAL           XERBLA, ZLACGV, ZLAR2V, ZLARGV, ZLARTG, ZLARTV,
  199:      $                   ZLASET, ZROT, ZSCAL
  200: *     ..
  201: *     .. Intrinsic Functions ..
  202:       INTRINSIC          ABS, DBLE, DCONJG, MAX, MIN
  203: *     ..
  204: *     .. External Functions ..
  205:       LOGICAL            LSAME
  206:       EXTERNAL           LSAME
  207: *     ..
  208: *     .. Executable Statements ..
  209: *
  210: *     Test the input parameters
  211: *
  212:       INITQ = LSAME( VECT, 'V' )
  213:       WANTQ = INITQ .OR. LSAME( VECT, 'U' )
  214:       UPPER = LSAME( UPLO, 'U' )
  215:       KD1 = KD + 1
  216:       KDM1 = KD - 1
  217:       INCX = LDAB - 1
  218:       IQEND = 1
  219: *
  220:       INFO = 0
  221:       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
  222:          INFO = -1
  223:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  224:          INFO = -2
  225:       ELSE IF( N.LT.0 ) THEN
  226:          INFO = -3
  227:       ELSE IF( KD.LT.0 ) THEN
  228:          INFO = -4
  229:       ELSE IF( LDAB.LT.KD1 ) THEN
  230:          INFO = -6
  231:       ELSE IF( LDQ.LT.MAX( 1, N ) .AND. WANTQ ) THEN
  232:          INFO = -10
  233:       END IF
  234:       IF( INFO.NE.0 ) THEN
  235:          CALL XERBLA( 'ZHBTRD', -INFO )
  236:          RETURN
  237:       END IF
  238: *
  239: *     Quick return if possible
  240: *
  241:       IF( N.EQ.0 )
  242:      $   RETURN
  243: *
  244: *     Initialize Q to the unit matrix, if needed
  245: *
  246:       IF( INITQ )
  247:      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
  248: *
  249: *     Wherever possible, plane rotations are generated and applied in
  250: *     vector operations of length NR over the index set J1:J2:KD1.
  251: *
  252: *     The real cosines and complex sines of the plane rotations are
  253: *     stored in the arrays D and WORK.
  254: *
  255:       INCA = KD1*LDAB
  256:       KDN = MIN( N-1, KD )
  257:       IF( UPPER ) THEN
  258: *
  259:          IF( KD.GT.1 ) THEN
  260: *
  261: *           Reduce to complex Hermitian tridiagonal form, working with
  262: *           the upper triangle
  263: *
  264:             NR = 0
  265:             J1 = KDN + 2
  266:             J2 = 1
  267: *
  268:             AB( KD1, 1 ) = DBLE( AB( KD1, 1 ) )
  269:             DO 90 I = 1, N - 2
  270: *
  271: *              Reduce i-th row of matrix to tridiagonal form
  272: *
  273:                DO 80 K = KDN + 1, 2, -1
  274:                   J1 = J1 + KDN
  275:                   J2 = J2 + KDN
  276: *
  277:                   IF( NR.GT.0 ) THEN
  278: *
  279: *                    generate plane rotations to annihilate nonzero
  280: *                    elements which have been created outside the band
  281: *
  282:                      CALL ZLARGV( NR, AB( 1, J1-1 ), INCA, WORK( J1 ),
  283:      $                            KD1, D( J1 ), KD1 )
  284: *
  285: *                    apply rotations from the right
  286: *
  287: *
  288: *                    Dependent on the the number of diagonals either
  289: *                    ZLARTV or ZROT is used
  290: *
  291:                      IF( NR.GE.2*KD-1 ) THEN
  292:                         DO 10 L = 1, KD - 1
  293:                            CALL ZLARTV( NR, AB( L+1, J1-1 ), INCA,
  294:      $                                  AB( L, J1 ), INCA, D( J1 ),
  295:      $                                  WORK( J1 ), KD1 )
  296:    10                   CONTINUE
  297: *
  298:                      ELSE
  299:                         JEND = J1 + ( NR-1 )*KD1
  300:                         DO 20 JINC = J1, JEND, KD1
  301:                            CALL ZROT( KDM1, AB( 2, JINC-1 ), 1,
  302:      $                                AB( 1, JINC ), 1, D( JINC ),
  303:      $                                WORK( JINC ) )
  304:    20                   CONTINUE
  305:                      END IF
  306:                   END IF
  307: *
  308: *
  309:                   IF( K.GT.2 ) THEN
  310:                      IF( K.LE.N-I+1 ) THEN
  311: *
  312: *                       generate plane rotation to annihilate a(i,i+k-1)
  313: *                       within the band
  314: *
  315:                         CALL ZLARTG( AB( KD-K+3, I+K-2 ),
  316:      $                               AB( KD-K+2, I+K-1 ), D( I+K-1 ),
  317:      $                               WORK( I+K-1 ), TEMP )
  318:                         AB( KD-K+3, I+K-2 ) = TEMP
  319: *
  320: *                       apply rotation from the right
  321: *
  322:                         CALL ZROT( K-3, AB( KD-K+4, I+K-2 ), 1,
  323:      $                             AB( KD-K+3, I+K-1 ), 1, D( I+K-1 ),
  324:      $                             WORK( I+K-1 ) )
  325:                      END IF
  326:                      NR = NR + 1
  327:                      J1 = J1 - KDN - 1
  328:                   END IF
  329: *
  330: *                 apply plane rotations from both sides to diagonal
  331: *                 blocks
  332: *
  333:                   IF( NR.GT.0 )
  334:      $               CALL ZLAR2V( NR, AB( KD1, J1-1 ), AB( KD1, J1 ),
  335:      $                            AB( KD, J1 ), INCA, D( J1 ),
  336:      $                            WORK( J1 ), KD1 )
  337: *
  338: *                 apply plane rotations from the left
  339: *
  340:                   IF( NR.GT.0 ) THEN
  341:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
  342:                      IF( 2*KD-1.LT.NR ) THEN
  343: *
  344: *                    Dependent on the the number of diagonals either
  345: *                    ZLARTV or ZROT is used
  346: *
  347:                         DO 30 L = 1, KD - 1
  348:                            IF( J2+L.GT.N ) THEN
  349:                               NRT = NR - 1
  350:                            ELSE
  351:                               NRT = NR
  352:                            END IF
  353:                            IF( NRT.GT.0 )
  354:      $                        CALL ZLARTV( NRT, AB( KD-L, J1+L ), INCA,
  355:      $                                     AB( KD-L+1, J1+L ), INCA,
  356:      $                                     D( J1 ), WORK( J1 ), KD1 )
  357:    30                   CONTINUE
  358:                      ELSE
  359:                         J1END = J1 + KD1*( NR-2 )
  360:                         IF( J1END.GE.J1 ) THEN
  361:                            DO 40 JIN = J1, J1END, KD1
  362:                               CALL ZROT( KD-1, AB( KD-1, JIN+1 ), INCX,
  363:      $                                   AB( KD, JIN+1 ), INCX,
  364:      $                                   D( JIN ), WORK( JIN ) )
  365:    40                      CONTINUE
  366:                         END IF
  367:                         LEND = MIN( KDM1, N-J2 )
  368:                         LAST = J1END + KD1
  369:                         IF( LEND.GT.0 )
  370:      $                     CALL ZROT( LEND, AB( KD-1, LAST+1 ), INCX,
  371:      $                                AB( KD, LAST+1 ), INCX, D( LAST ),
  372:      $                                WORK( LAST ) )
  373:                      END IF
  374:                   END IF
  375: *
  376:                   IF( WANTQ ) THEN
  377: *
  378: *                    accumulate product of plane rotations in Q
  379: *
  380:                      IF( INITQ ) THEN
  381: *
  382: *                 take advantage of the fact that Q was
  383: *                 initially the Identity matrix
  384: *
  385:                         IQEND = MAX( IQEND, J2 )
  386:                         I2 = MAX( 0, K-3 )
  387:                         IQAEND = 1 + I*KD
  388:                         IF( K.EQ.2 )
  389:      $                     IQAEND = IQAEND + KD
  390:                         IQAEND = MIN( IQAEND, IQEND )
  391:                         DO 50 J = J1, J2, KD1
  392:                            IBL = I - I2 / KDM1
  393:                            I2 = I2 + 1
  394:                            IQB = MAX( 1, J-IBL )
  395:                            NQ = 1 + IQAEND - IQB
  396:                            IQAEND = MIN( IQAEND+KD, IQEND )
  397:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
  398:      $                                1, D( J ), DCONJG( WORK( J ) ) )
  399:    50                   CONTINUE
  400:                      ELSE
  401: *
  402:                         DO 60 J = J1, J2, KD1
  403:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
  404:      $                                D( J ), DCONJG( WORK( J ) ) )
  405:    60                   CONTINUE
  406:                      END IF
  407: *
  408:                   END IF
  409: *
  410:                   IF( J2+KDN.GT.N ) THEN
  411: *
  412: *                    adjust J2 to keep within the bounds of the matrix
  413: *
  414:                      NR = NR - 1
  415:                      J2 = J2 - KDN - 1
  416:                   END IF
  417: *
  418:                   DO 70 J = J1, J2, KD1
  419: *
  420: *                    create nonzero element a(j-1,j+kd) outside the band
  421: *                    and store it in WORK
  422: *
  423:                      WORK( J+KD ) = WORK( J )*AB( 1, J+KD )
  424:                      AB( 1, J+KD ) = D( J )*AB( 1, J+KD )
  425:    70             CONTINUE
  426:    80          CONTINUE
  427:    90       CONTINUE
  428:          END IF
  429: *
  430:          IF( KD.GT.0 ) THEN
  431: *
  432: *           make off-diagonal elements real and copy them to E
  433: *
  434:             DO 100 I = 1, N - 1
  435:                T = AB( KD, I+1 )
  436:                ABST = ABS( T )
  437:                AB( KD, I+1 ) = ABST
  438:                E( I ) = ABST
  439:                IF( ABST.NE.ZERO ) THEN
  440:                   T = T / ABST
  441:                ELSE
  442:                   T = CONE
  443:                END IF
  444:                IF( I.LT.N-1 )
  445:      $            AB( KD, I+2 ) = AB( KD, I+2 )*T
  446:                IF( WANTQ ) THEN
  447:                   CALL ZSCAL( N, DCONJG( T ), Q( 1, I+1 ), 1 )
  448:                END IF
  449:   100       CONTINUE
  450:          ELSE
  451: *
  452: *           set E to zero if original matrix was diagonal
  453: *
  454:             DO 110 I = 1, N - 1
  455:                E( I ) = ZERO
  456:   110       CONTINUE
  457:          END IF
  458: *
  459: *        copy diagonal elements to D
  460: *
  461:          DO 120 I = 1, N
  462:             D( I ) = AB( KD1, I )
  463:   120    CONTINUE
  464: *
  465:       ELSE
  466: *
  467:          IF( KD.GT.1 ) THEN
  468: *
  469: *           Reduce to complex Hermitian tridiagonal form, working with
  470: *           the lower triangle
  471: *
  472:             NR = 0
  473:             J1 = KDN + 2
  474:             J2 = 1
  475: *
  476:             AB( 1, 1 ) = DBLE( AB( 1, 1 ) )
  477:             DO 210 I = 1, N - 2
  478: *
  479: *              Reduce i-th column of matrix to tridiagonal form
  480: *
  481:                DO 200 K = KDN + 1, 2, -1
  482:                   J1 = J1 + KDN
  483:                   J2 = J2 + KDN
  484: *
  485:                   IF( NR.GT.0 ) THEN
  486: *
  487: *                    generate plane rotations to annihilate nonzero
  488: *                    elements which have been created outside the band
  489: *
  490:                      CALL ZLARGV( NR, AB( KD1, J1-KD1 ), INCA,
  491:      $                            WORK( J1 ), KD1, D( J1 ), KD1 )
  492: *
  493: *                    apply plane rotations from one side
  494: *
  495: *
  496: *                    Dependent on the the number of diagonals either
  497: *                    ZLARTV or ZROT is used
  498: *
  499:                      IF( NR.GT.2*KD-1 ) THEN
  500:                         DO 130 L = 1, KD - 1
  501:                            CALL ZLARTV( NR, AB( KD1-L, J1-KD1+L ), INCA,
  502:      $                                  AB( KD1-L+1, J1-KD1+L ), INCA,
  503:      $                                  D( J1 ), WORK( J1 ), KD1 )
  504:   130                   CONTINUE
  505:                      ELSE
  506:                         JEND = J1 + KD1*( NR-1 )
  507:                         DO 140 JINC = J1, JEND, KD1
  508:                            CALL ZROT( KDM1, AB( KD, JINC-KD ), INCX,
  509:      $                                AB( KD1, JINC-KD ), INCX,
  510:      $                                D( JINC ), WORK( JINC ) )
  511:   140                   CONTINUE
  512:                      END IF
  513: *
  514:                   END IF
  515: *
  516:                   IF( K.GT.2 ) THEN
  517:                      IF( K.LE.N-I+1 ) THEN
  518: *
  519: *                       generate plane rotation to annihilate a(i+k-1,i)
  520: *                       within the band
  521: *
  522:                         CALL ZLARTG( AB( K-1, I ), AB( K, I ),
  523:      $                               D( I+K-1 ), WORK( I+K-1 ), TEMP )
  524:                         AB( K-1, I ) = TEMP
  525: *
  526: *                       apply rotation from the left
  527: *
  528:                         CALL ZROT( K-3, AB( K-2, I+1 ), LDAB-1,
  529:      $                             AB( K-1, I+1 ), LDAB-1, D( I+K-1 ),
  530:      $                             WORK( I+K-1 ) )
  531:                      END IF
  532:                      NR = NR + 1
  533:                      J1 = J1 - KDN - 1
  534:                   END IF
  535: *
  536: *                 apply plane rotations from both sides to diagonal
  537: *                 blocks
  538: *
  539:                   IF( NR.GT.0 )
  540:      $               CALL ZLAR2V( NR, AB( 1, J1-1 ), AB( 1, J1 ),
  541:      $                            AB( 2, J1-1 ), INCA, D( J1 ),
  542:      $                            WORK( J1 ), KD1 )
  543: *
  544: *                 apply plane rotations from the right
  545: *
  546: *
  547: *                    Dependent on the the number of diagonals either
  548: *                    ZLARTV or ZROT is used
  549: *
  550:                   IF( NR.GT.0 ) THEN
  551:                      CALL ZLACGV( NR, WORK( J1 ), KD1 )
  552:                      IF( NR.GT.2*KD-1 ) THEN
  553:                         DO 150 L = 1, KD - 1
  554:                            IF( J2+L.GT.N ) THEN
  555:                               NRT = NR - 1
  556:                            ELSE
  557:                               NRT = NR
  558:                            END IF
  559:                            IF( NRT.GT.0 )
  560:      $                        CALL ZLARTV( NRT, AB( L+2, J1-1 ), INCA,
  561:      $                                     AB( L+1, J1 ), INCA, D( J1 ),
  562:      $                                     WORK( J1 ), KD1 )
  563:   150                   CONTINUE
  564:                      ELSE
  565:                         J1END = J1 + KD1*( NR-2 )
  566:                         IF( J1END.GE.J1 ) THEN
  567:                            DO 160 J1INC = J1, J1END, KD1
  568:                               CALL ZROT( KDM1, AB( 3, J1INC-1 ), 1,
  569:      $                                   AB( 2, J1INC ), 1, D( J1INC ),
  570:      $                                   WORK( J1INC ) )
  571:   160                      CONTINUE
  572:                         END IF
  573:                         LEND = MIN( KDM1, N-J2 )
  574:                         LAST = J1END + KD1
  575:                         IF( LEND.GT.0 )
  576:      $                     CALL ZROT( LEND, AB( 3, LAST-1 ), 1,
  577:      $                                AB( 2, LAST ), 1, D( LAST ),
  578:      $                                WORK( LAST ) )
  579:                      END IF
  580:                   END IF
  581: *
  582: *
  583: *
  584:                   IF( WANTQ ) THEN
  585: *
  586: *                    accumulate product of plane rotations in Q
  587: *
  588:                      IF( INITQ ) THEN
  589: *
  590: *                 take advantage of the fact that Q was
  591: *                 initially the Identity matrix
  592: *
  593:                         IQEND = MAX( IQEND, J2 )
  594:                         I2 = MAX( 0, K-3 )
  595:                         IQAEND = 1 + I*KD
  596:                         IF( K.EQ.2 )
  597:      $                     IQAEND = IQAEND + KD
  598:                         IQAEND = MIN( IQAEND, IQEND )
  599:                         DO 170 J = J1, J2, KD1
  600:                            IBL = I - I2 / KDM1
  601:                            I2 = I2 + 1
  602:                            IQB = MAX( 1, J-IBL )
  603:                            NQ = 1 + IQAEND - IQB
  604:                            IQAEND = MIN( IQAEND+KD, IQEND )
  605:                            CALL ZROT( NQ, Q( IQB, J-1 ), 1, Q( IQB, J ),
  606:      $                                1, D( J ), WORK( J ) )
  607:   170                   CONTINUE
  608:                      ELSE
  609: *
  610:                         DO 180 J = J1, J2, KD1
  611:                            CALL ZROT( N, Q( 1, J-1 ), 1, Q( 1, J ), 1,
  612:      $                                D( J ), WORK( J ) )
  613:   180                   CONTINUE
  614:                      END IF
  615:                   END IF
  616: *
  617:                   IF( J2+KDN.GT.N ) THEN
  618: *
  619: *                    adjust J2 to keep within the bounds of the matrix
  620: *
  621:                      NR = NR - 1
  622:                      J2 = J2 - KDN - 1
  623:                   END IF
  624: *
  625:                   DO 190 J = J1, J2, KD1
  626: *
  627: *                    create nonzero element a(j+kd,j-1) outside the
  628: *                    band and store it in WORK
  629: *
  630:                      WORK( J+KD ) = WORK( J )*AB( KD1, J )
  631:                      AB( KD1, J ) = D( J )*AB( KD1, J )
  632:   190             CONTINUE
  633:   200          CONTINUE
  634:   210       CONTINUE
  635:          END IF
  636: *
  637:          IF( KD.GT.0 ) THEN
  638: *
  639: *           make off-diagonal elements real and copy them to E
  640: *
  641:             DO 220 I = 1, N - 1
  642:                T = AB( 2, I )
  643:                ABST = ABS( T )
  644:                AB( 2, I ) = ABST
  645:                E( I ) = ABST
  646:                IF( ABST.NE.ZERO ) THEN
  647:                   T = T / ABST
  648:                ELSE
  649:                   T = CONE
  650:                END IF
  651:                IF( I.LT.N-1 )
  652:      $            AB( 2, I+1 ) = AB( 2, I+1 )*T
  653:                IF( WANTQ ) THEN
  654:                   CALL ZSCAL( N, T, Q( 1, I+1 ), 1 )
  655:                END IF
  656:   220       CONTINUE
  657:          ELSE
  658: *
  659: *           set E to zero if original matrix was diagonal
  660: *
  661:             DO 230 I = 1, N - 1
  662:                E( I ) = ZERO
  663:   230       CONTINUE
  664:          END IF
  665: *
  666: *        copy diagonal elements to D
  667: *
  668:          DO 240 I = 1, N
  669:             D( I ) = AB( 1, I )
  670:   240    CONTINUE
  671:       END IF
  672: *
  673:       RETURN
  674: *
  675: *     End of ZHBTRD
  676: *
  677:       END

CVSweb interface <joel.bertrand@systella.fr>