176 lines
4.6 KiB
FortranFixed
176 lines
4.6 KiB
FortranFixed
SUBROUTINE SLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN,
|
|
$ DNM1, DNM2 )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER I0, N0, PP
|
|
REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2
|
|
* ..
|
|
* .. Array Arguments ..
|
|
REAL Z( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* SLASQ6 computes one dqd (shift equal to zero) transform in
|
|
* ping-pong form, with protection against underflow and overflow.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* I0 (input) INTEGER
|
|
* First index.
|
|
*
|
|
* N0 (input) INTEGER
|
|
* Last index.
|
|
*
|
|
* Z (input) REAL array, dimension ( 4*N )
|
|
* Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
|
|
* an extra argument.
|
|
*
|
|
* PP (input) INTEGER
|
|
* PP=0 for ping, PP=1 for pong.
|
|
*
|
|
* DMIN (output) REAL
|
|
* Minimum value of d.
|
|
*
|
|
* DMIN1 (output) REAL
|
|
* Minimum value of d, excluding D( N0 ).
|
|
*
|
|
* DMIN2 (output) REAL
|
|
* Minimum value of d, excluding D( N0 ) and D( N0-1 ).
|
|
*
|
|
* DN (output) REAL
|
|
* d(N0), the last value of d.
|
|
*
|
|
* DNM1 (output) REAL
|
|
* d(N0-1).
|
|
*
|
|
* DNM2 (output) REAL
|
|
* d(N0-2).
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameter ..
|
|
REAL ZERO
|
|
PARAMETER ( ZERO = 0.0E0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER J4, J4P2
|
|
REAL D, EMIN, SAFMIN, TEMP
|
|
* ..
|
|
* .. External Function ..
|
|
REAL SLAMCH
|
|
EXTERNAL SLAMCH
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
IF( ( N0-I0-1 ).LE.0 )
|
|
$ RETURN
|
|
*
|
|
SAFMIN = SLAMCH( 'Safe minimum' )
|
|
J4 = 4*I0 + PP - 3
|
|
EMIN = Z( J4+4 )
|
|
D = Z( J4 )
|
|
DMIN = D
|
|
*
|
|
IF( PP.EQ.0 ) THEN
|
|
DO 10 J4 = 4*I0, 4*( N0-3 ), 4
|
|
Z( J4-2 ) = D + Z( J4-1 )
|
|
IF( Z( J4-2 ).EQ.ZERO ) THEN
|
|
Z( J4 ) = ZERO
|
|
D = Z( J4+1 )
|
|
DMIN = D
|
|
EMIN = ZERO
|
|
ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND.
|
|
$ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN
|
|
TEMP = Z( J4+1 ) / Z( J4-2 )
|
|
Z( J4 ) = Z( J4-1 )*TEMP
|
|
D = D*TEMP
|
|
ELSE
|
|
Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
|
|
D = Z( J4+1 )*( D / Z( J4-2 ) )
|
|
END IF
|
|
DMIN = MIN( DMIN, D )
|
|
EMIN = MIN( EMIN, Z( J4 ) )
|
|
10 CONTINUE
|
|
ELSE
|
|
DO 20 J4 = 4*I0, 4*( N0-3 ), 4
|
|
Z( J4-3 ) = D + Z( J4 )
|
|
IF( Z( J4-3 ).EQ.ZERO ) THEN
|
|
Z( J4-1 ) = ZERO
|
|
D = Z( J4+2 )
|
|
DMIN = D
|
|
EMIN = ZERO
|
|
ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND.
|
|
$ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN
|
|
TEMP = Z( J4+2 ) / Z( J4-3 )
|
|
Z( J4-1 ) = Z( J4 )*TEMP
|
|
D = D*TEMP
|
|
ELSE
|
|
Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
|
|
D = Z( J4+2 )*( D / Z( J4-3 ) )
|
|
END IF
|
|
DMIN = MIN( DMIN, D )
|
|
EMIN = MIN( EMIN, Z( J4-1 ) )
|
|
20 CONTINUE
|
|
END IF
|
|
*
|
|
* Unroll last two steps.
|
|
*
|
|
DNM2 = D
|
|
DMIN2 = DMIN
|
|
J4 = 4*( N0-2 ) - PP
|
|
J4P2 = J4 + 2*PP - 1
|
|
Z( J4-2 ) = DNM2 + Z( J4P2 )
|
|
IF( Z( J4-2 ).EQ.ZERO ) THEN
|
|
Z( J4 ) = ZERO
|
|
DNM1 = Z( J4P2+2 )
|
|
DMIN = DNM1
|
|
EMIN = ZERO
|
|
ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
|
|
$ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
|
|
TEMP = Z( J4P2+2 ) / Z( J4-2 )
|
|
Z( J4 ) = Z( J4P2 )*TEMP
|
|
DNM1 = DNM2*TEMP
|
|
ELSE
|
|
Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
|
|
DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) )
|
|
END IF
|
|
DMIN = MIN( DMIN, DNM1 )
|
|
*
|
|
DMIN1 = DMIN
|
|
J4 = J4 + 4
|
|
J4P2 = J4 + 2*PP - 1
|
|
Z( J4-2 ) = DNM1 + Z( J4P2 )
|
|
IF( Z( J4-2 ).EQ.ZERO ) THEN
|
|
Z( J4 ) = ZERO
|
|
DN = Z( J4P2+2 )
|
|
DMIN = DN
|
|
EMIN = ZERO
|
|
ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND.
|
|
$ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN
|
|
TEMP = Z( J4P2+2 ) / Z( J4-2 )
|
|
Z( J4 ) = Z( J4P2 )*TEMP
|
|
DN = DNM1*TEMP
|
|
ELSE
|
|
Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
|
|
DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) )
|
|
END IF
|
|
DMIN = MIN( DMIN, DN )
|
|
*
|
|
Z( J4+2 ) = DN
|
|
Z( 4*N0-PP ) = EMIN
|
|
RETURN
|
|
*
|
|
* End of SLASQ6
|
|
*
|
|
END
|