SUBROUTINE SLARRE( 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 REAL TOL * .. * .. Array Arguments .. INTEGER ISPLIT( * ) REAL D( * ), E( * ), GERSCH( * ), W( * ), WOFF( * ), $ WORK( * ) * .. * Common block to return operation count * .. Common blocks .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. REAL ITCNT, OPS * .. * * Purpose * ======= * * Given the tridiagonal matrix T, SLARRE 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 * SSTEGR 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 SLASQ2). As * an added benefit, SLARRE 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) REAL 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) REAL 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) REAL * 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) REAL 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) REAL array, dimension (N) * The NSPLIT base points sigma_i. * * GERSCH (output) REAL array, dimension (2*N) * The n Gerschgorin intervals. * * WORK (input) REAL array, dimension (4*N???) * Workspace. * * INFO (output) INTEGER * Output error code from SLASQ2 * * Further Details * =============== * * Based on contributions by * Inderjit Dhillon, IBM Almaden, USA * Osni Marques, LBNL/NERSC, USA * * ===================================================================== * * .. Parameters .. REAL ZERO, ONE, TWO, FOUR, FOURTH PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, $ FOUR = 4.0E0, FOURTH = ONE / FOUR ) * .. * .. Local Scalars .. INTEGER CNT, I, IBEGIN, IEND, IN, J, JBLK, MAXCNT REAL DELTA, EPS, GL, GU, NRM, OFFD, S, SGNDEF, $ SIGMA, TAU, TMP1, WIDTH * .. * .. External Functions .. REAL SLAMCH EXTERNAL SLAMCH * .. * .. External Subroutines .. EXTERNAL SCOPY, SLASQ2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, REAL * .. * .. Executable Statements .. * INFO = 0 EPS = SLAMCH( '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 + REAL( 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 + REAL( 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 + REAL( 1 ) WORK( I ) = E( I )*E( I ) 30 CONTINUE OPS = OPS + REAL( 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 + REAL( 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 + REAL( 1 ) WORK( 3*IN ) = ONE DELTA = EPS TAU = SGNDEF*NRM 60 CONTINUE OPS = OPS + REAL( 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 + REAL( 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 + REAL( 1 ) DELTA = TWO*DELTA GO TO 60 END IF J = J - 1 80 CONTINUE * OPS = OPS + REAL( 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 SLASQ2( IN, WORK, INFO ) * OPS = OPS + REAL( 2 ) TAU = SGNDEF*WORK( IN ) WORK( 3*IN ) = ONE DELTA = TWO*EPS 100 CONTINUE OPS = OPS + REAL( 2 ) TAU = TAU*( ONE-DELTA ) * OPS = OPS + REAL( 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 + REAL( 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 + REAL( 1 ) DELTA = TWO*DELTA GO TO 100 END IF 120 CONTINUE * SIGMA = SIGMA + TAU CALL SCOPY( IN, WORK, 1, D( IBEGIN ), 1 ) CALL SCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 ) WOFF( JBLK ) = SIGMA * * Update the n Gerschgorin intervals * OPS = OPS + REAL( 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 + REAL( 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 SLASQ2( 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 SLARRE * END