dlarre.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

      SUBROUTINE DLARRE( N, D, E, TOL, NSPLIT, ISPLIT, M, W, WOFF,
     $                   GERSCH, WORK, INFO )
*
*  -- LAPACK auxiliary routine (instru to count ops, version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      INTEGER            INFO, M, N, NSPLIT
      DOUBLE PRECISION   TOL
*     ..
*     .. Array Arguments ..
      INTEGER            ISPLIT( * )
      DOUBLE PRECISION   D( * ), E( * ), GERSCH( * ), W( * ), WOFF( * ),
     $                   WORK( * )
*     ..
*     Common block to return operation count
*     .. Common blocks ..
      COMMON             / LATIME / OPS, ITCNT
*     ..
*     .. Scalars in Common ..
      DOUBLE PRECISION   ITCNT, OPS
*     ..
*
*  Purpose
*  =======
*
*  Given the tridiagonal matrix T, DLARRE sets "small" off-diagonal
*  elements to zero, and for each unreduced block T_i, it finds
*  (i) the numbers sigma_i
*  (ii) the base T_i - sigma_i I = L_i D_i L_i^T representations and
*  (iii) eigenvalues of each L_i D_i L_i^T.
*  The representations and eigenvalues found are then used by
*  DSTEGR to compute the eigenvectors of a symmetric tridiagonal
*  matrix. Currently, the base representations are limited to being
*  positive or negative definite, and the eigenvalues of the definite
*  matrices are found by the dqds algorithm (subroutine DLASQ2). As
*  an added benefit, DLARRE also outputs the n Gerschgorin
*  intervals for each L_i D_i L_i^T.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix.
*
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the n diagonal elements of the tridiagonal
*          matrix T.
*          On exit, the n diagonal elements of the diagonal
*          matrices D_i.
*
*  E       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, the (n-1) subdiagonal elements of the tridiagonal
*          matrix T; E(N) need not be set.
*          On exit, the subdiagonal elements of the unit bidiagonal
*          matrices L_i.
*
*  TOL     (input) DOUBLE PRECISION
*          The threshold for splitting. If on input |E(i)| < TOL, then
*          the matrix T is split into smaller blocks.
*
*  NSPLIT  (input) INTEGER
*          The number of blocks T splits into. 1 <= NSPLIT <= N.
*
*  ISPLIT  (output) INTEGER array, dimension (2*N)
*          The splitting points, at which T breaks up into submatrices.
*          The first submatrix consists of rows/columns 1 to ISPLIT(1),
*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
*          etc., and the NSPLIT-th consists of rows/columns
*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
*
*  M       (output) INTEGER
*          The total number of eigenvalues (of all the L_i D_i L_i^T)
*          found.
*
*  W       (output) DOUBLE PRECISION array, dimension (N)
*          The first M elements contain the eigenvalues. The
*          eigenvalues of each of the blocks, L_i D_i L_i^T, are
*          sorted in ascending order.
*
*  WOFF    (output) DOUBLE PRECISION array, dimension (N)
*          The NSPLIT base points sigma_i.
*
*  GERSCH  (output) DOUBLE PRECISION array, dimension (2*N)
*          The n Gerschgorin intervals.
*
*  WORK    (input) DOUBLE PRECISION array, dimension (4*N???)
*          Workspace.
*
*  INFO    (output) INTEGER
*          Output error code from DLASQ2
*
*  Further Details
*  ===============
*
*  Based on contributions by
*     Inderjit Dhillon, IBM Almaden, USA
*     Osni Marques, LBNL/NERSC, USA
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TWO, FOUR, FOURTH
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                   FOUR = 4.0D0, FOURTH = ONE / FOUR )
*     ..
*     .. Local Scalars ..
      INTEGER            CNT, I, IBEGIN, IEND, IN, J, JBLK, MAXCNT
      DOUBLE PRECISION   DELTA, EPS, GL, GU, NRM, OFFD, S, SGNDEF,
     $                   SIGMA, TAU, TMP1, WIDTH
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           DCOPY, DLASQ2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX, MIN
*     ..
*     .. Executable Statements ..
*
      INFO = 0
      EPS = DLAMCH( 'Precision' )
