206 lines
6.2 KiB
FortranFixed
206 lines
6.2 KiB
FortranFixed
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
|