265 lines
8.4 KiB
FortranFixed
265 lines
8.4 KiB
FortranFixed
SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX,
|
|
$ CTOT, W, S, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, K, LDQ, N, N1
|
|
REAL RHO
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER CTOT( * ), INDX( * )
|
|
REAL D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ),
|
|
$ S( * ), W( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* SLAED3 finds the roots of the secular equation, as defined by the
|
|
* values in D, W, and RHO, between 1 and K. It makes the
|
|
* appropriate calls to SLAED4 and then updates the eigenvectors by
|
|
* multiplying the matrix of eigenvectors of the pair of eigensystems
|
|
* being combined by the matrix of eigenvectors of the K-by-K system
|
|
* which is solved here.
|
|
*
|
|
* This code makes very mild assumptions about floating point
|
|
* arithmetic. It will work on machines with a guard digit in
|
|
* add/subtract, or on those binary machines without guard digits
|
|
* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
|
|
* It could conceivably fail on hexadecimal or decimal machines
|
|
* without guard digits, but we know of none.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* K (input) INTEGER
|
|
* The number of terms in the rational function to be solved by
|
|
* SLAED4. K >= 0.
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of rows and columns in the Q matrix.
|
|
* N >= K (deflation may result in N>K).
|
|
*
|
|
* N1 (input) INTEGER
|
|
* The location of the last eigenvalue in the leading submatrix.
|
|
* min(1,N) <= N1 <= N/2.
|
|
*
|
|
* D (output) REAL array, dimension (N)
|
|
* D(I) contains the updated eigenvalues for
|
|
* 1 <= I <= K.
|
|
*
|
|
* Q (output) REAL array, dimension (LDQ,N)
|
|
* Initially the first K columns are used as workspace.
|
|
* On output the columns 1 to K contain
|
|
* the updated eigenvectors.
|
|
*
|
|
* 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/output) 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. May be changed on output by
|
|
* having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
|
|
* Cray-2, or Cray C-90, as described above.
|
|
*
|
|
* Q2 (input) REAL array, dimension (LDQ2, N)
|
|
* The first K columns of this matrix contain the non-deflated
|
|
* eigenvectors for the split problem.
|
|
*
|
|
* INDX (input) INTEGER array, dimension (N)
|
|
* The permutation used to arrange the columns of the deflated
|
|
* Q matrix into three groups (see SLAED2).
|
|
* The rows of the eigenvectors found by SLAED4 must be likewise
|
|
* permuted before the matrix multiply can take place.
|
|
*
|
|
* CTOT (input) INTEGER array, dimension (4)
|
|
* A count of the total number of the various types of columns
|
|
* in Q, as described in INDX. The fourth column type is any
|
|
* column which has been deflated.
|
|
*
|
|
* W (input/output) REAL array, dimension (K)
|
|
* The first K elements of this array contain the components
|
|
* of the deflation-adjusted updating vector. Destroyed on
|
|
* output.
|
|
*
|
|
* S (workspace) REAL array, dimension (N1 + 1)*K
|
|
* Will contain the eigenvectors of the repaired matrix which
|
|
* will be 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
|
|
* Modified by Francoise Tisseur, University of Tennessee.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
REAL ONE, ZERO
|
|
PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER I, II, IQ2, J, N12, N2, N23
|
|
REAL TEMP
|
|
* ..
|
|
* .. External Functions ..
|
|
REAL SLAMC3, SNRM2
|
|
EXTERNAL SLAMC3, SNRM2
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL SCOPY, SGEMM, SLACPY, SLAED4, SLASET, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, SIGN, SQRT
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
*
|
|
IF( K.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.K ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
|
|
INFO = -6
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'SLAED3', -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, K
|
|
DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
|
|
10 CONTINUE
|
|
*
|
|
DO 20 J = 1, K
|
|
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 )
|
|
$ GO TO 110
|
|
IF( K.EQ.2 ) THEN
|
|
DO 30 J = 1, K
|
|
W( 1 ) = Q( 1, J )
|
|
W( 2 ) = Q( 2, J )
|
|
II = INDX( 1 )
|
|
Q( 1, J ) = W( II )
|
|
II = INDX( 2 )
|
|
Q( 2, J ) = W( II )
|
|
30 CONTINUE
|
|
GO TO 110
|
|
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 60 J = 1, K
|
|
DO 40 I = 1, J - 1
|
|
W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
|
|
40 CONTINUE
|
|
DO 50 I = J + 1, K
|
|
W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
DO 70 I = 1, K
|
|
W( I ) = SIGN( SQRT( -W( I ) ), S( I ) )
|
|
70 CONTINUE
|
|
*
|
|
* Compute eigenvectors of the modified rank-1 modification.
|
|
*
|
|
DO 100 J = 1, K
|
|
DO 80 I = 1, K
|
|
S( I ) = W( I ) / Q( I, J )
|
|
80 CONTINUE
|
|
TEMP = SNRM2( K, S, 1 )
|
|
DO 90 I = 1, K
|
|
II = INDX( I )
|
|
Q( I, J ) = S( II ) / TEMP
|
|
90 CONTINUE
|
|
100 CONTINUE
|
|
*
|
|
* Compute the updated eigenvectors.
|
|
*
|
|
110 CONTINUE
|
|
*
|
|
N2 = N - N1
|
|
N12 = CTOT( 1 ) + CTOT( 2 )
|
|
N23 = CTOT( 2 ) + CTOT( 3 )
|
|
*
|
|
CALL SLACPY( 'A', N23, K, Q( CTOT( 1 )+1, 1 ), LDQ, S, N23 )
|
|
IQ2 = N1*N12 + 1
|
|
IF( N23.NE.0 ) THEN
|
|
CALL SGEMM( 'N', 'N', N2, K, N23, ONE, Q2( IQ2 ), N2, S, N23,
|
|
$ ZERO, Q( N1+1, 1 ), LDQ )
|
|
ELSE
|
|
CALL SLASET( 'A', N2, K, ZERO, ZERO, Q( N1+1, 1 ), LDQ )
|
|
END IF
|
|
*
|
|
CALL SLACPY( 'A', N12, K, Q, LDQ, S, N12 )
|
|
IF( N12.NE.0 ) THEN
|
|
CALL SGEMM( 'N', 'N', N1, K, N12, ONE, Q2, N1, S, N12, ZERO, Q,
|
|
$ LDQ )
|
|
ELSE
|
|
CALL SLASET( 'A', N1, K, ZERO, ZERO, Q( 1, 1 ), LDQ )
|
|
END IF
|
|
*
|
|
*
|
|
120 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of SLAED3
|
|
*
|
|
END
|