*
*     Compute Splitting Points
*
      NSPLIT = 1
      DO 10 I = 1, N - 1
         IF( ABS( E( I ) ).LE.TOL ) THEN
            ISPLIT( NSPLIT ) = I
            NSPLIT = NSPLIT + 1
         END IF
   10 CONTINUE
      ISPLIT( NSPLIT ) = N
*
      IBEGIN = 1
      DO 170 JBLK = 1, NSPLIT
         IEND = ISPLIT( JBLK )
         IF( IBEGIN.EQ.IEND ) THEN
            W( IBEGIN ) = D( IBEGIN )
            WOFF( JBLK ) = ZERO
            IBEGIN = IEND + 1
            GO TO 170
         END IF
         IN = IEND - IBEGIN + 1
*
*        Form the n Gerschgorin intervals
*
         OPS = OPS + DBLE( 4 )
         GL = D( IBEGIN ) - ABS( E( IBEGIN ) )
         GU = D( IBEGIN ) + ABS( E( IBEGIN ) )
         GERSCH( 2*IBEGIN-1 ) = GL
         GERSCH( 2*IBEGIN ) = GU
         GERSCH( 2*IEND-1 ) = D( IEND ) - ABS( E( IEND-1 ) )
         GERSCH( 2*IEND ) = D( IEND ) + ABS( E( IEND-1 ) )
         GL = MIN( GERSCH( 2*IEND-1 ), GL )
         GU = MAX( GERSCH( 2*IEND ), GU )
         DO 20 I = IBEGIN + 1, IEND - 1
            OPS = OPS + DBLE( 3 )
            OFFD = ABS( E( I-1 ) ) + ABS( E( I ) )
            GERSCH( 2*I-1 ) = D( I ) - OFFD
            GL = MIN( GERSCH( 2*I-1 ), GL )
            GERSCH( 2*I ) = D( I ) + OFFD
            GU = MAX( GERSCH( 2*I ), GU )
   20    CONTINUE
         NRM = MAX( ABS( GL ), ABS( GU ) )
*
*        Find the number SIGMA where the base representation
*        T - sigma I = L D L^T is to be formed.
*
         WIDTH = GU - GL
         DO 30 I = IBEGIN, IEND - 1
            OPS = OPS + DBLE( 1 )
            WORK( I ) = E( I )*E( I )
   30    CONTINUE
         OPS = OPS + DBLE( 6 )
         DO 50 J = 1, 2
            IF( J.EQ.1 ) THEN
               TAU = GL + FOURTH*WIDTH
            ELSE
               TAU = GU - FOURTH*WIDTH
            END IF
            TMP1 = D( IBEGIN ) - TAU
            IF( TMP1.LT.ZERO ) THEN
               CNT = 1
            ELSE
               CNT = 0
            END IF
            DO 40 I = IBEGIN + 1, IEND
               OPS = OPS + DBLE( 3 )
               TMP1 = D( I ) - TAU - WORK( I-1 ) / TMP1
               IF( TMP1.LT.ZERO )
     $            CNT = CNT + 1
   40       CONTINUE
            IF( CNT.EQ.0 ) THEN
               GL = TAU
            ELSE IF( CNT.EQ.IN ) THEN
               GU = TAU
            END IF
            IF( J.EQ.1 ) THEN
               MAXCNT = CNT
               SIGMA = GL
               SGNDEF = ONE
            ELSE
               IF( IN-CNT.GT.MAXCNT ) THEN
                  SIGMA = GU
                  SGNDEF = -ONE
               END IF
            END IF
   50    CONTINUE
*
*        Find the base L D L^T representation
*
         OPS = OPS + DBLE( 1 )
         WORK( 3*IN ) = ONE
         DELTA = EPS
         TAU = SGNDEF*NRM
   60    CONTINUE
         OPS = OPS + DBLE( 3+5*( IN-1 ) )
         SIGMA = SIGMA - DELTA*TAU
         WORK( 1 ) = D( IBEGIN ) - SIGMA
         J = IBEGIN
         DO 70 I = 1, IN - 1
            WORK( 2*IN+I ) = ONE / WORK( 2*I-1 )
            TMP1 = E( J )*WORK( 2*IN+I )
            WORK( 2*I+1 ) = ( D( J+1 )-SIGMA ) - TMP1*E( J )
            WORK( 2*I ) = TMP1
            J = J + 1
   70    CONTINUE
         OPS = OPS + DBLE( IN )
         DO 80 I = IN, 1, -1
            TMP1 = SGNDEF*WORK( 2*I-1 )
            IF( TMP1.LT.ZERO .OR. WORK( 2*IN+I ).EQ.ZERO .OR. .NOT.
     $          ( TMP1.GT.ZERO .OR. TMP1.LT.ONE ) ) THEN
               OPS = OPS + DBLE( 1 )
               DELTA = TWO*DELTA
               GO TO 60
            END IF
            J = J - 1
   80    CONTINUE
