LAPACK: SRC/slasq5.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
00002      $                   DNM1, DNM2, IEEE )
00003 *
00004 *  -- LAPACK routine (version 3.2)                                    --
00005 *
00006 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
00007 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
00008 *  -- Berkeley                                                        --
00009 *  -- November 2008                                                   --
00010 *
00011 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00012 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00013 *
00014 *     .. Scalar Arguments ..
00015       LOGICAL            IEEE
00016       INTEGER            I0, N0, PP
00017       REAL               DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
00018 *     ..
00019 *     .. Array Arguments ..
00020       REAL               Z( * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  SLASQ5 computes one dqds transform in ping-pong form, one
00027 *  version for IEEE machines another for non IEEE machines.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  I0    (input) INTEGER
00033 *        First index.
00034 *
00035 *  N0    (input) INTEGER
00036 *        Last index.
00037 *
00038 *  Z     (input) REAL array, dimension ( 4*N )
00039 *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
00040 *        an extra argument.
00041 *
00042 *  PP    (input) INTEGER
00043 *        PP=0 for ping, PP=1 for pong.
00044 *
00045 *  TAU   (input) REAL
00046 *        This is the shift.
00047 *
00048 *  DMIN  (output) REAL
00049 *        Minimum value of d.
00050 *
00051 *  DMIN1 (output) REAL
00052 *        Minimum value of d, excluding D( N0 ).
00053 *
00054 *  DMIN2 (output) REAL
00055 *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
00056 *
00057 *  DN    (output) REAL
00058 *        d(N0), the last value of d.
00059 *
00060 *  DNM1  (output) REAL
00061 *        d(N0-1).
00062 *
00063 *  DNM2  (output) REAL
00064 *        d(N0-2).
00065 *
00066 *  IEEE  (input) LOGICAL
00067 *        Flag for IEEE or non IEEE arithmetic.
00068 *
00069 *  =====================================================================
00070 *
00071 *     .. Parameter ..
00072       REAL               ZERO
00073       PARAMETER          ( ZERO = 0.0E0 )
00074 *     ..
00075 *     .. Local Scalars ..
00076       INTEGER            J4, J4P2
00077       REAL               D, EMIN, TEMP
00078 *     ..
00079 *     .. Intrinsic Functions ..
00080       INTRINSIC          MIN
00081 *     ..
00082 *     .. Executable Statements ..
00083 *
00084       IF( ( N0-I0-1 ).LE.0 )
00085      $   RETURN
00086 *
00087       J4 = 4*I0 + PP - 3
00088       EMIN = Z( J4+4 )
00089       D = Z( J4 ) - TAU
00090       DMIN = D
00091       DMIN1 = -Z( J4 )
00092 *
00093       IF( IEEE ) THEN
00094 *
00095 *        Code for IEEE arithmetic.
00096 *
00097          IF( PP.EQ.0 ) THEN
00098             DO 10 J4 = 4*I0, 4*( N0-3 ), 4
00099                Z( J4-2 ) = D + Z( J4-1 )
00100                TEMP = Z( J4+1 ) / Z( J4-2 )
00101                D = D*TEMP - TAU
00102                DMIN = MIN( DMIN, D )
00103                Z( J4 ) = Z( J4-1 )*TEMP
00104                EMIN = MIN( Z( J4 ), EMIN )
00105    10       CONTINUE
00106          ELSE
00107             DO 20 J4 = 4*I0, 4*( N0-3 ), 4
00108                Z( J4-3 ) = D + Z( J4 )
00109                TEMP = Z( J4+2 ) / Z( J4-3 )
00110                D = D*TEMP - TAU
00111                DMIN = MIN( DMIN, D )
00112                Z( J4-1 ) = Z( J4 )*TEMP
00113                EMIN = MIN( Z( J4-1 ), EMIN )
00114    20       CONTINUE
00115          END IF
00116 *
00117 *        Unroll last two steps.
00118 *
00119          DNM2 = D
00120          DMIN2 = DMIN
00121          J4 = 4*( N0-2 ) - PP
00122          J4P2 = J4 + 2*PP - 1
00123          Z( J4-2 ) = DNM2 + Z( J4P2 )
00124          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00125          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00126          DMIN = MIN( DMIN, DNM1 )
00127 *
00128          DMIN1 = DMIN
00129          J4 = J4 + 4
00130          J4P2 = J4 + 2*PP - 1
00131          Z( J4-2 ) = DNM1 + Z( J4P2 )
00132          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00133          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00134          DMIN = MIN( DMIN, DN )
00135 *
00136       ELSE
00137 *
00138 *        Code for non IEEE arithmetic.
00139 *
00140          IF( PP.EQ.0 ) THEN
00141             DO 30 J4 = 4*I0, 4*( N0-3 ), 4
00142                Z( J4-2 ) = D + Z( J4-1 )
00143                IF( D.LT.ZERO ) THEN
00144                   RETURN
00145                ELSE
00146                   Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
00147                   D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
00148                END IF
00149                DMIN = MIN( DMIN, D )
00150                EMIN = MIN( EMIN, Z( J4 ) )
00151    30       CONTINUE
00152          ELSE
00153             DO 40 J4 = 4*I0, 4*( N0-3 ), 4
00154                Z( J4-3 ) = D + Z( J4 )
00155                IF( D.LT.ZERO ) THEN
00156                   RETURN
00157                ELSE
00158                   Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
00159                   D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
00160                END IF
00161                DMIN = MIN( DMIN, D )
00162                EMIN = MIN( EMIN, Z( J4-1 ) )
00163    40       CONTINUE
00164          END IF
00165 *
00166 *        Unroll last two steps.
00167 *
00168          DNM2 = D
00169          DMIN2 = DMIN
00170          J4 = 4*( N0-2 ) - PP
00171          J4P2 = J4 + 2*PP - 1
00172          Z( J4-2 ) = DNM2 + Z( J4P2 )
00173          IF( DNM2.LT.ZERO ) THEN
00174             RETURN
00175          ELSE
00176             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00177             DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
00178          END IF
00179          DMIN = MIN( DMIN, DNM1 )
00180 *
00181          DMIN1 = DMIN
00182          J4 = J4 + 4
00183          J4P2 = J4 + 2*PP - 1
00184          Z( J4-2 ) = DNM1 + Z( J4P2 )
00185          IF( DNM1.LT.ZERO ) THEN
00186             RETURN
00187          ELSE
00188             Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00189             DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
00190          END IF
00191          DMIN = MIN( DMIN, DN )
00192 *
00193       END IF
00194 *
00195       Z( J4+2 ) = DN
00196       Z( 4*N0-PP ) = EMIN
00197       RETURN
00198 *
00199 *     End of SLASQ5
00200 *
00201       END