00001 SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, 00002 $ DSIGMA, WORK, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.3.0) -- 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 2010 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER ICOMPQ, INFO, K, LDDIFR 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ), 00014 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ), 00015 $ Z( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLASD8 finds the square roots of the roots of the secular equation, 00022 * as defined by the values in DSIGMA and Z. It makes the appropriate 00023 * calls to DLASD4, and stores, for each element in D, the distance 00024 * to its two nearest poles (elements in DSIGMA). It also updates 00025 * the arrays VF and VL, the first and last components of all the 00026 * right singular vectors of the original bidiagonal matrix. 00027 * 00028 * DLASD8 is called from DLASD6. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * ICOMPQ (input) INTEGER 00034 * Specifies whether singular vectors are to be computed in 00035 * factored form in the calling routine: 00036 * = 0: Compute singular values only. 00037 * = 1: Compute singular vectors in factored form as well. 00038 * 00039 * K (input) INTEGER 00040 * The number of terms in the rational function to be solved 00041 * by DLASD4. K >= 1. 00042 * 00043 * D (output) DOUBLE PRECISION array, dimension ( K ) 00044 * On output, D contains the updated singular values. 00045 * 00046 * Z (input/output) DOUBLE PRECISION array, dimension ( K ) 00047 * On entry, the first K elements of this array contain the 00048 * components of the deflation-adjusted updating row vector. 00049 * On exit, Z is updated. 00050 * 00051 * VF (input/output) DOUBLE PRECISION array, dimension ( K ) 00052 * On entry, VF contains information passed through DBEDE8. 00053 * On exit, VF contains the first K components of the first 00054 * components of all right singular vectors of the bidiagonal 00055 * matrix. 00056 * 00057 * VL (input/output) DOUBLE PRECISION array, dimension ( K ) 00058 * On entry, VL contains information passed through DBEDE8. 00059 * On exit, VL contains the first K components of the last 00060 * components of all right singular vectors of the bidiagonal 00061 * matrix. 00062 * 00063 * DIFL (output) DOUBLE PRECISION array, dimension ( K ) 00064 * On exit, DIFL(I) = D(I) - DSIGMA(I). 00065 * 00066 * DIFR (output) DOUBLE PRECISION array, 00067 * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and 00068 * dimension ( K ) if ICOMPQ = 0. 00069 * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not 00070 * defined and will not be referenced. 00071 * 00072 * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 00073 * normalizing factors for the right singular vector matrix. 00074 * 00075 * LDDIFR (input) INTEGER 00076 * The leading dimension of DIFR, must be at least K. 00077 * 00078 * DSIGMA (input/output) DOUBLE PRECISION array, dimension ( K ) 00079 * On entry, the first K elements of this array contain the old 00080 * roots of the deflated updating problem. These are the poles 00081 * of the secular equation. 00082 * On exit, the elements of DSIGMA may be very slightly altered 00083 * in value. 00084 * 00085 * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K 00086 * 00087 * INFO (output) INTEGER 00088 * = 0: successful exit. 00089 * < 0: if INFO = -i, the i-th argument had an illegal value. 00090 * > 0: if INFO = 1, a singular value did not converge 00091 * 00092 * Further Details 00093 * =============== 00094 * 00095 * Based on contributions by 00096 * Ming Gu and Huan Ren, Computer Science Division, University of 00097 * California at Berkeley, USA 00098 * 00099 * ===================================================================== 00100 * 00101 * .. Parameters .. 00102 DOUBLE PRECISION ONE 00103 PARAMETER ( ONE = 1.0D+0 ) 00104 * .. 00105 * .. Local Scalars .. 00106 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J 00107 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP 00108 * .. 00109 * .. External Subroutines .. 00110 EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA 00111 * .. 00112 * .. External Functions .. 00113 DOUBLE PRECISION DDOT, DLAMC3, DNRM2 00114 EXTERNAL DDOT, DLAMC3, DNRM2 00115 * .. 00116 * .. Intrinsic Functions .. 00117 INTRINSIC ABS, SIGN, SQRT 00118 * .. 00119 * .. Executable Statements .. 00120 * 00121 * Test the input parameters. 00122 * 00123 INFO = 0 00124 * 00125 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00126 INFO = -1 00127 ELSE IF( K.LT.1 ) THEN 00128 INFO = -2 00129 ELSE IF( LDDIFR.LT.K ) THEN 00130 INFO = -9 00131 END IF 00132 IF( INFO.NE.0 ) THEN 00133 CALL XERBLA( 'DLASD8', -INFO ) 00134 RETURN 00135 END IF 00136 * 00137 * Quick return if possible 00138 * 00139 IF( K.EQ.1 ) THEN 00140 D( 1 ) = ABS( Z( 1 ) ) 00141 DIFL( 1 ) = D( 1 ) 00142 IF( ICOMPQ.EQ.1 ) THEN 00143 DIFL( 2 ) = ONE 00144 DIFR( 1, 2 ) = ONE 00145 END IF 00146 RETURN 00147 END IF 00148 * 00149 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can 00150 * be computed with high relative accuracy (barring over/underflow). 00151 * This is a problem on machines without a guard digit in 00152 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 00153 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), 00154 * which on any of these machines zeros out the bottommost 00155 * bit of DSIGMA(I) if it is 1; this makes the subsequent 00156 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation 00157 * occurs. On binary machines with a guard digit (almost all 00158 * machines) it does not change DSIGMA(I) at all. On hexadecimal 00159 * and decimal machines with a guard digit, it slightly 00160 * changes the bottommost bits of DSIGMA(I). It does not account 00161 * for hexadecimal or decimal machines without guard digits 00162 * (we know of none). We use a subroutine call to compute 00163 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating 00164 * this code. 00165 * 00166 DO 10 I = 1, K 00167 DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) 00168 10 CONTINUE 00169 * 00170 * Book keeping. 00171 * 00172 IWK1 = 1 00173 IWK2 = IWK1 + K 00174 IWK3 = IWK2 + K 00175 IWK2I = IWK2 - 1 00176 IWK3I = IWK3 - 1 00177 * 00178 * Normalize Z. 00179 * 00180 RHO = DNRM2( K, Z, 1 ) 00181 CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO ) 00182 RHO = RHO*RHO 00183 * 00184 * Initialize WORK(IWK3). 00185 * 00186 CALL DLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K ) 00187 * 00188 * Compute the updated singular values, the arrays DIFL, DIFR, 00189 * and the updated Z. 00190 * 00191 DO 40 J = 1, K 00192 CALL DLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ), 00193 $ WORK( IWK2 ), INFO ) 00194 * 00195 * If the root finder fails, the computation is terminated. 00196 * 00197 IF( INFO.NE.0 ) THEN 00198 CALL XERBLA( 'DLASD4', -INFO ) 00199 RETURN 00200 END IF 00201 WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J ) 00202 DIFL( J ) = -WORK( J ) 00203 DIFR( J, 1 ) = -WORK( J+1 ) 00204 DO 20 I = 1, J - 1 00205 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* 00206 $ WORK( IWK2I+I ) / ( DSIGMA( I )- 00207 $ DSIGMA( J ) ) / ( DSIGMA( I )+ 00208 $ DSIGMA( J ) ) 00209 20 CONTINUE 00210 DO 30 I = J + 1, K 00211 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* 00212 $ WORK( IWK2I+I ) / ( DSIGMA( I )- 00213 $ DSIGMA( J ) ) / ( DSIGMA( I )+ 00214 $ DSIGMA( J ) ) 00215 30 CONTINUE 00216 40 CONTINUE 00217 * 00218 * Compute updated Z. 00219 * 00220 DO 50 I = 1, K 00221 Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) ) 00222 50 CONTINUE 00223 * 00224 * Update VF and VL. 00225 * 00226 DO 80 J = 1, K 00227 DIFLJ = DIFL( J ) 00228 DJ = D( J ) 00229 DSIGJ = -DSIGMA( J ) 00230 IF( J.LT.K ) THEN 00231 DIFRJ = -DIFR( J, 1 ) 00232 DSIGJP = -DSIGMA( J+1 ) 00233 END IF 00234 WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ ) 00235 DO 60 I = 1, J - 1 00236 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ ) 00237 $ / ( DSIGMA( I )+DJ ) 00238 60 CONTINUE 00239 DO 70 I = J + 1, K 00240 WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ ) 00241 $ / ( DSIGMA( I )+DJ ) 00242 70 CONTINUE 00243 TEMP = DNRM2( K, WORK, 1 ) 00244 WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP 00245 WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP 00246 IF( ICOMPQ.EQ.1 ) THEN 00247 DIFR( J, 2 ) = TEMP 00248 END IF 00249 80 CONTINUE 00250 * 00251 CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 ) 00252 CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 ) 00253 * 00254 RETURN 00255 * 00256 * End of DLASD8 00257 * 00258 END 00259
1.7.2