![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAED9 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAED9 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed9.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed9.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed9.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAED9( 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 * REAL RHO 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), 00030 * $ W( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SLAED9 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 SLAED4 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 *> SLAED4. 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 SLAED9( 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 REAL RHO 00167 * .. 00168 * .. Array Arguments .. 00169 REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), 00170 $ W( * ) 00171 * .. 00172 * 00173 * ===================================================================== 00174 * 00175 * .. Local Scalars .. 00176 INTEGER I, J 00177 REAL TEMP 00178 * .. 00179 * .. External Functions .. 00180 REAL SLAMC3, SNRM2 00181 EXTERNAL SLAMC3, SNRM2 00182 * .. 00183 * .. External Subroutines .. 00184 EXTERNAL SCOPY, SLAED4, 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( 'SLAED9', -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 ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) 00238 10 CONTINUE 00239 * 00240 DO 20 J = KSTART, KSTOP 00241 CALL SLAED4( 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 SCOPY( K, W, 1, S, 1 ) 00261 * 00262 * Initialize W(I) = Q(I,I) 00263 * 00264 CALL SCOPY( 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 = SNRM2( 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 SLAED9 00293 * 00294 END