![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLASQ5 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLASQ5 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq5.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq5.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq5.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN, 00022 * DNM1, DNM2, IEEE, EPS ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL IEEE 00026 * INTEGER I0, N0, PP 00027 * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, SIGMA, EPS 00028 * .. 00029 * .. Array Arguments .. 00030 * DOUBLE PRECISION Z( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DLASQ5 computes one dqds transform in ping-pong form, one 00040 *> version for IEEE machines another for non IEEE machines. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] I0 00047 *> \verbatim 00048 *> I0 is INTEGER 00049 *> First index. 00050 *> \endverbatim 00051 *> 00052 *> \param[in] N0 00053 *> \verbatim 00054 *> N0 is INTEGER 00055 *> Last index. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] Z 00059 *> \verbatim 00060 *> Z is DOUBLE PRECISION array, dimension ( 4*N ) 00061 *> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid 00062 *> an extra argument. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] PP 00066 *> \verbatim 00067 *> PP is INTEGER 00068 *> PP=0 for ping, PP=1 for pong. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] TAU 00072 *> \verbatim 00073 *> TAU is DOUBLE PRECISION 00074 *> This is the shift. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] SIGMA 00078 *> \verbatim 00079 *> SIGMA is DOUBLE PRECISION 00080 *> This is the accumulated shift up to this step. 00081 *> \endverbatim 00082 *> 00083 *> \param[out] DMIN 00084 *> \verbatim 00085 *> DMIN is DOUBLE PRECISION 00086 *> Minimum value of d. 00087 *> \endverbatim 00088 *> 00089 *> \param[out] DMIN1 00090 *> \verbatim 00091 *> DMIN1 is DOUBLE PRECISION 00092 *> Minimum value of d, excluding D( N0 ). 00093 *> \endverbatim 00094 *> 00095 *> \param[out] DMIN2 00096 *> \verbatim 00097 *> DMIN2 is DOUBLE PRECISION 00098 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 00099 *> \endverbatim 00100 *> 00101 *> \param[out] DN 00102 *> \verbatim 00103 *> DN is DOUBLE PRECISION 00104 *> d(N0), the last value of d. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] DNM1 00108 *> \verbatim 00109 *> DNM1 is DOUBLE PRECISION 00110 *> d(N0-1). 00111 *> \endverbatim 00112 *> 00113 *> \param[out] DNM2 00114 *> \verbatim 00115 *> DNM2 is DOUBLE PRECISION 00116 *> d(N0-2). 00117 *> \endverbatim 00118 *> 00119 *> \param[in] IEEE 00120 *> \verbatim 00121 *> IEEE is LOGICAL 00122 *> Flag for IEEE or non IEEE arithmetic. 00123 *> \endverbatim 00124 * 00125 *> \param[in] EPS 00126 *> \verbatim 00127 *> EPS is DOUBLE PRECISION 00128 *> This is the value of epsilon used. 00129 *> \endverbatim 00130 *> 00131 * Authors: 00132 * ======== 00133 * 00134 *> \author Univ. of Tennessee 00135 *> \author Univ. of California Berkeley 00136 *> \author Univ. of Colorado Denver 00137 *> \author NAG Ltd. 00138 * 00139 *> \date April 2012 00140 * 00141 *> \ingroup auxOTHERcomputational 00142 * 00143 * ===================================================================== 00144 SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, 00145 $ DN, DNM1, DNM2, IEEE, EPS ) 00146 * 00147 * -- LAPACK computational routine (version 3.4.1) -- 00148 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00149 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00150 * April 2012 00151 * 00152 * .. Scalar Arguments .. 00153 LOGICAL IEEE 00154 INTEGER I0, N0, PP 00155 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU, 00156 $ SIGMA, EPS 00157 * .. 00158 * .. Array Arguments .. 00159 DOUBLE PRECISION Z( * ) 00160 * .. 00161 * 00162 * ===================================================================== 00163 * 00164 * .. Parameter .. 00165 DOUBLE PRECISION ZERO, HALF 00166 PARAMETER ( ZERO = 0.0D0, HALF = 0.5 ) 00167 * .. 00168 * .. Local Scalars .. 00169 INTEGER J4, J4P2 00170 DOUBLE PRECISION D, EMIN, TEMP, DTHRESH 00171 * .. 00172 * .. Intrinsic Functions .. 00173 INTRINSIC MIN 00174 * .. 00175 * .. Executable Statements .. 00176 * 00177 IF( ( N0-I0-1 ).LE.0 ) 00178 $ RETURN 00179 * 00180 DTHRESH = EPS*(SIGMA+TAU) 00181 IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO 00182 IF( TAU.NE.ZERO ) THEN 00183 J4 = 4*I0 + PP - 3 00184 EMIN = Z( J4+4 ) 00185 D = Z( J4 ) - TAU 00186 DMIN = D 00187 DMIN1 = -Z( J4 ) 00188 * 00189 IF( IEEE ) THEN 00190 * 00191 * Code for IEEE arithmetic. 00192 * 00193 IF( PP.EQ.0 ) THEN 00194 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 00195 Z( J4-2 ) = D + Z( J4-1 ) 00196 TEMP = Z( J4+1 ) / Z( J4-2 ) 00197 D = D*TEMP - TAU 00198 DMIN = MIN( DMIN, D ) 00199 Z( J4 ) = Z( J4-1 )*TEMP 00200 EMIN = MIN( Z( J4 ), EMIN ) 00201 10 CONTINUE 00202 ELSE 00203 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 00204 Z( J4-3 ) = D + Z( J4 ) 00205 TEMP = Z( J4+2 ) / Z( J4-3 ) 00206 D = D*TEMP - TAU 00207 DMIN = MIN( DMIN, D ) 00208 Z( J4-1 ) = Z( J4 )*TEMP 00209 EMIN = MIN( Z( J4-1 ), EMIN ) 00210 20 CONTINUE 00211 END IF 00212 * 00213 * Unroll last two steps. 00214 * 00215 DNM2 = D 00216 DMIN2 = DMIN 00217 J4 = 4*( N0-2 ) - PP 00218 J4P2 = J4 + 2*PP - 1 00219 Z( J4-2 ) = DNM2 + Z( J4P2 ) 00220 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00221 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 00222 DMIN = MIN( DMIN, DNM1 ) 00223 * 00224 DMIN1 = DMIN 00225 J4 = J4 + 4 00226 J4P2 = J4 + 2*PP - 1 00227 Z( J4-2 ) = DNM1 + Z( J4P2 ) 00228 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00229 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 00230 DMIN = MIN( DMIN, DN ) 00231 * 00232 ELSE 00233 * 00234 * Code for non IEEE arithmetic. 00235 * 00236 IF( PP.EQ.0 ) THEN 00237 DO 30 J4 = 4*I0, 4*( N0-3 ), 4 00238 Z( J4-2 ) = D + Z( J4-1 ) 00239 IF( D.LT.ZERO ) THEN 00240 RETURN 00241 ELSE 00242 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 00243 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 00244 END IF 00245 DMIN = MIN( DMIN, D ) 00246 EMIN = MIN( EMIN, Z( J4 ) ) 00247 30 CONTINUE 00248 ELSE 00249 DO 40 J4 = 4*I0, 4*( N0-3 ), 4 00250 Z( J4-3 ) = D + Z( J4 ) 00251 IF( D.LT.ZERO ) THEN 00252 RETURN 00253 ELSE 00254 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 00255 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 00256 END IF 00257 DMIN = MIN( DMIN, D ) 00258 EMIN = MIN( EMIN, Z( J4-1 ) ) 00259 40 CONTINUE 00260 END IF 00261 * 00262 * Unroll last two steps. 00263 * 00264 DNM2 = D 00265 DMIN2 = DMIN 00266 J4 = 4*( N0-2 ) - PP 00267 J4P2 = J4 + 2*PP - 1 00268 Z( J4-2 ) = DNM2 + Z( J4P2 ) 00269 IF( DNM2.LT.ZERO ) THEN 00270 RETURN 00271 ELSE 00272 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00273 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 00274 END IF 00275 DMIN = MIN( DMIN, DNM1 ) 00276 * 00277 DMIN1 = DMIN 00278 J4 = J4 + 4 00279 J4P2 = J4 + 2*PP - 1 00280 Z( J4-2 ) = DNM1 + Z( J4P2 ) 00281 IF( DNM1.LT.ZERO ) THEN 00282 RETURN 00283 ELSE 00284 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00285 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 00286 END IF 00287 DMIN = MIN( DMIN, DN ) 00288 * 00289 END IF 00290 ELSE 00291 * This is the version that sets d's to zero if they are small enough 00292 J4 = 4*I0 + PP - 3 00293 EMIN = Z( J4+4 ) 00294 D = Z( J4 ) - TAU 00295 DMIN = D 00296 DMIN1 = -Z( J4 ) 00297 IF( IEEE ) THEN 00298 * 00299 * Code for IEEE arithmetic. 00300 * 00301 IF( PP.EQ.0 ) THEN 00302 DO 50 J4 = 4*I0, 4*( N0-3 ), 4 00303 Z( J4-2 ) = D + Z( J4-1 ) 00304 TEMP = Z( J4+1 ) / Z( J4-2 ) 00305 D = D*TEMP - TAU 00306 IF( D.LT.DTHRESH ) D = ZERO 00307 DMIN = MIN( DMIN, D ) 00308 Z( J4 ) = Z( J4-1 )*TEMP 00309 EMIN = MIN( Z( J4 ), EMIN ) 00310 50 CONTINUE 00311 ELSE 00312 DO 60 J4 = 4*I0, 4*( N0-3 ), 4 00313 Z( J4-3 ) = D + Z( J4 ) 00314 TEMP = Z( J4+2 ) / Z( J4-3 ) 00315 D = D*TEMP - TAU 00316 IF( D.LT.DTHRESH ) D = ZERO 00317 DMIN = MIN( DMIN, D ) 00318 Z( J4-1 ) = Z( J4 )*TEMP 00319 EMIN = MIN( Z( J4-1 ), EMIN ) 00320 60 CONTINUE 00321 END IF 00322 * 00323 * Unroll last two steps. 00324 * 00325 DNM2 = D 00326 DMIN2 = DMIN 00327 J4 = 4*( N0-2 ) - PP 00328 J4P2 = J4 + 2*PP - 1 00329 Z( J4-2 ) = DNM2 + Z( J4P2 ) 00330 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00331 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 00332 DMIN = MIN( DMIN, DNM1 ) 00333 * 00334 DMIN1 = DMIN 00335 J4 = J4 + 4 00336 J4P2 = J4 + 2*PP - 1 00337 Z( J4-2 ) = DNM1 + Z( J4P2 ) 00338 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00339 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 00340 DMIN = MIN( DMIN, DN ) 00341 * 00342 ELSE 00343 * 00344 * Code for non IEEE arithmetic. 00345 * 00346 IF( PP.EQ.0 ) THEN 00347 DO 70 J4 = 4*I0, 4*( N0-3 ), 4 00348 Z( J4-2 ) = D + Z( J4-1 ) 00349 IF( D.LT.ZERO ) THEN 00350 RETURN 00351 ELSE 00352 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 00353 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 00354 END IF 00355 IF( D.LT.DTHRESH) D = ZERO 00356 DMIN = MIN( DMIN, D ) 00357 EMIN = MIN( EMIN, Z( J4 ) ) 00358 70 CONTINUE 00359 ELSE 00360 DO 80 J4 = 4*I0, 4*( N0-3 ), 4 00361 Z( J4-3 ) = D + Z( J4 ) 00362 IF( D.LT.ZERO ) THEN 00363 RETURN 00364 ELSE 00365 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 00366 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 00367 END IF 00368 IF( D.LT.DTHRESH) D = ZERO 00369 DMIN = MIN( DMIN, D ) 00370 EMIN = MIN( EMIN, Z( J4-1 ) ) 00371 80 CONTINUE 00372 END IF 00373 * 00374 * Unroll last two steps. 00375 * 00376 DNM2 = D 00377 DMIN2 = DMIN 00378 J4 = 4*( N0-2 ) - PP 00379 J4P2 = J4 + 2*PP - 1 00380 Z( J4-2 ) = DNM2 + Z( J4P2 ) 00381 IF( DNM2.LT.ZERO ) THEN 00382 RETURN 00383 ELSE 00384 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00385 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 00386 END IF 00387 DMIN = MIN( DMIN, DNM1 ) 00388 * 00389 DMIN1 = DMIN 00390 J4 = J4 + 4 00391 J4P2 = J4 + 2*PP - 1 00392 Z( J4-2 ) = DNM1 + Z( J4P2 ) 00393 IF( DNM1.LT.ZERO ) THEN 00394 RETURN 00395 ELSE 00396 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 00397 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 00398 END IF 00399 DMIN = MIN( DMIN, DN ) 00400 * 00401 END IF 00402 END IF 00403 * 00404 Z( J4+2 ) = DN 00405 Z( 4*N0-PP ) = EMIN 00406 RETURN 00407 * 00408 * End of DLASQ5 00409 * 00410 END