*
         OPS = OPS + DBLE( IN-1 )
         J = IBEGIN
         D( IBEGIN ) = WORK( 1 )
         WORK( 1 ) = ABS( WORK( 1 ) )
         DO 90 I = 1, IN - 1
            TMP1 = E( J )
            E( J ) = WORK( 2*I )
            WORK( 2*I ) = ABS( TMP1*WORK( 2*I ) )
            J = J + 1
            D( J ) = WORK( 2*I+1 )
            WORK( 2*I+1 ) = ABS( WORK( 2*I+1 ) )
   90    CONTINUE
*
         CALL DLASQ2( IN, WORK, INFO )
*
         OPS = OPS + DBLE( 2 )
         TAU = SGNDEF*WORK( IN )
         WORK( 3*IN ) = ONE
         DELTA = TWO*EPS
  100    CONTINUE
         OPS = OPS + DBLE( 2 )
         TAU = TAU*( ONE-DELTA )
*
         OPS = OPS + DBLE( 9*( IN-1 )+1 )
         S = -TAU
         J = IBEGIN
         DO 110 I = 1, IN - 1
            WORK( I ) = D( J ) + S
            WORK( 2*IN+I ) = ONE / WORK( I )
*           WORK( N+I ) = ( E( I ) * D( I ) ) / WORK( I )
            WORK( IN+I ) = ( E( J )*D( J ) )*WORK( 2*IN+I )
            S = S*WORK( IN+I )*E( J ) - TAU
            J = J + 1
  110    CONTINUE
         WORK( IN ) = D( IEND ) + S
*
*        Checking to see if all the diagonal elements of the new
*        L D L^T representation have the same sign
*
         OPS = OPS + DBLE( IN+1 )
         DO 120 I = IN, 1, -1
            TMP1 = SGNDEF*WORK( I )
            IF( TMP1.LT.ZERO .OR. WORK( 2*IN+I ).EQ.ZERO .OR. .NOT.
     $          ( TMP1.GT.ZERO .OR. TMP1.LT.ONE ) ) THEN
               OPS = OPS + DBLE( 1 )
               DELTA = TWO*DELTA
               GO TO 100
            END IF
  120    CONTINUE
*
         SIGMA = SIGMA + TAU
         CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 )
         CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 )
         WOFF( JBLK ) = SIGMA
*
*        Update the n Gerschgorin intervals
*
         OPS = OPS + DBLE( 2 )
         DO 130 I = IBEGIN, IEND
            GERSCH( 2*I-1 ) = GERSCH( 2*I-1 ) - SIGMA
            GERSCH( 2*I ) = GERSCH( 2*I ) - SIGMA
  130    CONTINUE
*
*        Compute the eigenvalues of L D L^T.
*
         J = IBEGIN
         OPS = OPS + DBLE( 2*( IN-1 ) )
         DO 140 I = 1, IN - 1
            WORK( 2*I-1 ) = ABS( D( J ) )
            WORK( 2*I ) = E( J )*E( J )*WORK( 2*I-1 )
            J = J + 1
  140    CONTINUE
         WORK( 2*IN-1 ) = ABS( D( IEND ) )
*
         CALL DLASQ2( IN, WORK, INFO )
*
         J = IBEGIN
         IF( SGNDEF.GT.ZERO ) THEN
            DO 150 I = 1, IN
               W( J ) = WORK( IN-I+1 )
               J = J + 1
  150       CONTINUE
         ELSE
            DO 160 I = 1, IN
               W( J ) = -WORK( I )
               J = J + 1
  160       CONTINUE
         END IF
         IBEGIN = IEND + 1
  170 CONTINUE
      M = N
*
      RETURN
*
*     End of DLARRE
*
      END