new file mode 100644
--- /dev/null
+++ b/libcruft/blas/dmach.f
@@ -0,0 +1,58 @@
+      DOUBLE PRECISION FUNCTION DMACH(JOB)
+      INTEGER JOB
+C
+C     SMACH COMPUTES MACHINE PARAMETERS OF FLOATING POINT
+C     ARITHMETIC FOR USE IN TESTING ONLY.  NOT REQUIRED BY
+C     LINPACK PROPER.
+C
+C     IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES,
+C     THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS.
+C     ASSUME THE COMPUTER HAS
+C
+C        B = BASE OF ARITHMETIC
+C        T = NUMBER OF BASE  B  DIGITS
+C        L = SMALLEST POSSIBLE EXPONENT
+C        U = LARGEST POSSIBLE EXPONENT
+C
+C     THEN
+C
+C        EPS = B**(1-T)
+C        TINY = 100.0*B**(-L+T)
+C        HUGE = 0.01*B**(U-T)
+C
+C     DMACH SAME AS SMACH EXCEPT T, L, U APPLY TO
+C     DOUBLE PRECISION.
+C
+C     CMACH SAME AS SMACH EXCEPT IF COMPLEX DIVISION
+C     IS DONE BY
+C
+C        1/(X+I*Y) = (X-I*Y)/(X**2+Y**2)
+C
+C     THEN
+C
+C        TINY = SQRT(TINY)
+C        HUGE = SQRT(HUGE)
+C
+C
+C     JOB IS 1, 2 OR 3 FOR EPSILON, TINY AND HUGE, RESPECTIVELY.
+C
+      DOUBLE PRECISION EPS,TINY,HUGE,S
+C
+      EPS = 1.0D0
+   10 EPS = EPS/2.0D0
+      S = 1.0D0 + EPS
+      IF (S .GT. 1.0D0) GO TO 10
+      EPS = 2.0D0*EPS
+C
+      S = 1.0D0
+   20 TINY = S
+      S = S/16.0D0
+      IF (S*1.0 .NE. 0.0D0) GO TO 20
+      TINY = (TINY/EPS)*100.0
+      HUGE = 1.0D0/TINY
+C
+      IF (JOB .EQ. 1) DMACH = EPS
+      IF (JOB .EQ. 2) DMACH = TINY
+      IF (JOB .EQ. 3) DMACH = HUGE
+      RETURN
+      END