LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slasq5.f
Go to the documentation of this file.
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
 All Files Functions