sync/ports
1
0
Fork 0
mirror of https://git.freebsd.org/ports.git synced 2025-06-15 17:50:31 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions
ports/math/lapack/files/patch-r1139
Maho Nakata 795bd78d81 Update with many small following changes. 
1. adding upstream patches and bump port revision.
2. adding precompiled man by default rather than building.
3. adding regression-test target, several tests fails, though.
4. prepare for xlapack and now blas is secondary port of this port.

Submitted by:	bf@
2012-01-22 12:02:17 +00:00

394 lines
12 KiB
Text
Raw Blame History

Index: SRC/slasq5.f
===================================================================
--- SRC/slasq5.f (revision 1138)
+++ SRC/slasq5.f (revision 1139)
@@ -129,8 +129,8 @@
*> \ingroup auxOTHERcomputational
*
* =====================================================================
- SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
- $ DNM1, DNM2, IEEE )
+ SUBROUTINE SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2,
+ $ DN, DNM1, DNM2, IEEE, EPS )
*
* -- LAPACK computational routine (version 3.4.0) --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -140,7 +140,8 @@
* .. Scalar Arguments ..
LOGICAL IEEE
INTEGER I0, N0, PP
- REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU
+ REAL DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU,
+ $ SIGMA, EPS
* ..
* .. Array Arguments ..
REAL Z( * )
@@ -149,12 +150,12 @@
* =====================================================================
*
* .. Parameter ..
- REAL ZERO
- PARAMETER ( ZERO = 0.0E0 )
+ REAL ZERO, HALF
+ PARAMETER ( ZERO = 0.0E0, HALF = 0.5 )
* ..
* .. Local Scalars ..
INTEGER J4, J4P2
- REAL D, EMIN, TEMP
+ REAL D, EMIN, TEMP, DTHRESH
* ..
* .. Intrinsic Functions ..
INTRINSIC MIN
@@ -164,114 +165,231 @@
IF( ( N0-I0-1 ).LE.0 )
$ RETURN
*
- J4 = 4*I0 + PP - 3
- EMIN = Z( J4+4 )
- D = Z( J4 ) - TAU
- DMIN = D
- DMIN1 = -Z( J4 )
-*
- IF( IEEE ) THEN
-*
-* Code for IEEE arithmetic.
-*
- IF( PP.EQ.0 ) THEN
- DO 10 J4 = 4*I0, 4*( N0-3 ), 4
- Z( J4-2 ) = D + Z( J4-1 )
- TEMP = Z( J4+1 ) / Z( J4-2 )
- D = D*TEMP - TAU
- DMIN = MIN( DMIN, D )
- Z( J4 ) = Z( J4-1 )*TEMP
- EMIN = MIN( Z( J4 ), EMIN )
- 10 CONTINUE
+ DTHRESH = EPS*(SIGMA+TAU)
+ IF( TAU.LT.DTHRESH*HALF ) TAU = ZERO
+ IF( TAU.NE.ZERO ) THEN
+ J4 = 4*I0 + PP - 3
+ EMIN = Z( J4+4 )
+ D = Z( J4 ) - TAU
+ DMIN = D
+ DMIN1 = -Z( J4 )
+*
+ IF( IEEE ) THEN
+*
+* Code for IEEE arithmetic.
+*
+ IF( PP.EQ.0 ) THEN
+ DO 10 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-2 ) = D + Z( J4-1 )
+ TEMP = Z( J4+1 ) / Z( J4-2 )
+ D = D*TEMP - TAU
+ DMIN = MIN( DMIN, D )
+ Z( J4 ) = Z( J4-1 )*TEMP
+ EMIN = MIN( Z( J4 ), EMIN )
+ 10 CONTINUE
+ ELSE
+ DO 20 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-3 ) = D + Z( J4 )
+ TEMP = Z( J4+2 ) / Z( J4-3 )
+ D = D*TEMP - TAU
+ DMIN = MIN( DMIN, D )
+ Z( J4-1 ) = Z( J4 )*TEMP
+ EMIN = MIN( Z( J4-1 ), EMIN )
+ 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 )
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
+ DMIN = MIN( DMIN, DNM1 )
+*
+ DMIN1 = DMIN
+ J4 = J4 + 4
+ J4P2 = J4 + 2*PP - 1
+ Z( J4-2 ) = DNM1 + Z( J4P2 )
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+ DMIN = MIN( DMIN, DN )
+*
ELSE
- DO 20 J4 = 4*I0, 4*( N0-3 ), 4
- Z( J4-3 ) = D + Z( J4 )
- TEMP = Z( J4+2 ) / Z( J4-3 )
- D = D*TEMP - TAU
- DMIN = MIN( DMIN, D )
- Z( J4-1 ) = Z( J4 )*TEMP
- EMIN = MIN( Z( J4-1 ), EMIN )
- 20 CONTINUE
+*
+* Code for non IEEE arithmetic.
+*
+ IF( PP.EQ.0 ) THEN
+ DO 30 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-2 ) = D + Z( J4-1 )
+ IF( D.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
+ D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, D )
+ EMIN = MIN( EMIN, Z( J4 ) )
+ 30 CONTINUE
+ ELSE
+ DO 40 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-3 ) = D + Z( J4 )
+ IF( D.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
+ D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, D )
+ EMIN = MIN( EMIN, Z( J4-1 ) )
+ 40 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( DNM2.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, DNM1 )
+*
+ DMIN1 = DMIN
+ J4 = J4 + 4
+ J4P2 = J4 + 2*PP - 1
+ Z( J4-2 ) = DNM1 + Z( J4P2 )
+ IF( DNM1.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, DN )
+*
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 )
- Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
- DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
- DMIN = MIN( DMIN, DNM1 )
-*
- DMIN1 = DMIN
- J4 = J4 + 4
- J4P2 = J4 + 2*PP - 1
- Z( J4-2 ) = DNM1 + Z( J4P2 )
- Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
- DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
- DMIN = MIN( DMIN, DN )
-*
ELSE
-*
-* Code for non IEEE arithmetic.
-*
- IF( PP.EQ.0 ) THEN
- DO 30 J4 = 4*I0, 4*( N0-3 ), 4
- Z( J4-2 ) = D + Z( J4-1 )
- IF( D.LT.ZERO ) THEN
- RETURN
- ELSE
- Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
- D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
- END IF
- DMIN = MIN( DMIN, D )
- EMIN = MIN( EMIN, Z( J4 ) )
- 30 CONTINUE
- ELSE
- DO 40 J4 = 4*I0, 4*( N0-3 ), 4
- Z( J4-3 ) = D + Z( J4 )
- IF( D.LT.ZERO ) THEN
- RETURN
- ELSE
- Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
- D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
- END IF
- DMIN = MIN( DMIN, D )
- EMIN = MIN( EMIN, Z( J4-1 ) )
- 40 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( DNM2.LT.ZERO ) THEN
- RETURN
- ELSE
+* This is the version that sets d's to zero if they are small enough
+ J4 = 4*I0 + PP - 3
+ EMIN = Z( J4+4 )
+ D = Z( J4 ) - TAU
+ DMIN = D
+ DMIN1 = -Z( J4 )
+ IF( IEEE ) THEN
+*
+* Code for IEEE arithmetic.
+*
+ IF( PP.EQ.0 ) THEN
+ DO 50 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-2 ) = D + Z( J4-1 )
+ TEMP = Z( J4+1 ) / Z( J4-2 )
+ D = D*TEMP - TAU
+ IF( D.LT.DTHRESH ) D = ZERO
+ DMIN = MIN( DMIN, D )
+ Z( J4 ) = Z( J4-1 )*TEMP
+ EMIN = MIN( Z( J4 ), EMIN )
+ 50 CONTINUE
+ ELSE
+ DO 60 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-3 ) = D + Z( J4 )
+ TEMP = Z( J4+2 ) / Z( J4-3 )
+ D = D*TEMP - TAU
+ IF( D.LT.DTHRESH ) D = ZERO
+ DMIN = MIN( DMIN, D )
+ Z( J4-1 ) = Z( J4 )*TEMP
+ EMIN = MIN( Z( J4-1 ), EMIN )
+ 60 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 )
Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
- END IF
- DMIN = MIN( DMIN, DNM1 )
-*
- DMIN1 = DMIN
- J4 = J4 + 4
- J4P2 = J4 + 2*PP - 1
- Z( J4-2 ) = DNM1 + Z( J4P2 )
- IF( DNM1.LT.ZERO ) THEN
- RETURN
- ELSE
+ DMIN = MIN( DMIN, DNM1 )
+*
+ DMIN1 = DMIN
+ J4 = J4 + 4
+ J4P2 = J4 + 2*PP - 1
+ Z( J4-2 ) = DNM1 + Z( J4P2 )
Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+ DMIN = MIN( DMIN, DN )
+*
+ ELSE
+*
+* Code for non IEEE arithmetic.
+*
+ IF( PP.EQ.0 ) THEN
+ DO 70 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-2 ) = D + Z( J4-1 )
+ IF( D.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) )
+ D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU
+ END IF
+ IF( D.LT.DTHRESH ) D = ZERO
+ DMIN = MIN( DMIN, D )
+ EMIN = MIN( EMIN, Z( J4 ) )
+ 70 CONTINUE
+ ELSE
+ DO 80 J4 = 4*I0, 4*( N0-3 ), 4
+ Z( J4-3 ) = D + Z( J4 )
+ IF( D.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) )
+ D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU
+ END IF
+ IF( D.LT.DTHRESH ) D = ZERO
+ DMIN = MIN( DMIN, D )
+ EMIN = MIN( EMIN, Z( J4-1 ) )
+ 80 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( DNM2.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, DNM1 )
+*
+ DMIN1 = DMIN
+ J4 = J4 + 4
+ J4P2 = J4 + 2*PP - 1
+ Z( J4-2 ) = DNM1 + Z( J4P2 )
+ IF( DNM1.LT.ZERO ) THEN
+ RETURN
+ ELSE
+ Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) )
+ DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU
+ END IF
+ DMIN = MIN( DMIN, DN )
+*
END IF
- DMIN = MIN( DMIN, DN )
-*
+*
END IF
-*
Z( J4+2 ) = DN
Z( 4*N0-PP ) = EMIN
RETURN
Index: SRC/slasq3.f
===================================================================
--- SRC/slasq3.f (revision 1138)
+++ SRC/slasq3.f (revision 1139)
@@ -331,15 +331,15 @@
*
70 CONTINUE
*
- CALL SLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN,
- $ DN1, DN2, IEEE )
+ CALL SLASQ5( I0, N0, Z, PP, TAU, SIGMA, DMIN, DMIN1, DMIN2, DN,
+ $ DN1, DN2, IEEE, EPS )
*
NDIV = NDIV + ( N0-I0+2 )
ITER = ITER + 1
*
* Check status.
*
- IF( DMIN.GE.ZERO .AND. DMIN1.GT.ZERO ) THEN
+ IF( DMIN.GE.ZERO .AND. DMIN1.GE.ZERO ) THEN
*
* Success.
*
Reference in a new issue View git blame Copy permalink