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