SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
$ S, LDS, INFO )
*
* -- LAPACK routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
REAL RHO
* ..
* .. Array Arguments ..
REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
$ W( * )
* ..
*
* Purpose
* =======
*
* SLAED9 finds the roots of the secular equation, as defined by the
* values in D, Z, and RHO, between KSTART and KSTOP. It makes the
* appropriate calls to SLAED4 and then stores the new matrix of
* eigenvectors for use in calculating the next level of Z vectors.
*
* Arguments
* =========
*
* K (input) INTEGER
* The number of terms in the rational function to be solved by
* SLAED4. K >= 0.
*
* KSTART (input) INTEGER
* KSTOP (input) INTEGER
* The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
* are to be computed. 1 <= KSTART <= KSTOP <= K.
*
* N (input) INTEGER
* The number of rows and columns in the Q matrix.
* N >= K (delation may result in N > K).
*
* D (output) REAL array, dimension (N)
* D(I) contains the updated eigenvalues
* for KSTART <= I <= KSTOP.
*
* Q (workspace) REAL array, dimension (LDQ,N)
*
* LDQ (input) INTEGER
* The leading dimension of the array Q. LDQ >= max( 1, N ).
*
* RHO (input) REAL
* The value of the parameter in the rank one update equation.
* RHO >= 0 required.
*
* DLAMDA (input) REAL array, dimension (K)
* The first K elements of this array contain the old roots
* of the deflated updating problem. These are the poles
* of the secular equation.
*
* W (input) REAL array, dimension (K)
* The first K elements of this array contain the components
* of the deflation-adjusted updating vector.
*
* S (output) REAL array, dimension (LDS, K)
* Will contain the eigenvectors of the repaired matrix which
* will be stored for subsequent Z vector calculation and
* multiplied by the previously accumulated eigenvectors
* to update the system.
*
* LDS (input) INTEGER
* The leading dimension of S. LDS >= max( 1, K ).
*
* INFO (output) INTEGER
* = 0: successful exit.
* < 0: if INFO = -i, the i-th argument had an illegal value.
* > 0: if INFO = 1, an eigenvalue did not converge
*
* Further Details
* ===============
*
* Based on contributions by
* Jeff Rutter, Computer Science Division, University of California
* at Berkeley, USA
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I, J
REAL TEMP
* ..
* .. External Functions ..
REAL SLAMC3, SNRM2
EXTERNAL SLAMC3, SNRM2
* ..
* .. External Subroutines ..
EXTERNAL SCOPY, SLAED4, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, SIGN, SQRT
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
*
IF( K.LT.0 ) THEN
INFO = -1
ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
INFO = -2
ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
$ THEN
INFO = -3
ELSE IF( N.LT.K ) THEN
INFO = -4
ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
INFO = -7
ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
INFO = -12
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SLAED9', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( K.EQ.0 )
$ RETURN
*
* Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
* be computed with high relative accuracy (barring over/underflow).
* This is a problem on machines without a guard digit in
* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
* The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
* which on any of these machines zeros out the bottommost
* bit of DLAMDA(I) if it is 1; this makes the subsequent
* subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
* occurs. On binary machines with a guard digit (almost all
* machines) it does not change DLAMDA(I) at all. On hexadecimal
* and decimal machines with a guard digit, it slightly
* changes the bottommost bits of DLAMDA(I). It does not account
* for hexadecimal or decimal machines without guard digits
* (we know of none). We use a subroutine call to compute
* 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
* this code.
*
DO 10 I = 1, N
DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
10 CONTINUE
*
DO 20 J = KSTART, KSTOP
CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
*
* If the zero finder fails, the computation is terminated.
*
IF( INFO.NE.0 )
$ GO TO 120
20 CONTINUE
*
IF( K.EQ.1 .OR. K.EQ.2 ) THEN
DO 40 I = 1, K
DO 30 J = 1, K
S( J, I ) = Q( J, I )
30 CONTINUE
40 CONTINUE
GO TO 120
END IF
*
* Compute updated W.
*
CALL SCOPY( K, W, 1, S, 1 )
*
* Initialize W(I) = Q(I,I)
*
CALL SCOPY( K, Q, LDQ+1, W, 1 )
DO 70 J = 1, K
DO 50 I = 1, J - 1
W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
50 CONTINUE
DO 60 I = J + 1, K
W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
60 CONTINUE
70 CONTINUE
DO 80 I = 1, K
W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
80 CONTINUE
*
* Compute eigenvectors of the modified rank-1 modification.
*
DO 110 J = 1, K
DO 90 I = 1, K
Q( I, J ) = W( I ) / Q( I, J )
90 CONTINUE
TEMP = SNRM2( K, Q( 1, J ), 1 )
DO 100 I = 1, K
S( I, J ) = Q( I, J ) / TEMP
100 CONTINUE
110 CONTINUE
*
120 CONTINUE
RETURN
*
* End of SLAED9
*
END