|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
