LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaed9.f
Go to the documentation of this file.
00001 *> \brief \b DLAED9
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAED9 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed9.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed9.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed9.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
00022 *                          S, LDS, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
00026 *       DOUBLE PRECISION   RHO
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
00030 *      $                   W( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLAED9 finds the roots of the secular equation, as defined by the
00040 *> values in D, Z, and RHO, between KSTART and KSTOP.  It makes the
00041 *> appropriate calls to DLAED4 and then stores the new matrix of
00042 *> eigenvectors for use in calculating the next level of Z vectors.
00043 *> \endverbatim
00044 *
00045 *  Arguments:
00046 *  ==========
00047 *
00048 *> \param[in] K
00049 *> \verbatim
00050 *>          K is INTEGER
00051 *>          The number of terms in the rational function to be solved by
00052 *>          DLAED4.  K >= 0.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] KSTART
00056 *> \verbatim
00057 *>          KSTART is INTEGER
00058 *> \endverbatim
00059 *>
00060 *> \param[in] KSTOP
00061 *> \verbatim
00062 *>          KSTOP is INTEGER
00063 *>          The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
00064 *>          are to be computed.  1 <= KSTART <= KSTOP <= K.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] N
00068 *> \verbatim
00069 *>          N is INTEGER
00070 *>          The number of rows and columns in the Q matrix.
00071 *>          N >= K (delation may result in N > K).
00072 *> \endverbatim
00073 *>
00074 *> \param[out] D
00075 *> \verbatim
00076 *>          D is DOUBLE PRECISION array, dimension (N)
00077 *>          D(I) contains the updated eigenvalues
00078 *>          for KSTART <= I <= KSTOP.
00079 *> \endverbatim
00080 *>
00081 *> \param[out] Q
00082 *> \verbatim
00083 *>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
00084 *> \endverbatim
00085 *>
00086 *> \param[in] LDQ
00087 *> \verbatim
00088 *>          LDQ is INTEGER
00089 *>          The leading dimension of the array Q.  LDQ >= max( 1, N ).
00090 *> \endverbatim
00091 *>
00092 *> \param[in] RHO
00093 *> \verbatim
00094 *>          RHO is DOUBLE PRECISION
00095 *>          The value of the parameter in the rank one update equation.
00096 *>          RHO >= 0 required.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] DLAMDA
00100 *> \verbatim
00101 *>          DLAMDA is DOUBLE PRECISION array, dimension (K)
00102 *>          The first K elements of this array contain the old roots
00103 *>          of the deflated updating problem.  These are the poles
00104 *>          of the secular equation.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] W
00108 *> \verbatim
00109 *>          W is DOUBLE PRECISION array, dimension (K)
00110 *>          The first K elements of this array contain the components
00111 *>          of the deflation-adjusted updating vector.
00112 *> \endverbatim
00113 *>
00114 *> \param[out] S
00115 *> \verbatim
00116 *>          S is DOUBLE PRECISION array, dimension (LDS, K)
00117 *>          Will contain the eigenvectors of the repaired matrix which
00118 *>          will be stored for subsequent Z vector calculation and
00119 *>          multiplied by the previously accumulated eigenvectors
00120 *>          to update the system.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDS
00124 *> \verbatim
00125 *>          LDS is INTEGER
00126 *>          The leading dimension of S.  LDS >= max( 1, K ).
00127 *> \endverbatim
00128 *>
00129 *> \param[out] INFO
00130 *> \verbatim
00131 *>          INFO is INTEGER
00132 *>          = 0:  successful exit.
00133 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00134 *>          > 0:  if INFO = 1, an eigenvalue did not converge
00135 *> \endverbatim
00136 *
00137 *  Authors:
00138 *  ========
00139 *
00140 *> \author Univ. of Tennessee 
00141 *> \author Univ. of California Berkeley 
00142 *> \author Univ. of Colorado Denver 
00143 *> \author NAG Ltd. 
00144 *
00145 *> \date November 2011
00146 *
00147 *> \ingroup auxOTHERcomputational
00148 *
00149 *> \par Contributors:
00150 *  ==================
00151 *>
00152 *> Jeff Rutter, Computer Science Division, University of California
00153 *> at Berkeley, USA
00154 *
00155 *  =====================================================================
00156       SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
00157      $                   S, LDS, INFO )
00158 *
00159 *  -- LAPACK computational routine (version 3.4.0) --
00160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00162 *     November 2011
00163 *
00164 *     .. Scalar Arguments ..
00165       INTEGER            INFO, K, KSTART, KSTOP, LDQ, LDS, N
00166       DOUBLE PRECISION   RHO
00167 *     ..
00168 *     .. Array Arguments ..
00169       DOUBLE PRECISION   D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
00170      $                   W( * )
00171 *     ..
00172 *
00173 *  =====================================================================
00174 *
00175 *     .. Local Scalars ..
00176       INTEGER            I, J
00177       DOUBLE PRECISION   TEMP
00178 *     ..
00179 *     .. External Functions ..
00180       DOUBLE PRECISION   DLAMC3, DNRM2
00181       EXTERNAL           DLAMC3, DNRM2
00182 *     ..
00183 *     .. External Subroutines ..
00184       EXTERNAL           DCOPY, DLAED4, XERBLA
00185 *     ..
00186 *     .. Intrinsic Functions ..
00187       INTRINSIC          MAX, SIGN, SQRT
00188 *     ..
00189 *     .. Executable Statements ..
00190 *
00191 *     Test the input parameters.
00192 *
00193       INFO = 0
00194 *
00195       IF( K.LT.0 ) THEN
00196          INFO = -1
00197       ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
00198          INFO = -2
00199       ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
00200      $          THEN
00201          INFO = -3
00202       ELSE IF( N.LT.K ) THEN
00203          INFO = -4
00204       ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
00205          INFO = -7
00206       ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
00207          INFO = -12
00208       END IF
00209       IF( INFO.NE.0 ) THEN
00210          CALL XERBLA( 'DLAED9', -INFO )
00211          RETURN
00212       END IF
00213 *
00214 *     Quick return if possible
00215 *
00216       IF( K.EQ.0 )
00217      $   RETURN
00218 *
00219 *     Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
00220 *     be computed with high relative accuracy (barring over/underflow).
00221 *     This is a problem on machines without a guard digit in
00222 *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
00223 *     The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
00224 *     which on any of these machines zeros out the bottommost
00225 *     bit of DLAMDA(I) if it is 1; this makes the subsequent
00226 *     subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
00227 *     occurs. On binary machines with a guard digit (almost all
00228 *     machines) it does not change DLAMDA(I) at all. On hexadecimal
00229 *     and decimal machines with a guard digit, it slightly
00230 *     changes the bottommost bits of DLAMDA(I). It does not account
00231 *     for hexadecimal or decimal machines without guard digits
00232 *     (we know of none). We use a subroutine call to compute
00233 *     2*DLAMBDA(I) to prevent optimizing compilers from eliminating
00234 *     this code.
00235 *
00236       DO 10 I = 1, N
00237          DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
00238    10 CONTINUE
00239 *
00240       DO 20 J = KSTART, KSTOP
00241          CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
00242 *
00243 *        If the zero finder fails, the computation is terminated.
00244 *
00245          IF( INFO.NE.0 )
00246      $      GO TO 120
00247    20 CONTINUE
00248 *
00249       IF( K.EQ.1 .OR. K.EQ.2 ) THEN
00250          DO 40 I = 1, K
00251             DO 30 J = 1, K
00252                S( J, I ) = Q( J, I )
00253    30       CONTINUE
00254    40    CONTINUE
00255          GO TO 120
00256       END IF
00257 *
00258 *     Compute updated W.
00259 *
00260       CALL DCOPY( K, W, 1, S, 1 )
00261 *
00262 *     Initialize W(I) = Q(I,I)
00263 *
00264       CALL DCOPY( K, Q, LDQ+1, W, 1 )
00265       DO 70 J = 1, K
00266          DO 50 I = 1, J - 1
00267             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00268    50    CONTINUE
00269          DO 60 I = J + 1, K
00270             W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
00271    60    CONTINUE
00272    70 CONTINUE
00273       DO 80 I = 1, K
00274          W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
00275    80 CONTINUE
00276 *
00277 *     Compute eigenvectors of the modified rank-1 modification.
00278 *
00279       DO 110 J = 1, K
00280          DO 90 I = 1, K
00281             Q( I, J ) = W( I ) / Q( I, J )
00282    90    CONTINUE
00283          TEMP = DNRM2( K, Q( 1, J ), 1 )
00284          DO 100 I = 1, K
00285             S( I, J ) = Q( I, J ) / TEMP
00286   100    CONTINUE
00287   110 CONTINUE
00288 *
00289   120 CONTINUE
00290       RETURN
00291 *
00292 *     End of DLAED9
00293 *
00294       END
 All Files Functions