d7/d95/slaed9_8f_source.html

00001       SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
00002      $                   S, LDS, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
00011       REAL               RHO
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
00015      $                   W( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SLAED9 finds the roots of the secular equation, as defined by the
00022 *  values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
00023 *  appropriate calls to SLAED4 and then stores the new matrix of
00024 *  eigenvectors for use in calculating the next level of Z vectors.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  K       (input) INTEGER
00030 *          The number of terms in the rational function to be solved by
00031 *          SLAED4.  K >= 0.
00032 *
00033 *  KSTART  (input) INTEGER
00034 *  KSTOP   (input) INTEGER
00035 *          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
00036 *          are to be computed.  1 <= KSTART <= KSTOP <= K.
00037 *
00038 *  N       (input) INTEGER
00039 *          The number of rows and columns in the Q matrix.
00040 *          N >= K (delation may result in N > K).
00041 *
00042 *  D       (output) REAL array, dimension (N)
00043 *          D(I) contains the updated eigenvalues
00044 *          for KSTART <= I <= KSTOP.
00045 *
00046 *  Q       (workspace) REAL array, dimension (LDQ,N)
00047 *
00048 *  LDQ     (input) INTEGER
00049 *          The leading dimension of the array Q.  LDQ >= max( 1, N ).
00050 *
00051 *  RHO     (input) REAL
00052 *          The value of the parameter in the rank one update equation.
00053 *          RHO >= 0 required.
00054 *
00055 *  DLAMDA  (input) REAL array, dimension (K)
00056 *          The first K elements of this array contain the old roots
00057 *          of the deflated updating problem.  These are the poles
00058 *          of the secular equation.
00059 *
00060 *  W       (input) REAL array, dimension (K)
00061 *          The first K elements of this array contain the components
00062 *          of the deflation-adjusted updating vector.
00063 *
00064 *  S       (output) REAL array, dimension (LDS, K)
00065 *          Will contain the eigenvectors of the repaired matrix which
00066 *          will be stored for subsequent Z vector calculation and
00067 *          multiplied by the previously accumulated eigenvectors
00068 *          to update the system.
00069 *
00070 *  LDS     (input) INTEGER
00071 *          The leading dimension of S.  LDS >= max( 1, K ).
00072 *
00073 *  INFO    (output) INTEGER
00074 *          = 0:  successful exit.
00075 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00076 *          > 0:  if INFO = 1, an eigenvalue did not converge
00077 *
00078 *  Further Details
00079 *  ===============
00080 *
00081 *  Based on contributions by
00082 *     Jeff Rutter, Computer Science Division, University of California
00083 *     at Berkeley, USA
00084 *
00085 *  =====================================================================
00086 *
00087 *     .. Local Scalars ..
00088       INTEGER            I, J
00089       REAL               TEMP
00090 *     ..
00091 *     .. External Functions ..
00092       REAL               SLAMC3, SNRM2
00093       EXTERNAL           SLAMC3, SNRM2
00094 *     ..
00095 *     .. External Subroutines ..
00096       EXTERNAL           SCOPY, SLAED4, XERBLA
00097 *     ..
00098 *     .. Intrinsic Functions ..
00099       INTRINSIC          MAX, SIGN, SQRT
00100 *     ..
00101 *     .. Executable Statements ..
00102 *
00103 *     Test the input parameters.
00104 *
00105       INFO = 0
00106 *
00107       IF( K.LT.0 ) THEN
00108          INFO = -1
00109       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
00110          INFO = -2
00111       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
00112      $          THEN
00113          INFO = -3
00114       ELSE IF( N.LT.K ) THEN
00115          INFO = -4
00116       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
00117          INFO = -7
00118       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
00119          INFO = -12
00120       END IF
00121       IF( INFO.NE.0 ) THEN
00122          CALL XERBLA( 'SLAED9', -INFO )
00123          RETURN
00124       END IF
00125 *
00126 *     Quick return if possible
00127 *
00128       IF( K.EQ.0 )
00129      $   RETURN
00130 *
00131 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
00132 *     be computed with high relative accuracy (barring over/underflow).
00133 *     This is a problem on machines without a guard digit in
00134 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00135 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
00136 *     which on any of these machines zeros out the bottommost
00137 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
00138 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
00139 *     occurs. On binary machines with a guard digit (almost all
00140 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
00141 *     and decimal machines with a guard digit, it slightly
00142 *     changes the bottommost bits of DLAMDA(I). It does not account
00143 *     for hexadecimal or decimal machines without guard digits
00144 *     (we know of none). We use a subroutine call to compute
00145 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00146 *     this code.
00147 *
00148       DO 10 I = 1, N
00149          DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
00150    10 CONTINUE
00151 *
00152       DO 20 J = KSTART, KSTOP
00153          CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
00154 *
00155 *        If the zero finder fails, the computation is terminated.
00156 *
00157          IF( INFO.NE.0 )
00158      $      GO TO 120
00159    20 CONTINUE
00160 *
00161       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
00162          DO 40 I = 1, K
00163             DO 30 J = 1, K
00164                S( J, I ) = Q( J, I )
00165    30       CONTINUE
00166    40    CONTINUE
00167          GO TO 120
00168       END IF
00169 *
00170 *     Compute updated W.
00171 *
00172       CALL SCOPY( K, W, 1, S, 1 )
00173 *
00174 *     Initialize W(I) = Q(I,I)
00175 *
00176       CALL SCOPY( K, Q, LDQ+1, W, 1 )
00177       DO 70 J = 1, K
00178          DO 50 I = 1, J - 1
00179             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00180    50    CONTINUE
00181          DO 60 I = J + 1, K
00182             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00183    60    CONTINUE
00184    70 CONTINUE
00185       DO 80 I = 1, K
00186          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
00187    80 CONTINUE
00188 *
00189 *     Compute eigenvectors of the modified rank-1 modification.
00190 *
00191       DO 110 J = 1, K
00192          DO 90 I = 1, K
00193             Q( I, J ) = W( I ) / Q( I, J )
00194    90    CONTINUE
00195          TEMP = SNRM2( K, Q( 1, J ), 1 )
00196          DO 100 I = 1, K
00197             S( I, J ) = Q( I, J ) / TEMP
00198   100    CONTINUE
00199   110 CONTINUE
00200 *
00201   120 CONTINUE
00202       RETURN
00203 *
00204 *     End of SLAED9
00205 *
00206       END