![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLASD8 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLASD8 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd8.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd8.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd8.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, 00022 * DSIGMA, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER ICOMPQ, INFO, K, LDDIFR 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL D( * ), DIFL( * ), DIFR( LDDIFR, * ), 00029 * $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ), 00030 * $ Z( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> SLASD8 finds the square roots of the roots of the secular equation, 00040 *> as defined by the values in DSIGMA and Z. It makes the appropriate 00041 *> calls to SLASD4, and stores, for each element in D, the distance 00042 *> to its two nearest poles (elements in DSIGMA). It also updates 00043 *> the arrays VF and VL, the first and last components of all the 00044 *> right singular vectors of the original bidiagonal matrix. 00045 *> 00046 *> SLASD8 is called from SLASD6. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] ICOMPQ 00053 *> \verbatim 00054 *> ICOMPQ is INTEGER 00055 *> Specifies whether singular vectors are to be computed in 00056 *> factored form in the calling routine: 00057 *> = 0: Compute singular values only. 00058 *> = 1: Compute singular vectors in factored form as well. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] K 00062 *> \verbatim 00063 *> K is INTEGER 00064 *> The number of terms in the rational function to be solved 00065 *> by SLASD4. K >= 1. 00066 *> \endverbatim 00067 *> 00068 *> \param[out] D 00069 *> \verbatim 00070 *> D is REAL array, dimension ( K ) 00071 *> On output, D contains the updated singular values. 00072 *> \endverbatim 00073 *> 00074 *> \param[in,out] Z 00075 *> \verbatim 00076 *> Z is REAL array, dimension ( K ) 00077 *> On entry, the first K elements of this array contain the 00078 *> components of the deflation-adjusted updating row vector. 00079 *> On exit, Z is updated. 00080 *> \endverbatim 00081 *> 00082 *> \param[in,out] VF 00083 *> \verbatim 00084 *> VF is REAL array, dimension ( K ) 00085 *> On entry, VF contains information passed through DBEDE8. 00086 *> On exit, VF contains the first K components of the first 00087 *> components of all right singular vectors of the bidiagonal 00088 *> matrix. 00089 *> \endverbatim 00090 *> 00091 *> \param[in,out] VL 00092 *> \verbatim 00093 *> VL is REAL array, dimension ( K ) 00094 *> On entry, VL contains information passed through DBEDE8. 00095 *> On exit, VL contains the first K components of the last 00096 *> components of all right singular vectors of the bidiagonal 00097 *> matrix. 00098 *> \endverbatim 00099 *> 00100 *> \param[out] DIFL 00101 *> \verbatim 00102 *> DIFL is REAL array, dimension ( K ) 00103 *> On exit, DIFL(I) = D(I) - DSIGMA(I). 00104 *> \endverbatim 00105 *> 00106 *> \param[out] DIFR 00107 *> \verbatim 00108 *> DIFR is REAL array, 00109 *> dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and 00110 *> dimension ( K ) if ICOMPQ = 0. 00111 *> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not 00112 *> defined and will not be referenced. 00113 *> 00114 *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 00115 *> normalizing factors for the right singular vector matrix. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] LDDIFR 00119 *> \verbatim 00120 *> LDDIFR is INTEGER 00121 *> The leading dimension of DIFR, must be at least K. 00122 *> \endverbatim 00123 *> 00124 *> \param[in,out] DSIGMA 00125 *> \verbatim 00126 *> DSIGMA is REAL array, dimension ( K ) 00127 *> On entry, the first K elements of this array contain the old 00128 *> roots of the deflated updating problem. These are the poles 00129 *> of the secular equation. 00130 *> On exit, the elements of DSIGMA may be very slightly altered 00131 *> in value. 00132 *> \endverbatim 00133 *> 00134 *> \param[out] WORK 00135 *> \verbatim 00136 *> WORK is REAL array, dimension at least 3 * K 00137 *> \endverbatim 00138 *> 00139 *> \param[out] INFO 00140 *> \verbatim 00141 *> INFO is INTEGER 00142 *> = 0: successful exit. 00143 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00144 *> > 0: if INFO = 1, a singular value did not converge 00145 *> \endverbatim 00146 * 00147 * Authors: 00148 * ======== 00149 * 00150 *> \author Univ. of Tennessee 00151 *> \author Univ. of California Berkeley 00152 *> \author Univ. of Colorado Denver 00153 *> \author NAG Ltd. 00154 * 00155 *> \date November 2011 00156 * 00157 *> \ingroup auxOTHERauxiliary 00158 * 00159 *> \par Contributors: 00160 * ================== 00161 *> 00162 *> Ming Gu and Huan Ren, Computer Science Division, University of 00163 *> California at Berkeley, USA 00164 *> 00165 * ===================================================================== 00166 SUBROUTINE SLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR, 00167 $ DSIGMA, WORK, INFO ) 00168 * 00169 * -- LAPACK auxiliary routine (version 3.4.0) -- 00170 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00172 * November 2011 00173 * 00174 * .. Scalar Arguments .. 00175 INTEGER ICOMPQ, INFO, K, LDDIFR 00176 * .. 00177 * .. Array Arguments .. 00178 REAL D( * ), DIFL( * ), DIFR( LDDIFR, * ), 00179 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ), 00180 $ Z( * ) 00181 * .. 00182 * 00183 * ===================================================================== 00184 * 00185 * .. Parameters .. 00186 REAL ONE 00187 PARAMETER ( ONE = 1.0E+0 ) 00188 * .. 00189 * .. Local Scalars .. 00190 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J 00191 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP 00192 * .. 00193 * .. External Subroutines .. 00194 EXTERNAL SCOPY, SLASCL, SLASD4, SLASET, XERBLA 00195 * .. 00196 * .. External Functions .. 00197 REAL SDOT, SLAMC3, SNRM2 00198 EXTERNAL SDOT, SLAMC3, SNRM2 00199 * .. 00200 * .. Intrinsic Functions .. 00201 INTRINSIC ABS, SIGN, SQRT 00202 * .. 00203 * .. Executable Statements .. 00204 * 00205 * Test the input parameters. 00206 * 00207 INFO = 0 00208 * 00209 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00210 INFO = -1 00211 ELSE IF( K.LT.1 ) THEN 00212 INFO = -2 00213 ELSE IF( LDDIFR.LT.K ) THEN 00214 INFO = -9 00215 END IF 00216 IF( INFO.NE.0 ) THEN 00217 CALL XERBLA( 'SLASD8', -INFO ) 00218 RETURN 00219 END IF 00220 * 00221 * Quick return if possible 00222 * 00223 IF( K.EQ.1 ) THEN 00224 D( 1 ) = ABS( Z( 1 ) ) 00225 DIFL( 1 ) = D( 1 ) 00226 IF( ICOMPQ.EQ.1 ) THEN 00227 DIFL( 2 ) = ONE 00228 DIFR( 1, 2 ) = ONE 00229 END IF 00230 RETURN 00231 END IF 00232 * 00233 * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can 00234 * be computed with high relative accuracy (barring over/underflow). 00235 * This is a problem on machines without a guard digit in 00236 * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). 00237 * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), 00238 * which on any of these machines zeros out the bottommost 00239 * bit of DSIGMA(I) if it is 1; this makes the subsequent 00240 * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation 00241 * occurs. On binary machines with a guard digit (almost all 00242 * machines) it does not change DSIGMA(I) at all. On hexadecimal 00243 * and decimal machines with a guard digit, it slightly 00244 * changes the bottommost bits of DSIGMA(I). It does not account 00245 * for hexadecimal or decimal machines without guard digits 00246 * (we know of none). We use a subroutine call to compute 00247 * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating 00248 * this code. 00249 * 00250 DO 10 I = 1, K 00251 DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) 00252 10 CONTINUE 00253 * 00254 * Book keeping. 00255 * 00256 IWK1 = 1 00257 IWK2 = IWK1 + K 00258 IWK3 = IWK2 + K 00259 IWK2I = IWK2 - 1 00260 IWK3I = IWK3 - 1 00261 * 00262 * Normalize Z. 00263 * 00264 RHO = SNRM2( K, Z, 1 ) 00265 CALL SLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO ) 00266 RHO = RHO*RHO 00267 * 00268 * Initialize WORK(IWK3). 00269 * 00270 CALL SLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K ) 00271 * 00272 * Compute the updated singular values, the arrays DIFL, DIFR, 00273 * and the updated Z. 00274 * 00275 DO 40 J = 1, K 00276 CALL SLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ), 00277 $ WORK( IWK2 ), INFO ) 00278 * 00279 * If the root finder fails, the computation is terminated. 00280 * 00281 IF( INFO.NE.0 ) THEN 00282 CALL XERBLA( 'SLASD4', -INFO ) 00283 RETURN 00284 END IF 00285 WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J ) 00286 DIFL( J ) = -WORK( J ) 00287 DIFR( J, 1 ) = -WORK( J+1 ) 00288 DO 20 I = 1, J - 1 00289 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* 00290 $ WORK( IWK2I+I ) / ( DSIGMA( I )- 00291 $ DSIGMA( J ) ) / ( DSIGMA( I )+ 00292 $ DSIGMA( J ) ) 00293 20 CONTINUE 00294 DO 30 I = J + 1, K 00295 WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )* 00296 $ WORK( IWK2I+I ) / ( DSIGMA( I )- 00297 $ DSIGMA( J ) ) / ( DSIGMA( I )+ 00298 $ DSIGMA( J ) ) 00299 30 CONTINUE 00300 40 CONTINUE 00301 * 00302 * Compute updated Z. 00303 * 00304 DO 50 I = 1, K 00305 Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) ) 00306 50 CONTINUE 00307 * 00308 * Update VF and VL. 00309 * 00310 DO 80 J = 1, K 00311 DIFLJ = DIFL( J ) 00312 DJ = D( J ) 00313 DSIGJ = -DSIGMA( J ) 00314 IF( J.LT.K ) THEN 00315 DIFRJ = -DIFR( J, 1 ) 00316 DSIGJP = -DSIGMA( J+1 ) 00317 END IF 00318 WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ ) 00319 DO 60 I = 1, J - 1 00320 WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ ) 00321 $ / ( DSIGMA( I )+DJ ) 00322 60 CONTINUE 00323 DO 70 I = J + 1, K 00324 WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ ) 00325 $ / ( DSIGMA( I )+DJ ) 00326 70 CONTINUE 00327 TEMP = SNRM2( K, WORK, 1 ) 00328 WORK( IWK2I+J ) = SDOT( K, WORK, 1, VF, 1 ) / TEMP 00329 WORK( IWK3I+J ) = SDOT( K, WORK, 1, VL, 1 ) / TEMP 00330 IF( ICOMPQ.EQ.1 ) THEN 00331 DIFR( J, 2 ) = TEMP 00332 END IF 00333 80 CONTINUE 00334 * 00335 CALL SCOPY( K, WORK( IWK2 ), 1, VF, 1 ) 00336 CALL SCOPY( K, WORK( IWK3 ), 1, VL, 1 ) 00337 * 00338 RETURN 00339 * 00340 * End of SLASD8 00341 * 00342 END 00343