LAPACK: SRC/dlasq6.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
00002      $                   DNM1, DNM2 )
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       INTEGER            I0, N0, PP
00016       DOUBLE PRECISION   DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
00017 *     ..
00018 *     .. Array Arguments ..
00019       DOUBLE PRECISION   Z( * )
00020 *     ..
00021 *
00022 *  Purpose
00023 *  =======
00024 *
00025 *  DLASQ6 computes one dqd (shift equal to zero) transform in
00026 *  ping-pong form, with protection against underflow and overflow.
00027 *
00028 *  Arguments
00029 *  =========
00030 *
00031 *  I0    (input) INTEGER
00032 *        First index.
00033 *
00034 *  N0    (input) INTEGER
00035 *        Last index.
00036 *
00037 *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
00038 *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
00039 *        an extra argument.
00040 *
00041 *  PP    (input) INTEGER
00042 *        PP=0 for ping, PP=1 for pong.
00043 *
00044 *  DMIN  (output) DOUBLE PRECISION
00045 *        Minimum value of d.
00046 *
00047 *  DMIN1 (output) DOUBLE PRECISION
00048 *        Minimum value of d, excluding D( N0 ).
00049 *
00050 *  DMIN2 (output) DOUBLE PRECISION
00051 *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
00052 *
00053 *  DN    (output) DOUBLE PRECISION
00054 *        d(N0), the last value of d.
00055 *
00056 *  DNM1  (output) DOUBLE PRECISION
00057 *        d(N0-1).
00058 *
00059 *  DNM2  (output) DOUBLE PRECISION
00060 *        d(N0-2).
00061 *
00062 *  =====================================================================
00063 *
00064 *     .. Parameter ..
00065       DOUBLE PRECISION   ZERO
00066       PARAMETER          ( ZERO = 0.0D0 )
00067 *     ..
00068 *     .. Local Scalars ..
00069       INTEGER            J4, J4P2
00070       DOUBLE PRECISION   D, EMIN, SAFMIN, TEMP
00071 *     ..
00072 *     .. External Function ..
00073       DOUBLE PRECISION   DLAMCH
00074       EXTERNAL           DLAMCH
00075 *     ..
00076 *     .. Intrinsic Functions ..
00077       INTRINSIC          MIN
00078 *     ..
00079 *     .. Executable Statements ..
00080 *
00081       IF( ( N0-I0-1 ).LE.0 )
00082      $   RETURN
00083 *
00084       SAFMIN = DLAMCH( 'Safe minimum' )
00085       J4 = 4*I0 + PP - 3
00086       EMIN = Z( J4+4 ) 
00087       D = Z( J4 )
00088       DMIN = D
00089 *
00090       IF( PP.EQ.0 ) THEN
00091          DO 10 J4 = 4*I0, 4*( N0-3 ), 4
00092             Z( J4-2 ) = D + Z( J4-1 ) 
00093             IF( Z( J4-2 ).EQ.ZERO ) THEN
00094                Z( J4 ) = ZERO
00095                D = Z( J4+1 )
00096                DMIN = D
00097                EMIN = ZERO
00098             ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
00099      $               SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
00100                TEMP = Z( J4+1 ) / Z( J4-2 )
00101                Z( J4 ) = Z( J4-1 )*TEMP
00102                D = D*TEMP
00103             ELSE 
00104                Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
00105                D = Z( J4+1 )*( D / Z( J4-2 ) )
00106             END IF
00107             DMIN = MIN( DMIN, D )
00108             EMIN = MIN( EMIN, Z( J4 ) )
00109    10    CONTINUE
00110       ELSE
00111          DO 20 J4 = 4*I0, 4*( N0-3 ), 4
00112             Z( J4-3 ) = D + Z( J4 ) 
00113             IF( Z( J4-3 ).EQ.ZERO ) THEN
00114                Z( J4-1 ) = ZERO
00115                D = Z( J4+2 )
00116                DMIN = D
00117                EMIN = ZERO
00118             ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
00119      $               SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
00120                TEMP = Z( J4+2 ) / Z( J4-3 )
00121                Z( J4-1 ) = Z( J4 )*TEMP
00122                D = D*TEMP
00123             ELSE 
00124                Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
00125                D = Z( J4+2 )*( D / Z( J4-3 ) )
00126             END IF
00127             DMIN = MIN( DMIN, D )
00128             EMIN = MIN( EMIN, Z( J4-1 ) )
00129    20    CONTINUE
00130       END IF
00131 *
00132 *     Unroll last two steps. 
00133 *
00134       DNM2 = D
00135       DMIN2 = DMIN
00136       J4 = 4*( N0-2 ) - PP
00137       J4P2 = J4 + 2*PP - 1
00138       Z( J4-2 ) = DNM2 + Z( J4P2 )
00139       IF( Z( J4-2 ).EQ.ZERO ) THEN
00140          Z( J4 ) = ZERO
00141          DNM1 = Z( J4P2+2 )
00142          DMIN = DNM1
00143          EMIN = ZERO
00144       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
00145      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
00146          TEMP = Z( J4P2+2 ) / Z( J4-2 )
00147          Z( J4 ) = Z( J4P2 )*TEMP
00148          DNM1 = DNM2*TEMP
00149       ELSE
00150          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00151          DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
00152       END IF
00153       DMIN = MIN( DMIN, DNM1 )
00154 *
00155       DMIN1 = DMIN
00156       J4 = J4 + 4
00157       J4P2 = J4 + 2*PP - 1
00158       Z( J4-2 ) = DNM1 + Z( J4P2 )
00159       IF( Z( J4-2 ).EQ.ZERO ) THEN
00160          Z( J4 ) = ZERO
00161          DN = Z( J4P2+2 )
00162          DMIN = DN
00163          EMIN = ZERO
00164       ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
00165      $         SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
00166          TEMP = Z( J4P2+2 ) / Z( J4-2 )
00167          Z( J4 ) = Z( J4P2 )*TEMP
00168          DN = DNM1*TEMP
00169       ELSE
00170          Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
00171          DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
00172       END IF
00173       DMIN = MIN( DMIN, DN )
00174 *
00175       Z( J4+2 ) = DN
00176       Z( 4*N0-PP ) = EMIN
00177       RETURN
00178 *
00179 *     End of DLASQ6
00180 *
00181       END