Time has come to have the fortran90 slamch and dlamch in the lapack package.

Jason (Riedy) wrote our ( ... his? :) ) ideas about it three years ago:
http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
Piotr (Luszczek) has written two subroutines, tested them on few platforms, collected the result on his webpage, and sent emails to the lapackers a few times. See:
http://www.cs.utk.edu/~luszczek/lapack/lamch.html
Theses slamch.f and dlamch.f subroutines were taken from PLASMA-2.1.0. 

Change to the LAPACK library:
* move the current LAPACK subroutine slamch.f (resp dlamch.f) as
 slamchf77.f (resp. dlamchf77.f),
* take the new slamch.f subroutines (resp. dlamch.f), remove the PLASMA
 header, have a LAPACK header, and insert the new routines in the
 library.

Minor:
* I would leave these routines compiled with the NOOPT flag.

Problem:
* CLAPACK: no idea how CLAPACK's going to handle this. CLAPACK can rely on
 IEEE arithmetic, can relay on float.h, or can rely on the previous
 xlamch.f

4 files changed, 1777 insertions, 1465 deletions

diff --git a/INSTALL/dlamch.f b/INSTALL/dlamch.f
index 8f830e87..aa90112b 100644
--- a/INSTALL/dlamch.f
+++ b/INSTALL/dlamch.f
@@ -1,8 +1,12 @@
       DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
 *
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
+*  -- LAPACK auxiliary routine (version 3.2.2) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     Based on LAPACK DLAMCH but with Fortran 95 query functions
+*     See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
+*     and  http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
+*     July 2010
 *
 *     .. Scalar Arguments ..
       CHARACTER          CMACH
@@ -49,526 +53,68 @@
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 *     ..
 *     .. Local Scalars ..
-      LOGICAL            FIRST, LRND
-      INTEGER            BETA, IMAX, IMIN, IT
-      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
-     $                   RND, SFMIN, SMALL, T
+      DOUBLE PRECISION   RND, EPS, SFMIN, SMALL, RMACH
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     ..
-*     .. External Subroutines ..
-      EXTERNAL           DLAMC2
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
-     $                   EMAX, RMAX, PREC
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. /
+*     .. Intrinsic Functions ..
+      INTRINSIC          DIGITS, EPSILON, HUGE, MAXEXPONENT,
+     $                   MINEXPONENT, RADIX, TINY
 *     ..
 *     .. Executable Statements ..
 *
-      IF( FIRST ) THEN
-         CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
-         BASE = BETA
-         T = IT
-         IF( LRND ) THEN
-            RND = ONE
-            EPS = ( BASE**( 1-IT ) ) / 2
-         ELSE
-            RND = ZERO
-            EPS = BASE**( 1-IT )
-         END IF
-         PREC = EPS*BASE
-         EMIN = IMIN
-         EMAX = IMAX
-         SFMIN = RMIN
-         SMALL = ONE / RMAX
-         IF( SMALL.GE.SFMIN ) THEN
 *
-*           Use SMALL plus a bit, to avoid the possibility of rounding
-*           causing overflow when computing  1/sfmin.
+*     Assume rounding, not chopping. Always.
 *
-            SFMIN = SMALL*( ONE+EPS )
-         END IF
+      RND = ONE
+*
+      IF( ONE.EQ.RND ) THEN
+         EPS = EPSILON(ZERO) * 0.5
+      ELSE
+         EPS = EPSILON(ZERO)
       END IF
 *
       IF( LSAME( CMACH, 'E' ) ) THEN
          RMACH = EPS
       ELSE IF( LSAME( CMACH, 'S' ) ) THEN
+         SFMIN = TINY(ZERO)
+         SMALL = ONE / HUGE(ZERO)
+         IF( SMALL.GE.SFMIN ) THEN
+*
+*           Use SMALL plus a bit, to avoid the possibility of rounding
+*           causing overflow when computing  1/sfmin.
+*
+            SFMIN = SMALL*( ONE+EPS )
+         END IF
          RMACH = SFMIN
       ELSE IF( LSAME( CMACH, 'B' ) ) THEN
-         RMACH = BASE
+         RMACH = RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'P' ) ) THEN
-         RMACH = PREC
+         RMACH = EPS * RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'N' ) ) THEN
-         RMACH = T
+         RMACH = DIGITS(ZERO)
       ELSE IF( LSAME( CMACH, 'R' ) ) THEN
          RMACH = RND
       ELSE IF( LSAME( CMACH, 'M' ) ) THEN
-         RMACH = EMIN
+         RMACH = MINEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'U' ) ) THEN
-         RMACH = RMIN
+         RMACH = tiny(zero)
       ELSE IF( LSAME( CMACH, 'L' ) ) THEN
-         RMACH = EMAX
+         RMACH = MAXEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'O' ) ) THEN
-         RMACH = RMAX
+         RMACH = HUGE(ZERO)
+      ELSE
+         RMACH = ZERO
       END IF
 *
       DLAMCH = RMACH
-      FIRST  = .FALSE.
       RETURN
 *
 *     End of DLAMCH
 *
       END
-*
-************************************************************************
-*
-      SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            IEEE1, RND
-      INTEGER            BETA, T
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
-*  IEEE1.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (output) INTEGER
-*          The base of the machine.
-*
-*  T       (output) INTEGER
-*          The number of ( BETA ) digits in the mantissa.
-*
-*  RND     (output) LOGICAL
-*          Specifies whether proper rounding  ( RND = .TRUE. )  or
-*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
-*          be a reliable guide to the way in which the machine performs
-*          its arithmetic.
-*
-*  IEEE1   (output) LOGICAL
-*          Specifies whether rounding appears to be done in the IEEE
-*          'round to nearest' style.
-*
-*  Further Details
-*  ===============
-*
-*  The routine is based on the routine  ENVRON  by Malcolm and
-*  incorporates suggestions by Gentleman and Marovich. See
-*
-*     Malcolm M. A. (1972) Algorithms to reveal properties of
-*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
-*
-*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
-*        that reveal properties of floating point arithmetic units.
-*        Comms. of the ACM, 17, 276-277.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      LOGICAL            FIRST, LIEEE1, LRND
-      INTEGER            LBETA, LT
-      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMC3
-      EXTERNAL           DLAMC3
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. /
-*     ..
-*     .. Executable Statements ..
-*
-      IF( FIRST ) THEN
-         ONE = 1
-*
-*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
-*        IEEE1, T and RND.
-*
-*        Throughout this routine  we use the function  DLAMC3  to ensure
-*        that relevant values are  stored and not held in registers,  or
-*        are not affected by optimizers.
-*
-*        Compute  a = 2.0**m  with the  smallest positive integer m such
-*        that
-*
-*           fl( a + 1.0 ) = a.
-*
-         A = 1
-         C = 1
-*
-*+       WHILE( C.EQ.ONE )LOOP
-   10    CONTINUE
-         IF( C.EQ.ONE ) THEN
-            A = 2*A
-            C = DLAMC3( A, ONE )
-            C = DLAMC3( C, -A )
-            GO TO 10
-         END IF
-*+       END WHILE
-*
-*        Now compute  b = 2.0**m  with the smallest positive integer m
-*        such that
-*
-*           fl( a + b ) .gt. a.
-*
-         B = 1
-         C = DLAMC3( A, B )
-*
-*+       WHILE( C.EQ.A )LOOP
-   20    CONTINUE
-         IF( C.EQ.A ) THEN
-            B = 2*B
-            C = DLAMC3( A, B )
-            GO TO 20
-         END IF
-*+       END WHILE
-*
-*        Now compute the base.  a and c  are neighbouring floating point
-*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
-*        their difference is beta. Adding 0.25 to c is to ensure that it
-*        is truncated to beta and not ( beta - 1 ).
-*
-         QTR = ONE / 4
-         SAVEC = C
-         C = DLAMC3( C, -A )
-         LBETA = C + QTR
-*
-*        Now determine whether rounding or chopping occurs,  by adding a
-*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
-*
-         B = LBETA
-         F = DLAMC3( B / 2, -B / 100 )
-         C = DLAMC3( F, A )
-         IF( C.EQ.A ) THEN
-            LRND = .TRUE.
-         ELSE
-            LRND = .FALSE.
-         END IF
-         F = DLAMC3( B / 2, B / 100 )
-         C = DLAMC3( F, A )
-         IF( ( LRND ) .AND. ( C.EQ.A ) )
-     $      LRND = .FALSE.
-*
-*        Try and decide whether rounding is done in the  IEEE  'round to
-*        nearest' style. B/2 is half a unit in the last place of the two
-*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
-*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
-*        A, but adding B/2 to SAVEC should change SAVEC.
-*
-         T1 = DLAMC3( B / 2, A )
-         T2 = DLAMC3( B / 2, SAVEC )
-         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
-*
-*        Now find  the  mantissa, t.  It should  be the  integer part of
-*        log to the base beta of a,  however it is safer to determine  t
-*        by powering.  So we find t as the smallest positive integer for
-*        which
-*
-*           fl( beta**t + 1.0 ) = 1.0.
-*
-         LT = 0
-         A = 1
-         C = 1
-*
-*+       WHILE( C.EQ.ONE )LOOP
-   30    CONTINUE
-         IF( C.EQ.ONE ) THEN
-            LT = LT + 1
-            A = A*LBETA
-            C = DLAMC3( A, ONE )
-            C = DLAMC3( C, -A )
-            GO TO 30
-         END IF
-*+       END WHILE
-*
-      END IF
-*
-      BETA = LBETA
-      T = LT
-      RND = LRND
-      IEEE1 = LIEEE1
-      FIRST = .FALSE.
-      RETURN
-*
-*     End of DLAMC1
-*
-      END
-*
-************************************************************************
-*
-      SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            RND
-      INTEGER            BETA, EMAX, EMIN, T
-      DOUBLE PRECISION   EPS, RMAX, RMIN
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DLAMC2 determines the machine parameters specified in its argument
-*  list.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (output) INTEGER
-*          The base of the machine.
-*
-*  T       (output) INTEGER
-*          The number of ( BETA ) digits in the mantissa.
-*
-*  RND     (output) LOGICAL
-*          Specifies whether proper rounding  ( RND = .TRUE. )  or
-*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
-*          be a reliable guide to the way in which the machine performs
-*          its arithmetic.
-*
-*  EPS     (output) DOUBLE PRECISION
-*          The smallest positive number such that
-*
-*             fl( 1.0 - EPS ) .LT. 1.0,
-*
-*          where fl denotes the computed value.
-*
-*  EMIN    (output) INTEGER
-*          The minimum exponent before (gradual) underflow occurs.
-*
-*  RMIN    (output) DOUBLE PRECISION
-*          The smallest normalized number for the machine, given by
-*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
-*          of BETA.
-*
-*  EMAX    (output) INTEGER
-*          The maximum exponent before overflow occurs.
-*
-*  RMAX    (output) DOUBLE PRECISION
-*          The largest positive number for the machine, given by
-*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
-*          value of BETA.
-*
-*  Further Details
-*  ===============
-*
-*  The computation of  EPS  is based on a routine PARANOIA by
-*  W. Kahan of the University of California at Berkeley.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
-      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
-     $                   NGNMIN, NGPMIN
-      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
-     $                   SIXTH, SMALL, THIRD, TWO, ZERO
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMC3
-      EXTERNAL           DLAMC3
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           DLAMC1, DLAMC4, DLAMC5
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, MAX, MIN
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
-     $                   LRMIN, LT
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
-*     ..
-*     .. Executable Statements ..
-*
-      IF( FIRST ) THEN
-         ZERO = 0
-         ONE = 1
-         TWO = 2
-*
-*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
-*        BETA, T, RND, EPS, EMIN and RMIN.
-*
-*        Throughout this routine  we use the function  DLAMC3  to ensure
-*        that relevant values are stored  and not held in registers,  or
-*        are not affected by optimizers.
-*
-*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
-*
-         CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
-*
-*        Start to find EPS.
-*
-         B = LBETA
-         A = B**( -LT )
-         LEPS = A
-*
-*        Try some tricks to see whether or not this is the correct  EPS.
-*
-         B = TWO / 3
-         HALF = ONE / 2
-         SIXTH = DLAMC3( B, -HALF )
-         THIRD = DLAMC3( SIXTH, SIXTH )
-         B = DLAMC3( THIRD, -HALF )
-         B = DLAMC3( B, SIXTH )
-         B = ABS( B )
-         IF( B.LT.LEPS )
-     $      B = LEPS
-*
-         LEPS = 1
-*
-*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
-   10    CONTINUE
-         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
-            LEPS = B
-            C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
-            C = DLAMC3( HALF, -C )
-            B = DLAMC3( HALF, C )
-            C = DLAMC3( HALF, -B )
-            B = DLAMC3( HALF, C )
-            GO TO 10
-         END IF
-*+       END WHILE
-*
-         IF( A.LT.LEPS )
-     $      LEPS = A
-*
-*        Computation of EPS complete.
-*
-*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
-*        Keep dividing  A by BETA until (gradual) underflow occurs. This
-*        is detected when we cannot recover the previous A.
-*
-         RBASE = ONE / LBETA
-         SMALL = ONE
-         DO 20 I = 1, 3
-            SMALL = DLAMC3( SMALL*RBASE, ZERO )
-   20    CONTINUE
-         A = DLAMC3( ONE, SMALL )
-         CALL DLAMC4( NGPMIN, ONE, LBETA )
-         CALL DLAMC4( NGNMIN, -ONE, LBETA )
-         CALL DLAMC4( GPMIN, A, LBETA )
-         CALL DLAMC4( GNMIN, -A, LBETA )
-         IEEE = .FALSE.
-*
-         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
-            IF( NGPMIN.EQ.GPMIN ) THEN
-               LEMIN = NGPMIN
-*            ( Non twos-complement machines, no gradual underflow;
-*              e.g.,  VAX )
-            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
-               LEMIN = NGPMIN - 1 + LT
-               IEEE = .TRUE.
-*            ( Non twos-complement machines, with gradual underflow;
-*              e.g., IEEE standard followers )
-            ELSE
-               LEMIN = MIN( NGPMIN, GPMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
-            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
-               LEMIN = MAX( NGPMIN, NGNMIN )
-*            ( Twos-complement machines, no gradual underflow;
-*              e.g., CYBER 205 )
-            ELSE
-               LEMIN = MIN( NGPMIN, NGNMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
-     $            ( GPMIN.EQ.GNMIN ) ) THEN
-            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
-               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
-*            ( Twos-complement machines with gradual underflow;
-*              no known machine )
-            ELSE
-               LEMIN = MIN( NGPMIN, NGNMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE
-            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
-*         ( A guess; no known machine )
-            IWARN = .TRUE.
-         END IF
-         FIRST = .FALSE.
-***
-* Comment out this if block if EMIN is ok
-         IF( IWARN ) THEN
-            FIRST = .TRUE.
-            WRITE( 6, FMT = 9999 )LEMIN
-         END IF
-***
-*
-*        Assume IEEE arithmetic if we found denormalised  numbers above,
-*        or if arithmetic seems to round in the  IEEE style,  determined
-*        in routine DLAMC1. A true IEEE machine should have both  things
-*        true; however, faulty machines may have one or the other.
-*
-         IEEE = IEEE .OR. LIEEE1
-*
-*        Compute  RMIN by successive division by  BETA. We could compute
-*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
-*        this computation.
-*
-         LRMIN = 1
-         DO 30 I = 1, 1 - LEMIN
-            LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
-   30    CONTINUE
-*
-*        Finally, call DLAMC5 to compute EMAX and RMAX.
-*
-         CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
-      END IF
-*
-      BETA = LBETA
-      T = LT
-      RND = LRND
-      EPS = LEPS
-      EMIN = LEMIN
-      RMIN = LRMIN
-      EMAX = LEMAX
-      RMAX = LRMAX
-*
-      RETURN
-*
- 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
-     $      '  EMIN = ', I8, /
-     $      ' If, after inspection, the value EMIN looks',
-     $      ' acceptable please comment out ',
-     $      / ' the IF block as marked within the code of routine',
-     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
-*
-*     End of DLAMC2
-*
-      END
-*
 ************************************************************************
 *
       DOUBLE PRECISION FUNCTION DLAMC3( A, B )
@@ -608,245 +154,3 @@
       END
 *
 ************************************************************************
-*
-      SUBROUTINE DLAMC4( EMIN, START, BASE )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      INTEGER            BASE, EMIN
-      DOUBLE PRECISION   START
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DLAMC4 is a service routine for DLAMC2.
-*
-*  Arguments
-*  =========
-*
-*  EMIN    (output) INTEGER 
-*          The minimum exponent before (gradual) underflow, computed by
-*          setting A = START and dividing by BASE until the previous A
-*          can not be recovered.
-*
-*  START   (input) DOUBLE PRECISION
-*          The starting point for determining EMIN.
-*
-*  BASE    (input) INTEGER
-*          The base of the machine.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      INTEGER            I
-      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMC3
-      EXTERNAL           DLAMC3
-*     ..
-*     .. Executable Statements ..
-*
-      A = START
-      ONE = 1
-      RBASE = ONE / BASE
-      ZERO = 0
-      EMIN = 1
-      B1 = DLAMC3( A*RBASE, ZERO )
-      C1 = A
-      C2 = A
-      D1 = A
-      D2 = A
-*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
-*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
-   10 CONTINUE
-      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
-     $    ( D2.EQ.A ) ) THEN
-         EMIN = EMIN - 1
-         A = B1
-         B1 = DLAMC3( A / BASE, ZERO )
-         C1 = DLAMC3( B1*BASE, ZERO )
-         D1 = ZERO
-         DO 20 I = 1, BASE
-            D1 = D1 + B1
-   20    CONTINUE
-         B2 = DLAMC3( A*RBASE, ZERO )
-         C2 = DLAMC3( B2 / RBASE, ZERO )
-         D2 = ZERO
-         DO 30 I = 1, BASE
-            D2 = D2 + B2
-   30    CONTINUE
-         GO TO 10
-      END IF
-*+    END WHILE
-*
-      RETURN
-*
-*     End of DLAMC4
-*
-      END
-*
-************************************************************************
-*
-      SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            IEEE
-      INTEGER            BETA, EMAX, EMIN, P
-      DOUBLE PRECISION   RMAX
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
-*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
-*  approximately to a power of 2.  It will fail on machines where this
-*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
-*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
-*  too large (i.e. too close to zero), probably with overflow.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (input) INTEGER
-*          The base of floating-point arithmetic.
-*
-*  P       (input) INTEGER
-*          The number of base BETA digits in the mantissa of a
-*          floating-point value.
-*
-*  EMIN    (input) INTEGER
-*          The minimum exponent before (gradual) underflow.
-*
-*  IEEE    (input) LOGICAL
-*          A logical flag specifying whether or not the arithmetic
-*          system is thought to comply with the IEEE standard.
-*
-*  EMAX    (output) INTEGER
-*          The largest exponent before overflow
-*
-*  RMAX    (output) DOUBLE PRECISION
-*          The largest machine floating-point number.
-*
-* =====================================================================
-*
-*     .. Parameters ..
-      DOUBLE PRECISION   ZERO, ONE
-      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
-*     ..
-*     .. Local Scalars ..
-      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
-      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
-*     ..
-*     .. External Functions ..
-      DOUBLE PRECISION   DLAMC3
-      EXTERNAL           DLAMC3
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          MOD
-*     ..
-*     .. Executable Statements ..
-*
-*     First compute LEXP and UEXP, two powers of 2 that bound
-*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
-*     approximately to the bound that is closest to abs(EMIN).
-*     (EMAX is the exponent of the required number RMAX).
-*
-      LEXP = 1
-      EXBITS = 1
-   10 CONTINUE
-      TRY = LEXP*2
-      IF( TRY.LE.( -EMIN ) ) THEN
-         LEXP = TRY
-         EXBITS = EXBITS + 1
-         GO TO 10
-      END IF
-      IF( LEXP.EQ.-EMIN ) THEN
-         UEXP = LEXP
-      ELSE
-         UEXP = TRY
-         EXBITS = EXBITS + 1
-      END IF
-*
-*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
-*     than or equal to EMIN. EXBITS is the number of bits needed to
-*     store the exponent.
-*
-      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
-         EXPSUM = 2*LEXP
-      ELSE
-         EXPSUM = 2*UEXP
-      END IF
-*
-*     EXPSUM is the exponent range, approximately equal to
-*     EMAX - EMIN + 1 .
-*
-      EMAX = EXPSUM + EMIN - 1
-      NBITS = 1 + EXBITS + P
-*
-*     NBITS is the total number of bits needed to store a
-*     floating-point number.
-*
-      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
-*
-*        Either there are an odd number of bits used to store a
-*        floating-point number, which is unlikely, or some bits are
-*        not used in the representation of numbers, which is possible,
-*        (e.g. Cray machines) or the mantissa has an implicit bit,
-*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
-*        most likely. We have to assume the last alternative.
-*        If this is true, then we need to reduce EMAX by one because
-*        there must be some way of representing zero in an implicit-bit
-*        system. On machines like Cray, we are reducing EMAX by one
-*        unnecessarily.
-*
-         EMAX = EMAX - 1
-      END IF
-*
-      IF( IEEE ) THEN
-*
-*        Assume we are on an IEEE machine which reserves one exponent
-*        for infinity and NaN.
-*
-         EMAX = EMAX - 1
-      END IF
-*
-*     Now create RMAX, the largest machine number, which should
-*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
-*
-*     First compute 1.0 - BETA**(-P), being careful that the
-*     result is less than 1.0 .
-*
-      RECBAS = ONE / BETA
-      Z = BETA - ONE
-      Y = ZERO
-      DO 20 I = 1, P
-         Z = Z*RECBAS
-         IF( Y.LT.ONE )
-     $      OLDY = Y
-         Y = DLAMC3( Y, Z )
-   20 CONTINUE
-      IF( Y.GE.ONE )
-     $   Y = OLDY
-*
-*     Now multiply by BETA**EMAX to get RMAX.
-*
-      DO 30 I = 1, EMAX
-         Y = DLAMC3( Y*BETA, ZERO )
-   30 CONTINUE
-*
-      RMAX = Y
-      RETURN
-*
-*     End of DLAMC5
-*
-      END
diff --git a/INSTALL/dlamchf77.f b/INSTALL/dlamchf77.f
new file mode 100644
index 00000000..8f830e87
--- /dev/null
+++ b/INSTALL/dlamchf77.f
@@ -0,0 +1,852 @@
+      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          CMACH
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMCH determines double precision machine parameters.
+*
+*  Arguments
+*  =========
+*
+*  CMACH   (input) CHARACTER*1
+*          Specifies the value to be returned by DLAMCH:
+*          = 'E' or 'e',   DLAMCH := eps
+*          = 'S' or 's ,   DLAMCH := sfmin
+*          = 'B' or 'b',   DLAMCH := base
+*          = 'P' or 'p',   DLAMCH := eps*base
+*          = 'N' or 'n',   DLAMCH := t
+*          = 'R' or 'r',   DLAMCH := rnd
+*          = 'M' or 'm',   DLAMCH := emin
+*          = 'U' or 'u',   DLAMCH := rmin
+*          = 'L' or 'l',   DLAMCH := emax
+*          = 'O' or 'o',   DLAMCH := rmax
+*
+*          where
+*
+*          eps   = relative machine precision
+*          sfmin = safe minimum, such that 1/sfmin does not overflow
+*          base  = base of the machine
+*          prec  = eps*base
+*          t     = number of (base) digits in the mantissa
+*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
+*          emin  = minimum exponent before (gradual) underflow
+*          rmin  = underflow threshold - base**(emin-1)
+*          emax  = largest exponent before overflow
+*          rmax  = overflow threshold  - (base**emax)*(1-eps)
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LRND
+      INTEGER            BETA, IMAX, IMIN, IT
+      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
+     $                   RND, SFMIN, SMALL, T
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLAMC2
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
+     $                   EMAX, RMAX, PREC
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
+         BASE = BETA
+         T = IT
+         IF( LRND ) THEN
+            RND = ONE
+            EPS = ( BASE**( 1-IT ) ) / 2
+         ELSE
+            RND = ZERO
+            EPS = BASE**( 1-IT )
+         END IF
+         PREC = EPS*BASE
+         EMIN = IMIN
+         EMAX = IMAX
+         SFMIN = RMIN
+         SMALL = ONE / RMAX
+         IF( SMALL.GE.SFMIN ) THEN
+*
+*           Use SMALL plus a bit, to avoid the possibility of rounding
+*           causing overflow when computing  1/sfmin.
+*
+            SFMIN = SMALL*( ONE+EPS )
+         END IF
+      END IF
+*
+      IF( LSAME( CMACH, 'E' ) ) THEN
+         RMACH = EPS
+      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
+         RMACH = SFMIN
+      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
+         RMACH = BASE
+      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
+         RMACH = PREC
+      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
+         RMACH = T
+      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
+         RMACH = RND
+      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
+         RMACH = EMIN
+      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
+         RMACH = RMIN
+      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
+         RMACH = EMAX
+      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
+         RMACH = RMAX
+      END IF
+*
+      DLAMCH = RMACH
+      FIRST  = .FALSE.
+      RETURN
+*
+*     End of DLAMCH
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE1, RND
+      INTEGER            BETA, T
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
+*  IEEE1.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  IEEE1   (output) LOGICAL
+*          Specifies whether rounding appears to be done in the IEEE
+*          'round to nearest' style.
+*
+*  Further Details
+*  ===============
+*
+*  The routine is based on the routine  ENVRON  by Malcolm and
+*  incorporates suggestions by Gentleman and Marovich. See
+*
+*     Malcolm M. A. (1972) Algorithms to reveal properties of
+*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
+*
+*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
+*        that reveal properties of floating point arithmetic units.
+*        Comms. of the ACM, 17, 276-277.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LIEEE1, LRND
+      INTEGER            LBETA, LT
+      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMC3
+      EXTERNAL           DLAMC3
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         ONE = 1
+*
+*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
+*        IEEE1, T and RND.
+*
+*        Throughout this routine  we use the function  DLAMC3  to ensure
+*        that relevant values are  stored and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        Compute  a = 2.0**m  with the  smallest positive integer m such
+*        that
+*
+*           fl( a + 1.0 ) = a.
+*
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   10    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            A = 2*A
+            C = DLAMC3( A, ONE )
+            C = DLAMC3( C, -A )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+*        Now compute  b = 2.0**m  with the smallest positive integer m
+*        such that
+*
+*           fl( a + b ) .gt. a.
+*
+         B = 1
+         C = DLAMC3( A, B )
+*
+*+       WHILE( C.EQ.A )LOOP
+   20    CONTINUE
+         IF( C.EQ.A ) THEN
+            B = 2*B
+            C = DLAMC3( A, B )
+            GO TO 20
+         END IF
+*+       END WHILE
+*
+*        Now compute the base.  a and c  are neighbouring floating point
+*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
+*        their difference is beta. Adding 0.25 to c is to ensure that it
+*        is truncated to beta and not ( beta - 1 ).
+*
+         QTR = ONE / 4
+         SAVEC = C
+         C = DLAMC3( C, -A )
+         LBETA = C + QTR
+*
+*        Now determine whether rounding or chopping occurs,  by adding a
+*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
+*
+         B = LBETA
+         F = DLAMC3( B / 2, -B / 100 )
+         C = DLAMC3( F, A )
+         IF( C.EQ.A ) THEN
+            LRND = .TRUE.
+         ELSE
+            LRND = .FALSE.
+         END IF
+         F = DLAMC3( B / 2, B / 100 )
+         C = DLAMC3( F, A )
+         IF( ( LRND ) .AND. ( C.EQ.A ) )
+     $      LRND = .FALSE.
+*
+*        Try and decide whether rounding is done in the  IEEE  'round to
+*        nearest' style. B/2 is half a unit in the last place of the two
+*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
+*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
+*        A, but adding B/2 to SAVEC should change SAVEC.
+*
+         T1 = DLAMC3( B / 2, A )
+         T2 = DLAMC3( B / 2, SAVEC )
+         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
+*
+*        Now find  the  mantissa, t.  It should  be the  integer part of
+*        log to the base beta of a,  however it is safer to determine  t
+*        by powering.  So we find t as the smallest positive integer for
+*        which
+*
+*           fl( beta**t + 1.0 ) = 1.0.
+*
+         LT = 0
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   30    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            LT = LT + 1
+            A = A*LBETA
+            C = DLAMC3( A, ONE )
+            C = DLAMC3( C, -A )
+            GO TO 30
+         END IF
+*+       END WHILE
+*
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      IEEE1 = LIEEE1
+      FIRST = .FALSE.
+      RETURN
+*
+*     End of DLAMC1
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            RND
+      INTEGER            BETA, EMAX, EMIN, T
+      DOUBLE PRECISION   EPS, RMAX, RMIN
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC2 determines the machine parameters specified in its argument
+*  list.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  EPS     (output) DOUBLE PRECISION
+*          The smallest positive number such that
+*
+*             fl( 1.0 - EPS ) .LT. 1.0,
+*
+*          where fl denotes the computed value.
+*
+*  EMIN    (output) INTEGER
+*          The minimum exponent before (gradual) underflow occurs.
+*
+*  RMIN    (output) DOUBLE PRECISION
+*          The smallest normalized number for the machine, given by
+*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
+*          of BETA.
+*
+*  EMAX    (output) INTEGER
+*          The maximum exponent before overflow occurs.
+*
+*  RMAX    (output) DOUBLE PRECISION
+*          The largest positive number for the machine, given by
+*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
+*          value of BETA.
+*
+*  Further Details
+*  ===============
+*
+*  The computation of  EPS  is based on a routine PARANOIA by
+*  W. Kahan of the University of California at Berkeley.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
+      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
+     $                   NGNMIN, NGPMIN
+      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
+     $                   SIXTH, SMALL, THIRD, TWO, ZERO
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMC3
+      EXTERNAL           DLAMC3
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DLAMC1, DLAMC4, DLAMC5
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
+     $                   LRMIN, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         ZERO = 0
+         ONE = 1
+         TWO = 2
+*
+*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
+*        BETA, T, RND, EPS, EMIN and RMIN.
+*
+*        Throughout this routine  we use the function  DLAMC3  to ensure
+*        that relevant values are stored  and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
+*
+         CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
+*
+*        Start to find EPS.
+*
+         B = LBETA
+         A = B**( -LT )
+         LEPS = A
+*
+*        Try some tricks to see whether or not this is the correct  EPS.
+*
+         B = TWO / 3
+         HALF = ONE / 2
+         SIXTH = DLAMC3( B, -HALF )
+         THIRD = DLAMC3( SIXTH, SIXTH )
+         B = DLAMC3( THIRD, -HALF )
+         B = DLAMC3( B, SIXTH )
+         B = ABS( B )
+         IF( B.LT.LEPS )
+     $      B = LEPS
+*
+         LEPS = 1
+*
+*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
+   10    CONTINUE
+         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
+            LEPS = B
+            C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
+            C = DLAMC3( HALF, -C )
+            B = DLAMC3( HALF, C )
+            C = DLAMC3( HALF, -B )
+            B = DLAMC3( HALF, C )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+         IF( A.LT.LEPS )
+     $      LEPS = A
+*
+*        Computation of EPS complete.
+*
+*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
+*        Keep dividing  A by BETA until (gradual) underflow occurs. This
+*        is detected when we cannot recover the previous A.
+*
+         RBASE = ONE / LBETA
+         SMALL = ONE
+         DO 20 I = 1, 3
+            SMALL = DLAMC3( SMALL*RBASE, ZERO )
+   20    CONTINUE
+         A = DLAMC3( ONE, SMALL )
+         CALL DLAMC4( NGPMIN, ONE, LBETA )
+         CALL DLAMC4( NGNMIN, -ONE, LBETA )
+         CALL DLAMC4( GPMIN, A, LBETA )
+         CALL DLAMC4( GNMIN, -A, LBETA )
+         IEEE = .FALSE.
+*
+         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( NGPMIN.EQ.GPMIN ) THEN
+               LEMIN = NGPMIN
+*            ( Non twos-complement machines, no gradual underflow;
+*              e.g.,  VAX )
+            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
+               LEMIN = NGPMIN - 1 + LT
+               IEEE = .TRUE.
+*            ( Non twos-complement machines, with gradual underflow;
+*              e.g., IEEE standard followers )
+            ELSE
+               LEMIN = MIN( NGPMIN, GPMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
+            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN )
+*            ( Twos-complement machines, no gradual underflow;
+*              e.g., CYBER 205 )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
+     $            ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
+*            ( Twos-complement machines with gradual underflow;
+*              no known machine )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE
+            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
+*         ( A guess; no known machine )
+            IWARN = .TRUE.
+         END IF
+         FIRST = .FALSE.
+***
+* Comment out this if block if EMIN is ok
+         IF( IWARN ) THEN
+            FIRST = .TRUE.
+            WRITE( 6, FMT = 9999 )LEMIN
+         END IF
+***
+*
+*        Assume IEEE arithmetic if we found denormalised  numbers above,
+*        or if arithmetic seems to round in the  IEEE style,  determined
+*        in routine DLAMC1. A true IEEE machine should have both  things
+*        true; however, faulty machines may have one or the other.
+*
+         IEEE = IEEE .OR. LIEEE1
+*
+*        Compute  RMIN by successive division by  BETA. We could compute
+*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
+*        this computation.
+*
+         LRMIN = 1
+         DO 30 I = 1, 1 - LEMIN
+            LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
+   30    CONTINUE
+*
+*        Finally, call DLAMC5 to compute EMAX and RMAX.
+*
+         CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      EPS = LEPS
+      EMIN = LEMIN
+      RMIN = LRMIN
+      EMAX = LEMAX
+      RMAX = LRMAX
+*
+      RETURN
+*
+ 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
+     $      '  EMIN = ', I8, /
+     $      ' If, after inspection, the value EMIN looks',
+     $      ' acceptable please comment out ',
+     $      / ' the IF block as marked within the code of routine',
+     $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
+*
+*     End of DLAMC2
+*
+      END
+*
+************************************************************************
+*
+      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   A, B
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing
+*  the addition of  A  and  B ,  for use in situations where optimizers
+*  might hold one of these in a register.
+*
+*  Arguments
+*  =========
+*
+*  A       (input) DOUBLE PRECISION
+*  B       (input) DOUBLE PRECISION
+*          The values A and B.
+*
+* =====================================================================
+*
+*     .. Executable Statements ..
+*
+      DLAMC3 = A + B
+*
+      RETURN
+*
+*     End of DLAMC3
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE DLAMC4( EMIN, START, BASE )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            BASE, EMIN
+      DOUBLE PRECISION   START
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC4 is a service routine for DLAMC2.
+*
+*  Arguments
+*  =========
+*
+*  EMIN    (output) INTEGER 
+*          The minimum exponent before (gradual) underflow, computed by
+*          setting A = START and dividing by BASE until the previous A
+*          can not be recovered.
+*
+*  START   (input) DOUBLE PRECISION
+*          The starting point for determining EMIN.
+*
+*  BASE    (input) INTEGER
+*          The base of the machine.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER            I
+      DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMC3
+      EXTERNAL           DLAMC3
+*     ..
+*     .. Executable Statements ..
+*
+      A = START
+      ONE = 1
+      RBASE = ONE / BASE
+      ZERO = 0
+      EMIN = 1
+      B1 = DLAMC3( A*RBASE, ZERO )
+      C1 = A
+      C2 = A
+      D1 = A
+      D2 = A
+*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
+*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
+   10 CONTINUE
+      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
+     $    ( D2.EQ.A ) ) THEN
+         EMIN = EMIN - 1
+         A = B1
+         B1 = DLAMC3( A / BASE, ZERO )
+         C1 = DLAMC3( B1*BASE, ZERO )
+         D1 = ZERO
+         DO 20 I = 1, BASE
+            D1 = D1 + B1
+   20    CONTINUE
+         B2 = DLAMC3( A*RBASE, ZERO )
+         C2 = DLAMC3( B2 / RBASE, ZERO )
+         D2 = ZERO
+         DO 30 I = 1, BASE
+            D2 = D2 + B2
+   30    CONTINUE
+         GO TO 10
+      END IF
+*+    END WHILE
+*
+      RETURN
+*
+*     End of DLAMC4
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE
+      INTEGER            BETA, EMAX, EMIN, P
+      DOUBLE PRECISION   RMAX
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
+*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
+*  approximately to a power of 2.  It will fail on machines where this
+*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
+*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
+*  too large (i.e. too close to zero), probably with overflow.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (input) INTEGER
+*          The base of floating-point arithmetic.
+*
+*  P       (input) INTEGER
+*          The number of base BETA digits in the mantissa of a
+*          floating-point value.
+*
+*  EMIN    (input) INTEGER
+*          The minimum exponent before (gradual) underflow.
+*
+*  IEEE    (input) LOGICAL
+*          A logical flag specifying whether or not the arithmetic
+*          system is thought to comply with the IEEE standard.
+*
+*  EMAX    (output) INTEGER
+*          The largest exponent before overflow
+*
+*  RMAX    (output) DOUBLE PRECISION
+*          The largest machine floating-point number.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
+      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMC3
+      EXTERNAL           DLAMC3
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MOD
+*     ..
+*     .. Executable Statements ..
+*
+*     First compute LEXP and UEXP, two powers of 2 that bound
+*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
+*     approximately to the bound that is closest to abs(EMIN).
+*     (EMAX is the exponent of the required number RMAX).
+*
+      LEXP = 1
+      EXBITS = 1
+   10 CONTINUE
+      TRY = LEXP*2
+      IF( TRY.LE.( -EMIN ) ) THEN
+         LEXP = TRY
+         EXBITS = EXBITS + 1
+         GO TO 10
+      END IF
+      IF( LEXP.EQ.-EMIN ) THEN
+         UEXP = LEXP
+      ELSE
+         UEXP = TRY
+         EXBITS = EXBITS + 1
+      END IF
+*
+*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
+*     than or equal to EMIN. EXBITS is the number of bits needed to
+*     store the exponent.
+*
+      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
+         EXPSUM = 2*LEXP
+      ELSE
+         EXPSUM = 2*UEXP
+      END IF
+*
+*     EXPSUM is the exponent range, approximately equal to
+*     EMAX - EMIN + 1 .
+*
+      EMAX = EXPSUM + EMIN - 1
+      NBITS = 1 + EXBITS + P
+*
+*     NBITS is the total number of bits needed to store a
+*     floating-point number.
+*
+      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
+*
+*        Either there are an odd number of bits used to store a
+*        floating-point number, which is unlikely, or some bits are
+*        not used in the representation of numbers, which is possible,
+*        (e.g. Cray machines) or the mantissa has an implicit bit,
+*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
+*        most likely. We have to assume the last alternative.
+*        If this is true, then we need to reduce EMAX by one because
+*        there must be some way of representing zero in an implicit-bit
+*        system. On machines like Cray, we are reducing EMAX by one
+*        unnecessarily.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+      IF( IEEE ) THEN
+*
+*        Assume we are on an IEEE machine which reserves one exponent
+*        for infinity and NaN.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+*     Now create RMAX, the largest machine number, which should
+*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
+*
+*     First compute 1.0 - BETA**(-P), being careful that the
+*     result is less than 1.0 .
+*
+      RECBAS = ONE / BETA
+      Z = BETA - ONE
+      Y = ZERO
+      DO 20 I = 1, P
+         Z = Z*RECBAS
+         IF( Y.LT.ONE )
+     $      OLDY = Y
+         Y = DLAMC3( Y, Z )
+   20 CONTINUE
+      IF( Y.GE.ONE )
+     $   Y = OLDY
+*
+*     Now multiply by BETA**EMAX to get RMAX.
+*
+      DO 30 I = 1, EMAX
+         Y = DLAMC3( Y*BETA, ZERO )
+   30 CONTINUE
+*
+      RMAX = Y
+      RETURN
+*
+*     End of DLAMC5
+*
+      END
diff --git a/INSTALL/slamch.f b/INSTALL/slamch.f
index 5bb2c7be..902f3220 100644
--- a/INSTALL/slamch.f
+++ b/INSTALL/slamch.f
@@ -1,8 +1,12 @@
       REAL             FUNCTION SLAMCH( CMACH )
 *
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
+*  -- LAPACK auxiliary routine (version 3.2.2) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     Based on LAPACK DLAMCH but with Fortran 95 query functions
+*     See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
+*     and  http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
+*     July 2010
 *
 *     .. Scalar Arguments ..
       CHARACTER          CMACH
@@ -49,526 +53,68 @@
       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 *     ..
 *     .. Local Scalars ..
-      LOGICAL            FIRST, LRND
-      INTEGER            BETA, IMAX, IMIN, IT
-      REAL               BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
-     $                   RND, SFMIN, SMALL, T
+      REAL               RND, EPS, SFMIN, SMALL, RMACH
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
       EXTERNAL           LSAME
 *     ..
-*     .. External Subroutines ..
-      EXTERNAL           SLAMC2
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
-     $                   EMAX, RMAX, PREC
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. /
+*     .. Intrinsic Functions ..
+      INTRINSIC          DIGITS, EPSILON, HUGE, MAXEXPONENT,
+     $                   MINEXPONENT, RADIX, TINY
 *     ..
 *     .. Executable Statements ..
 *
-      IF( FIRST ) THEN
-         CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
-         BASE = BETA
-         T = IT
-         IF( LRND ) THEN
-            RND = ONE
-            EPS = ( BASE**( 1-IT ) ) / 2
-         ELSE
-            RND = ZERO
-            EPS = BASE**( 1-IT )
-         END IF
-         PREC = EPS*BASE
-         EMIN = IMIN
-         EMAX = IMAX
-         SFMIN = RMIN
-         SMALL = ONE / RMAX
-         IF( SMALL.GE.SFMIN ) THEN
 *
-*           Use SMALL plus a bit, to avoid the possibility of rounding
-*           causing overflow when computing  1/sfmin.
+*     Assume rounding, not chopping. Always.
 *
-            SFMIN = SMALL*( ONE+EPS )
-         END IF
+      RND = ONE
+*
+      IF( ONE.EQ.RND ) THEN
+         EPS = EPSILON(ZERO) * 0.5
+      ELSE
+         EPS = EPSILON(ZERO)
       END IF
 *
       IF( LSAME( CMACH, 'E' ) ) THEN
          RMACH = EPS
       ELSE IF( LSAME( CMACH, 'S' ) ) THEN
+         SFMIN = TINY(ZERO)
+         SMALL = ONE / HUGE(ZERO)
+         IF( SMALL.GE.SFMIN ) THEN
+*
+*           Use SMALL plus a bit, to avoid the possibility of rounding
+*           causing overflow when computing  1/sfmin.
+*
+            SFMIN = SMALL*( ONE+EPS )
+         END IF
          RMACH = SFMIN
       ELSE IF( LSAME( CMACH, 'B' ) ) THEN
-         RMACH = BASE
+         RMACH = RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'P' ) ) THEN
-         RMACH = PREC
+         RMACH = EPS * RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'N' ) ) THEN
-         RMACH = T
+         RMACH = DIGITS(ZERO)
       ELSE IF( LSAME( CMACH, 'R' ) ) THEN
          RMACH = RND
       ELSE IF( LSAME( CMACH, 'M' ) ) THEN
-         RMACH = EMIN
+         RMACH = MINEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'U' ) ) THEN
-         RMACH = RMIN
+         RMACH = tiny(zero)
       ELSE IF( LSAME( CMACH, 'L' ) ) THEN
-         RMACH = EMAX
+         RMACH = MAXEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'O' ) ) THEN
-         RMACH = RMAX
+         RMACH = HUGE(ZERO)
+      ELSE
+         RMACH = ZERO
       END IF
 *
       SLAMCH = RMACH
-      FIRST  = .FALSE.
       RETURN
 *
 *     End of SLAMCH
 *
       END
-*
-************************************************************************
-*
-      SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            IEEE1, RND
-      INTEGER            BETA, T
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SLAMC1 determines the machine parameters given by BETA, T, RND, and
-*  IEEE1.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (output) INTEGER
-*          The base of the machine.
-*
-*  T       (output) INTEGER
-*          The number of ( BETA ) digits in the mantissa.
-*
-*  RND     (output) LOGICAL
-*          Specifies whether proper rounding  ( RND = .TRUE. )  or
-*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
-*          be a reliable guide to the way in which the machine performs
-*          its arithmetic.
-*
-*  IEEE1   (output) LOGICAL
-*          Specifies whether rounding appears to be done in the IEEE
-*          'round to nearest' style.
-*
-*  Further Details
-*  ===============
-*
-*  The routine is based on the routine  ENVRON  by Malcolm and
-*  incorporates suggestions by Gentleman and Marovich. See
-*
-*     Malcolm M. A. (1972) Algorithms to reveal properties of
-*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
-*
-*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
-*        that reveal properties of floating point arithmetic units.
-*        Comms. of the ACM, 17, 276-277.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      LOGICAL            FIRST, LIEEE1, LRND
-      INTEGER            LBETA, LT
-      REAL               A, B, C, F, ONE, QTR, SAVEC, T1, T2
-*     ..
-*     .. External Functions ..
-      REAL               SLAMC3
-      EXTERNAL           SLAMC3
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. /
-*     ..
-*     .. Executable Statements ..
-*
-      IF( FIRST ) THEN
-         ONE = 1
-*
-*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
-*        IEEE1, T and RND.
-*
-*        Throughout this routine  we use the function  SLAMC3  to ensure
-*        that relevant values are  stored and not held in registers,  or
-*        are not affected by optimizers.
-*
-*        Compute  a = 2.0**m  with the  smallest positive integer m such
-*        that
-*
-*           fl( a + 1.0 ) = a.
-*
-         A = 1
-         C = 1
-*
-*+       WHILE( C.EQ.ONE )LOOP
-   10    CONTINUE
-         IF( C.EQ.ONE ) THEN
-            A = 2*A
-            C = SLAMC3( A, ONE )
-            C = SLAMC3( C, -A )
-            GO TO 10
-         END IF
-*+       END WHILE
-*
-*        Now compute  b = 2.0**m  with the smallest positive integer m
-*        such that
-*
-*           fl( a + b ) .gt. a.
-*
-         B = 1
-         C = SLAMC3( A, B )
-*
-*+       WHILE( C.EQ.A )LOOP
-   20    CONTINUE
-         IF( C.EQ.A ) THEN
-            B = 2*B
-            C = SLAMC3( A, B )
-            GO TO 20
-         END IF
-*+       END WHILE
-*
-*        Now compute the base.  a and c  are neighbouring floating point
-*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
-*        their difference is beta. Adding 0.25 to c is to ensure that it
-*        is truncated to beta and not ( beta - 1 ).
-*
-         QTR = ONE / 4
-         SAVEC = C
-         C = SLAMC3( C, -A )
-         LBETA = C + QTR
-*
-*        Now determine whether rounding or chopping occurs,  by adding a
-*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
-*
-         B = LBETA
-         F = SLAMC3( B / 2, -B / 100 )
-         C = SLAMC3( F, A )
-         IF( C.EQ.A ) THEN
-            LRND = .TRUE.
-         ELSE
-            LRND = .FALSE.
-         END IF
-         F = SLAMC3( B / 2, B / 100 )
-         C = SLAMC3( F, A )
-         IF( ( LRND ) .AND. ( C.EQ.A ) )
-     $      LRND = .FALSE.
-*
-*        Try and decide whether rounding is done in the  IEEE  'round to
-*        nearest' style. B/2 is half a unit in the last place of the two
-*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
-*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
-*        A, but adding B/2 to SAVEC should change SAVEC.
-*
-         T1 = SLAMC3( B / 2, A )
-         T2 = SLAMC3( B / 2, SAVEC )
-         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
-*
-*        Now find  the  mantissa, t.  It should  be the  integer part of
-*        log to the base beta of a,  however it is safer to determine  t
-*        by powering.  So we find t as the smallest positive integer for
-*        which
-*
-*           fl( beta**t + 1.0 ) = 1.0.
-*
-         LT = 0
-         A = 1
-         C = 1
-*
-*+       WHILE( C.EQ.ONE )LOOP
-   30    CONTINUE
-         IF( C.EQ.ONE ) THEN
-            LT = LT + 1
-            A = A*LBETA
-            C = SLAMC3( A, ONE )
-            C = SLAMC3( C, -A )
-            GO TO 30
-         END IF
-*+       END WHILE
-*
-      END IF
-*
-      BETA = LBETA
-      T = LT
-      RND = LRND
-      IEEE1 = LIEEE1
-      FIRST = .FALSE.
-      RETURN
-*
-*     End of SLAMC1
-*
-      END
-*
-************************************************************************
-*
-      SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            RND
-      INTEGER            BETA, EMAX, EMIN, T
-      REAL               EPS, RMAX, RMIN
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SLAMC2 determines the machine parameters specified in its argument
-*  list.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (output) INTEGER
-*          The base of the machine.
-*
-*  T       (output) INTEGER
-*          The number of ( BETA ) digits in the mantissa.
-*
-*  RND     (output) LOGICAL
-*          Specifies whether proper rounding  ( RND = .TRUE. )  or
-*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
-*          be a reliable guide to the way in which the machine performs
-*          its arithmetic.
-*
-*  EPS     (output) REAL
-*          The smallest positive number such that
-*
-*             fl( 1.0 - EPS ) .LT. 1.0,
-*
-*          where fl denotes the computed value.
-*
-*  EMIN    (output) INTEGER
-*          The minimum exponent before (gradual) underflow occurs.
-*
-*  RMIN    (output) REAL
-*          The smallest normalized number for the machine, given by
-*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
-*          of BETA.
-*
-*  EMAX    (output) INTEGER
-*          The maximum exponent before overflow occurs.
-*
-*  RMAX    (output) REAL
-*          The largest positive number for the machine, given by
-*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
-*          value of BETA.
-*
-*  Further Details
-*  ===============
-*
-*  The computation of  EPS  is based on a routine PARANOIA by
-*  W. Kahan of the University of California at Berkeley.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
-      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
-     $                   NGNMIN, NGPMIN
-      REAL               A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
-     $                   SIXTH, SMALL, THIRD, TWO, ZERO
-*     ..
-*     .. External Functions ..
-      REAL               SLAMC3
-      EXTERNAL           SLAMC3
-*     ..
-*     .. External Subroutines ..
-      EXTERNAL           SLAMC1, SLAMC4, SLAMC5
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          ABS, MAX, MIN
-*     ..
-*     .. Save statement ..
-      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
-     $                   LRMIN, LT
-*     ..
-*     .. Data statements ..
-      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
-*     ..
-*     .. Executable Statements ..
-*
-      IF( FIRST ) THEN
-         ZERO = 0
-         ONE = 1
-         TWO = 2
-*
-*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
-*        BETA, T, RND, EPS, EMIN and RMIN.
-*
-*        Throughout this routine  we use the function  SLAMC3  to ensure
-*        that relevant values are stored  and not held in registers,  or
-*        are not affected by optimizers.
-*
-*        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
-*
-         CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
-*
-*        Start to find EPS.
-*
-         B = LBETA
-         A = B**( -LT )
-         LEPS = A
-*
-*        Try some tricks to see whether or not this is the correct  EPS.
-*
-         B = TWO / 3
-         HALF = ONE / 2
-         SIXTH = SLAMC3( B, -HALF )
-         THIRD = SLAMC3( SIXTH, SIXTH )
-         B = SLAMC3( THIRD, -HALF )
-         B = SLAMC3( B, SIXTH )
-         B = ABS( B )
-         IF( B.LT.LEPS )
-     $      B = LEPS
-*
-         LEPS = 1
-*
-*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
-   10    CONTINUE
-         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
-            LEPS = B
-            C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
-            C = SLAMC3( HALF, -C )
-            B = SLAMC3( HALF, C )
-            C = SLAMC3( HALF, -B )
-            B = SLAMC3( HALF, C )
-            GO TO 10
-         END IF
-*+       END WHILE
-*
-         IF( A.LT.LEPS )
-     $      LEPS = A
-*
-*        Computation of EPS complete.
-*
-*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
-*        Keep dividing  A by BETA until (gradual) underflow occurs. This
-*        is detected when we cannot recover the previous A.
-*
-         RBASE = ONE / LBETA
-         SMALL = ONE
-         DO 20 I = 1, 3
-            SMALL = SLAMC3( SMALL*RBASE, ZERO )
-   20    CONTINUE
-         A = SLAMC3( ONE, SMALL )
-         CALL SLAMC4( NGPMIN, ONE, LBETA )
-         CALL SLAMC4( NGNMIN, -ONE, LBETA )
-         CALL SLAMC4( GPMIN, A, LBETA )
-         CALL SLAMC4( GNMIN, -A, LBETA )
-         IEEE = .FALSE.
-*
-         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
-            IF( NGPMIN.EQ.GPMIN ) THEN
-               LEMIN = NGPMIN
-*            ( Non twos-complement machines, no gradual underflow;
-*              e.g.,  VAX )
-            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
-               LEMIN = NGPMIN - 1 + LT
-               IEEE = .TRUE.
-*            ( Non twos-complement machines, with gradual underflow;
-*              e.g., IEEE standard followers )
-            ELSE
-               LEMIN = MIN( NGPMIN, GPMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
-            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
-               LEMIN = MAX( NGPMIN, NGNMIN )
-*            ( Twos-complement machines, no gradual underflow;
-*              e.g., CYBER 205 )
-            ELSE
-               LEMIN = MIN( NGPMIN, NGNMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
-     $            ( GPMIN.EQ.GNMIN ) ) THEN
-            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
-               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
-*            ( Twos-complement machines with gradual underflow;
-*              no known machine )
-            ELSE
-               LEMIN = MIN( NGPMIN, NGNMIN )
-*            ( A guess; no known machine )
-               IWARN = .TRUE.
-            END IF
-*
-         ELSE
-            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
-*         ( A guess; no known machine )
-            IWARN = .TRUE.
-         END IF
-         FIRST = .FALSE.
-***
-* Comment out this if block if EMIN is ok
-         IF( IWARN ) THEN
-            FIRST = .TRUE.
-            WRITE( 6, FMT = 9999 )LEMIN
-         END IF
-***
-*
-*        Assume IEEE arithmetic if we found denormalised  numbers above,
-*        or if arithmetic seems to round in the  IEEE style,  determined
-*        in routine SLAMC1. A true IEEE machine should have both  things
-*        true; however, faulty machines may have one or the other.
-*
-         IEEE = IEEE .OR. LIEEE1
-*
-*        Compute  RMIN by successive division by  BETA. We could compute
-*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
-*        this computation.
-*
-         LRMIN = 1
-         DO 30 I = 1, 1 - LEMIN
-            LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
-   30    CONTINUE
-*
-*        Finally, call SLAMC5 to compute EMAX and RMAX.
-*
-         CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
-      END IF
-*
-      BETA = LBETA
-      T = LT
-      RND = LRND
-      EPS = LEPS
-      EMIN = LEMIN
-      RMIN = LRMIN
-      EMAX = LEMAX
-      RMAX = LRMAX
-*
-      RETURN
-*
- 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
-     $      '  EMIN = ', I8, /
-     $      ' If, after inspection, the value EMIN looks',
-     $      ' acceptable please comment out ',
-     $      / ' the IF block as marked within the code of routine',
-     $      ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
-*
-*     End of SLAMC2
-*
-      END
-*
 ************************************************************************
 *
       REAL             FUNCTION SLAMC3( A, B )
@@ -608,246 +154,3 @@
       END
 *
 ************************************************************************
-*
-      SUBROUTINE SLAMC4( EMIN, START, BASE )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      INTEGER            BASE
-      INTEGER            EMIN
-      REAL               START
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SLAMC4 is a service routine for SLAMC2.
-*
-*  Arguments
-*  =========
-*
-*  EMIN    (output) INTEGER 
-*          The minimum exponent before (gradual) underflow, computed by
-*          setting A = START and dividing by BASE until the previous A
-*          can not be recovered.
-*
-*  START   (input) REAL
-*          The starting point for determining EMIN.
-*
-*  BASE    (input) INTEGER
-*          The base of the machine.
-*
-* =====================================================================
-*
-*     .. Local Scalars ..
-      INTEGER            I
-      REAL               A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
-*     ..
-*     .. External Functions ..
-      REAL               SLAMC3
-      EXTERNAL           SLAMC3
-*     ..
-*     .. Executable Statements ..
-*
-      A = START
-      ONE = 1
-      RBASE = ONE / BASE
-      ZERO = 0
-      EMIN = 1
-      B1 = SLAMC3( A*RBASE, ZERO )
-      C1 = A
-      C2 = A
-      D1 = A
-      D2 = A
-*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
-*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
-   10 CONTINUE
-      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
-     $    ( D2.EQ.A ) ) THEN
-         EMIN = EMIN - 1
-         A = B1
-         B1 = SLAMC3( A / BASE, ZERO )
-         C1 = SLAMC3( B1*BASE, ZERO )
-         D1 = ZERO
-         DO 20 I = 1, BASE
-            D1 = D1 + B1
-   20    CONTINUE
-         B2 = SLAMC3( A*RBASE, ZERO )
-         C2 = SLAMC3( B2 / RBASE, ZERO )
-         D2 = ZERO
-         DO 30 I = 1, BASE
-            D2 = D2 + B2
-   30    CONTINUE
-         GO TO 10
-      END IF
-*+    END WHILE
-*
-      RETURN
-*
-*     End of SLAMC4
-*
-      END
-*
-************************************************************************
-*
-      SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
-*
-*  -- LAPACK auxiliary routine (version 3.2) --
-*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
-*     November 2006
-*
-*     .. Scalar Arguments ..
-      LOGICAL            IEEE
-      INTEGER            BETA, EMAX, EMIN, P
-      REAL               RMAX
-*     ..
-*
-*  Purpose
-*  =======
-*
-*  SLAMC5 attempts to compute RMAX, the largest machine floating-point
-*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
-*  approximately to a power of 2.  It will fail on machines where this
-*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
-*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
-*  too large (i.e. too close to zero), probably with overflow.
-*
-*  Arguments
-*  =========
-*
-*  BETA    (input) INTEGER
-*          The base of floating-point arithmetic.
-*
-*  P       (input) INTEGER
-*          The number of base BETA digits in the mantissa of a
-*          floating-point value.
-*
-*  EMIN    (input) INTEGER
-*          The minimum exponent before (gradual) underflow.
-*
-*  IEEE    (input) LOGICAL
-*          A logical flag specifying whether or not the arithmetic
-*          system is thought to comply with the IEEE standard.
-*
-*  EMAX    (output) INTEGER
-*          The largest exponent before overflow
-*
-*  RMAX    (output) REAL
-*          The largest machine floating-point number.
-*
-* =====================================================================
-*
-*     .. Parameters ..
-      REAL               ZERO, ONE
-      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
-*     ..
-*     .. Local Scalars ..
-      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
-      REAL               OLDY, RECBAS, Y, Z
-*     ..
-*     .. External Functions ..
-      REAL               SLAMC3
-      EXTERNAL           SLAMC3
-*     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC          MOD
-*     ..
-*     .. Executable Statements ..
-*
-*     First compute LEXP and UEXP, two powers of 2 that bound
-*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
-*     approximately to the bound that is closest to abs(EMIN).
-*     (EMAX is the exponent of the required number RMAX).
-*
-      LEXP = 1
-      EXBITS = 1
-   10 CONTINUE
-      TRY = LEXP*2
-      IF( TRY.LE.( -EMIN ) ) THEN
-         LEXP = TRY
-         EXBITS = EXBITS + 1
-         GO TO 10
-      END IF
-      IF( LEXP.EQ.-EMIN ) THEN
-         UEXP = LEXP
-      ELSE
-         UEXP = TRY
-         EXBITS = EXBITS + 1
-      END IF
-*
-*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
-*     than or equal to EMIN. EXBITS is the number of bits needed to
-*     store the exponent.
-*
-      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
-         EXPSUM = 2*LEXP
-      ELSE
-         EXPSUM = 2*UEXP
-      END IF
-*
-*     EXPSUM is the exponent range, approximately equal to
-*     EMAX - EMIN + 1 .
-*
-      EMAX = EXPSUM + EMIN - 1
-      NBITS = 1 + EXBITS + P
-*
-*     NBITS is the total number of bits needed to store a
-*     floating-point number.
-*
-      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
-*
-*        Either there are an odd number of bits used to store a
-*        floating-point number, which is unlikely, or some bits are
-*        not used in the representation of numbers, which is possible,
-*        (e.g. Cray machines) or the mantissa has an implicit bit,
-*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
-*        most likely. We have to assume the last alternative.
-*        If this is true, then we need to reduce EMAX by one because
-*        there must be some way of representing zero in an implicit-bit
-*        system. On machines like Cray, we are reducing EMAX by one
-*        unnecessarily.
-*
-         EMAX = EMAX - 1
-      END IF
-*
-      IF( IEEE ) THEN
-*
-*        Assume we are on an IEEE machine which reserves one exponent
-*        for infinity and NaN.
-*
-         EMAX = EMAX - 1
-      END IF
-*
-*     Now create RMAX, the largest machine number, which should
-*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
-*
-*     First compute 1.0 - BETA**(-P), being careful that the
-*     result is less than 1.0 .
-*
-      RECBAS = ONE / BETA
-      Z = BETA - ONE
-      Y = ZERO
-      DO 20 I = 1, P
-         Z = Z*RECBAS
-         IF( Y.LT.ONE )
-     $      OLDY = Y
-         Y = SLAMC3( Y, Z )
-   20 CONTINUE
-      IF( Y.GE.ONE )
-     $   Y = OLDY
-*
-*     Now multiply by BETA**EMAX to get RMAX.
-*
-      DO 30 I = 1, EMAX
-         Y = SLAMC3( Y*BETA, ZERO )
-   30 CONTINUE
-*
-      RMAX = Y
-      RETURN
-*
-*     End of SLAMC5
-*
-      END
diff --git a/INSTALL/slamchf77.f b/INSTALL/slamchf77.f
new file mode 100644
index 00000000..5bb2c7be
--- /dev/null
+++ b/INSTALL/slamchf77.f
@@ -0,0 +1,853 @@
+      REAL             FUNCTION SLAMCH( CMACH )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      CHARACTER          CMACH
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMCH determines single precision machine parameters.
+*
+*  Arguments
+*  =========
+*
+*  CMACH   (input) CHARACTER*1
+*          Specifies the value to be returned by SLAMCH:
+*          = 'E' or 'e',   SLAMCH := eps
+*          = 'S' or 's ,   SLAMCH := sfmin
+*          = 'B' or 'b',   SLAMCH := base
+*          = 'P' or 'p',   SLAMCH := eps*base
+*          = 'N' or 'n',   SLAMCH := t
+*          = 'R' or 'r',   SLAMCH := rnd
+*          = 'M' or 'm',   SLAMCH := emin
+*          = 'U' or 'u',   SLAMCH := rmin
+*          = 'L' or 'l',   SLAMCH := emax
+*          = 'O' or 'o',   SLAMCH := rmax
+*
+*          where
+*
+*          eps   = relative machine precision
+*          sfmin = safe minimum, such that 1/sfmin does not overflow
+*          base  = base of the machine
+*          prec  = eps*base
+*          t     = number of (base) digits in the mantissa
+*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
+*          emin  = minimum exponent before (gradual) underflow
+*          rmin  = underflow threshold - base**(emin-1)
+*          emax  = largest exponent before overflow
+*          rmax  = overflow threshold  - (base**emax)*(1-eps)
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LRND
+      INTEGER            BETA, IMAX, IMIN, IT
+      REAL               BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
+     $                   RND, SFMIN, SMALL, T
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLAMC2
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
+     $                   EMAX, RMAX, PREC
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
+         BASE = BETA
+         T = IT
+         IF( LRND ) THEN
+            RND = ONE
+            EPS = ( BASE**( 1-IT ) ) / 2
+         ELSE
+            RND = ZERO
+            EPS = BASE**( 1-IT )
+         END IF
+         PREC = EPS*BASE
+         EMIN = IMIN
+         EMAX = IMAX
+         SFMIN = RMIN
+         SMALL = ONE / RMAX
+         IF( SMALL.GE.SFMIN ) THEN
+*
+*           Use SMALL plus a bit, to avoid the possibility of rounding
+*           causing overflow when computing  1/sfmin.
+*
+            SFMIN = SMALL*( ONE+EPS )
+         END IF
+      END IF
+*
+      IF( LSAME( CMACH, 'E' ) ) THEN
+         RMACH = EPS
+      ELSE IF( LSAME( CMACH, 'S' ) ) THEN
+         RMACH = SFMIN
+      ELSE IF( LSAME( CMACH, 'B' ) ) THEN
+         RMACH = BASE
+      ELSE IF( LSAME( CMACH, 'P' ) ) THEN
+         RMACH = PREC
+      ELSE IF( LSAME( CMACH, 'N' ) ) THEN
+         RMACH = T
+      ELSE IF( LSAME( CMACH, 'R' ) ) THEN
+         RMACH = RND
+      ELSE IF( LSAME( CMACH, 'M' ) ) THEN
+         RMACH = EMIN
+      ELSE IF( LSAME( CMACH, 'U' ) ) THEN
+         RMACH = RMIN
+      ELSE IF( LSAME( CMACH, 'L' ) ) THEN
+         RMACH = EMAX
+      ELSE IF( LSAME( CMACH, 'O' ) ) THEN
+         RMACH = RMAX
+      END IF
+*
+      SLAMCH = RMACH
+      FIRST  = .FALSE.
+      RETURN
+*
+*     End of SLAMCH
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE1, RND
+      INTEGER            BETA, T
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMC1 determines the machine parameters given by BETA, T, RND, and
+*  IEEE1.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  IEEE1   (output) LOGICAL
+*          Specifies whether rounding appears to be done in the IEEE
+*          'round to nearest' style.
+*
+*  Further Details
+*  ===============
+*
+*  The routine is based on the routine  ENVRON  by Malcolm and
+*  incorporates suggestions by Gentleman and Marovich. See
+*
+*     Malcolm M. A. (1972) Algorithms to reveal properties of
+*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
+*
+*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
+*        that reveal properties of floating point arithmetic units.
+*        Comms. of the ACM, 17, 276-277.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, LIEEE1, LRND
+      INTEGER            LBETA, LT
+      REAL               A, B, C, F, ONE, QTR, SAVEC, T1, T2
+*     ..
+*     .. External Functions ..
+      REAL               SLAMC3
+      EXTERNAL           SLAMC3
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         ONE = 1
+*
+*        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
+*        IEEE1, T and RND.
+*
+*        Throughout this routine  we use the function  SLAMC3  to ensure
+*        that relevant values are  stored and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        Compute  a = 2.0**m  with the  smallest positive integer m such
+*        that
+*
+*           fl( a + 1.0 ) = a.
+*
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   10    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            A = 2*A
+            C = SLAMC3( A, ONE )
+            C = SLAMC3( C, -A )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+*        Now compute  b = 2.0**m  with the smallest positive integer m
+*        such that
+*
+*           fl( a + b ) .gt. a.
+*
+         B = 1
+         C = SLAMC3( A, B )
+*
+*+       WHILE( C.EQ.A )LOOP
+   20    CONTINUE
+         IF( C.EQ.A ) THEN
+            B = 2*B
+            C = SLAMC3( A, B )
+            GO TO 20
+         END IF
+*+       END WHILE
+*
+*        Now compute the base.  a and c  are neighbouring floating point
+*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
+*        their difference is beta. Adding 0.25 to c is to ensure that it
+*        is truncated to beta and not ( beta - 1 ).
+*
+         QTR = ONE / 4
+         SAVEC = C
+         C = SLAMC3( C, -A )
+         LBETA = C + QTR
+*
+*        Now determine whether rounding or chopping occurs,  by adding a
+*        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
+*
+         B = LBETA
+         F = SLAMC3( B / 2, -B / 100 )
+         C = SLAMC3( F, A )
+         IF( C.EQ.A ) THEN
+            LRND = .TRUE.
+         ELSE
+            LRND = .FALSE.
+         END IF
+         F = SLAMC3( B / 2, B / 100 )
+         C = SLAMC3( F, A )
+         IF( ( LRND ) .AND. ( C.EQ.A ) )
+     $      LRND = .FALSE.
+*
+*        Try and decide whether rounding is done in the  IEEE  'round to
+*        nearest' style. B/2 is half a unit in the last place of the two
+*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
+*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
+*        A, but adding B/2 to SAVEC should change SAVEC.
+*
+         T1 = SLAMC3( B / 2, A )
+         T2 = SLAMC3( B / 2, SAVEC )
+         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
+*
+*        Now find  the  mantissa, t.  It should  be the  integer part of
+*        log to the base beta of a,  however it is safer to determine  t
+*        by powering.  So we find t as the smallest positive integer for
+*        which
+*
+*           fl( beta**t + 1.0 ) = 1.0.
+*
+         LT = 0
+         A = 1
+         C = 1
+*
+*+       WHILE( C.EQ.ONE )LOOP
+   30    CONTINUE
+         IF( C.EQ.ONE ) THEN
+            LT = LT + 1
+            A = A*LBETA
+            C = SLAMC3( A, ONE )
+            C = SLAMC3( C, -A )
+            GO TO 30
+         END IF
+*+       END WHILE
+*
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      IEEE1 = LIEEE1
+      FIRST = .FALSE.
+      RETURN
+*
+*     End of SLAMC1
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            RND
+      INTEGER            BETA, EMAX, EMIN, T
+      REAL               EPS, RMAX, RMIN
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMC2 determines the machine parameters specified in its argument
+*  list.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (output) INTEGER
+*          The base of the machine.
+*
+*  T       (output) INTEGER
+*          The number of ( BETA ) digits in the mantissa.
+*
+*  RND     (output) LOGICAL
+*          Specifies whether proper rounding  ( RND = .TRUE. )  or
+*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
+*          be a reliable guide to the way in which the machine performs
+*          its arithmetic.
+*
+*  EPS     (output) REAL
+*          The smallest positive number such that
+*
+*             fl( 1.0 - EPS ) .LT. 1.0,
+*
+*          where fl denotes the computed value.
+*
+*  EMIN    (output) INTEGER
+*          The minimum exponent before (gradual) underflow occurs.
+*
+*  RMIN    (output) REAL
+*          The smallest normalized number for the machine, given by
+*          BASE**( EMIN - 1 ), where  BASE  is the floating point value
+*          of BETA.
+*
+*  EMAX    (output) INTEGER
+*          The maximum exponent before overflow occurs.
+*
+*  RMAX    (output) REAL
+*          The largest positive number for the machine, given by
+*          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
+*          value of BETA.
+*
+*  Further Details
+*  ===============
+*
+*  The computation of  EPS  is based on a routine PARANOIA by
+*  W. Kahan of the University of California at Berkeley.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
+      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
+     $                   NGNMIN, NGPMIN
+      REAL               A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
+     $                   SIXTH, SMALL, THIRD, TWO, ZERO
+*     ..
+*     .. External Functions ..
+      REAL               SLAMC3
+      EXTERNAL           SLAMC3
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SLAMC1, SLAMC4, SLAMC5
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, MAX, MIN
+*     ..
+*     .. Save statement ..
+      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
+     $                   LRMIN, LT
+*     ..
+*     .. Data statements ..
+      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
+*     ..
+*     .. Executable Statements ..
+*
+      IF( FIRST ) THEN
+         ZERO = 0
+         ONE = 1
+         TWO = 2
+*
+*        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
+*        BETA, T, RND, EPS, EMIN and RMIN.
+*
+*        Throughout this routine  we use the function  SLAMC3  to ensure
+*        that relevant values are stored  and not held in registers,  or
+*        are not affected by optimizers.
+*
+*        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
+*
+         CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
+*
+*        Start to find EPS.
+*
+         B = LBETA
+         A = B**( -LT )
+         LEPS = A
+*
+*        Try some tricks to see whether or not this is the correct  EPS.
+*
+         B = TWO / 3
+         HALF = ONE / 2
+         SIXTH = SLAMC3( B, -HALF )
+         THIRD = SLAMC3( SIXTH, SIXTH )
+         B = SLAMC3( THIRD, -HALF )
+         B = SLAMC3( B, SIXTH )
+         B = ABS( B )
+         IF( B.LT.LEPS )
+     $      B = LEPS
+*
+         LEPS = 1
+*
+*+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
+   10    CONTINUE
+         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
+            LEPS = B
+            C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
+            C = SLAMC3( HALF, -C )
+            B = SLAMC3( HALF, C )
+            C = SLAMC3( HALF, -B )
+            B = SLAMC3( HALF, C )
+            GO TO 10
+         END IF
+*+       END WHILE
+*
+         IF( A.LT.LEPS )
+     $      LEPS = A
+*
+*        Computation of EPS complete.
+*
+*        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
+*        Keep dividing  A by BETA until (gradual) underflow occurs. This
+*        is detected when we cannot recover the previous A.
+*
+         RBASE = ONE / LBETA
+         SMALL = ONE
+         DO 20 I = 1, 3
+            SMALL = SLAMC3( SMALL*RBASE, ZERO )
+   20    CONTINUE
+         A = SLAMC3( ONE, SMALL )
+         CALL SLAMC4( NGPMIN, ONE, LBETA )
+         CALL SLAMC4( NGNMIN, -ONE, LBETA )
+         CALL SLAMC4( GPMIN, A, LBETA )
+         CALL SLAMC4( GNMIN, -A, LBETA )
+         IEEE = .FALSE.
+*
+         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( NGPMIN.EQ.GPMIN ) THEN
+               LEMIN = NGPMIN
+*            ( Non twos-complement machines, no gradual underflow;
+*              e.g.,  VAX )
+            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
+               LEMIN = NGPMIN - 1 + LT
+               IEEE = .TRUE.
+*            ( Non twos-complement machines, with gradual underflow;
+*              e.g., IEEE standard followers )
+            ELSE
+               LEMIN = MIN( NGPMIN, GPMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
+            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN )
+*            ( Twos-complement machines, no gradual underflow;
+*              e.g., CYBER 205 )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
+     $            ( GPMIN.EQ.GNMIN ) ) THEN
+            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
+               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
+*            ( Twos-complement machines with gradual underflow;
+*              no known machine )
+            ELSE
+               LEMIN = MIN( NGPMIN, NGNMIN )
+*            ( A guess; no known machine )
+               IWARN = .TRUE.
+            END IF
+*
+         ELSE
+            LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
+*         ( A guess; no known machine )
+            IWARN = .TRUE.
+         END IF
+         FIRST = .FALSE.
+***
+* Comment out this if block if EMIN is ok
+         IF( IWARN ) THEN
+            FIRST = .TRUE.
+            WRITE( 6, FMT = 9999 )LEMIN
+         END IF
+***
+*
+*        Assume IEEE arithmetic if we found denormalised  numbers above,
+*        or if arithmetic seems to round in the  IEEE style,  determined
+*        in routine SLAMC1. A true IEEE machine should have both  things
+*        true; however, faulty machines may have one or the other.
+*
+         IEEE = IEEE .OR. LIEEE1
+*
+*        Compute  RMIN by successive division by  BETA. We could compute
+*        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
+*        this computation.
+*
+         LRMIN = 1
+         DO 30 I = 1, 1 - LEMIN
+            LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
+   30    CONTINUE
+*
+*        Finally, call SLAMC5 to compute EMAX and RMAX.
+*
+         CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
+      END IF
+*
+      BETA = LBETA
+      T = LT
+      RND = LRND
+      EPS = LEPS
+      EMIN = LEMIN
+      RMIN = LRMIN
+      EMAX = LEMAX
+      RMAX = LRMAX
+*
+      RETURN
+*
+ 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
+     $      '  EMIN = ', I8, /
+     $      ' If, after inspection, the value EMIN looks',
+     $      ' acceptable please comment out ',
+     $      / ' the IF block as marked within the code of routine',
+     $      ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
+*
+*     End of SLAMC2
+*
+      END
+*
+************************************************************************
+*
+      REAL             FUNCTION SLAMC3( A, B )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      REAL               A, B
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMC3  is intended to force  A  and  B  to be stored prior to doing
+*  the addition of  A  and  B ,  for use in situations where optimizers
+*  might hold one of these in a register.
+*
+*  Arguments
+*  =========
+*
+*  A       (input) REAL
+*  B       (input) REAL
+*          The values A and B.
+*
+* =====================================================================
+*
+*     .. Executable Statements ..
+*
+      SLAMC3 = A + B
+*
+      RETURN
+*
+*     End of SLAMC3
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE SLAMC4( EMIN, START, BASE )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      INTEGER            BASE
+      INTEGER            EMIN
+      REAL               START
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMC4 is a service routine for SLAMC2.
+*
+*  Arguments
+*  =========
+*
+*  EMIN    (output) INTEGER 
+*          The minimum exponent before (gradual) underflow, computed by
+*          setting A = START and dividing by BASE until the previous A
+*          can not be recovered.
+*
+*  START   (input) REAL
+*          The starting point for determining EMIN.
+*
+*  BASE    (input) INTEGER
+*          The base of the machine.
+*
+* =====================================================================
+*
+*     .. Local Scalars ..
+      INTEGER            I
+      REAL               A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
+*     ..
+*     .. External Functions ..
+      REAL               SLAMC3
+      EXTERNAL           SLAMC3
+*     ..
+*     .. Executable Statements ..
+*
+      A = START
+      ONE = 1
+      RBASE = ONE / BASE
+      ZERO = 0
+      EMIN = 1
+      B1 = SLAMC3( A*RBASE, ZERO )
+      C1 = A
+      C2 = A
+      D1 = A
+      D2 = A
+*+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
+*    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
+   10 CONTINUE
+      IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
+     $    ( D2.EQ.A ) ) THEN
+         EMIN = EMIN - 1
+         A = B1
+         B1 = SLAMC3( A / BASE, ZERO )
+         C1 = SLAMC3( B1*BASE, ZERO )
+         D1 = ZERO
+         DO 20 I = 1, BASE
+            D1 = D1 + B1
+   20    CONTINUE
+         B2 = SLAMC3( A*RBASE, ZERO )
+         C2 = SLAMC3( B2 / RBASE, ZERO )
+         D2 = ZERO
+         DO 30 I = 1, BASE
+            D2 = D2 + B2
+   30    CONTINUE
+         GO TO 10
+      END IF
+*+    END WHILE
+*
+      RETURN
+*
+*     End of SLAMC4
+*
+      END
+*
+************************************************************************
+*
+      SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
+*
+*  -- LAPACK auxiliary routine (version 3.2) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     .. Scalar Arguments ..
+      LOGICAL            IEEE
+      INTEGER            BETA, EMAX, EMIN, P
+      REAL               RMAX
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SLAMC5 attempts to compute RMAX, the largest machine floating-point
+*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
+*  approximately to a power of 2.  It will fail on machines where this
+*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
+*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
+*  too large (i.e. too close to zero), probably with overflow.
+*
+*  Arguments
+*  =========
+*
+*  BETA    (input) INTEGER
+*          The base of floating-point arithmetic.
+*
+*  P       (input) INTEGER
+*          The number of base BETA digits in the mantissa of a
+*          floating-point value.
+*
+*  EMIN    (input) INTEGER
+*          The minimum exponent before (gradual) underflow.
+*
+*  IEEE    (input) LOGICAL
+*          A logical flag specifying whether or not the arithmetic
+*          system is thought to comply with the IEEE standard.
+*
+*  EMAX    (output) INTEGER
+*          The largest exponent before overflow
+*
+*  RMAX    (output) REAL
+*          The largest machine floating-point number.
+*
+* =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
+      REAL               OLDY, RECBAS, Y, Z
+*     ..
+*     .. External Functions ..
+      REAL               SLAMC3
+      EXTERNAL           SLAMC3
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MOD
+*     ..
+*     .. Executable Statements ..
+*
+*     First compute LEXP and UEXP, two powers of 2 that bound
+*     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
+*     approximately to the bound that is closest to abs(EMIN).
+*     (EMAX is the exponent of the required number RMAX).
+*
+      LEXP = 1
+      EXBITS = 1
+   10 CONTINUE
+      TRY = LEXP*2
+      IF( TRY.LE.( -EMIN ) ) THEN
+         LEXP = TRY
+         EXBITS = EXBITS + 1
+         GO TO 10
+      END IF
+      IF( LEXP.EQ.-EMIN ) THEN
+         UEXP = LEXP
+      ELSE
+         UEXP = TRY
+         EXBITS = EXBITS + 1
+      END IF
+*
+*     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
+*     than or equal to EMIN. EXBITS is the number of bits needed to
+*     store the exponent.
+*
+      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
+         EXPSUM = 2*LEXP
+      ELSE
+         EXPSUM = 2*UEXP
+      END IF
+*
+*     EXPSUM is the exponent range, approximately equal to
+*     EMAX - EMIN + 1 .
+*
+      EMAX = EXPSUM + EMIN - 1
+      NBITS = 1 + EXBITS + P
+*
+*     NBITS is the total number of bits needed to store a
+*     floating-point number.
+*
+      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
+*
+*        Either there are an odd number of bits used to store a
+*        floating-point number, which is unlikely, or some bits are
+*        not used in the representation of numbers, which is possible,
+*        (e.g. Cray machines) or the mantissa has an implicit bit,
+*        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
+*        most likely. We have to assume the last alternative.
+*        If this is true, then we need to reduce EMAX by one because
+*        there must be some way of representing zero in an implicit-bit
+*        system. On machines like Cray, we are reducing EMAX by one
+*        unnecessarily.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+      IF( IEEE ) THEN
+*
+*        Assume we are on an IEEE machine which reserves one exponent
+*        for infinity and NaN.
+*
+         EMAX = EMAX - 1
+      END IF
+*
+*     Now create RMAX, the largest machine number, which should
+*     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
+*
+*     First compute 1.0 - BETA**(-P), being careful that the
+*     result is less than 1.0 .
+*
+      RECBAS = ONE / BETA
+      Z = BETA - ONE
+      Y = ZERO
+      DO 20 I = 1, P
+         Z = Z*RECBAS
+         IF( Y.LT.ONE )
+     $      OLDY = Y
+         Y = SLAMC3( Y, Z )
+   20 CONTINUE
+      IF( Y.GE.ONE )
+     $   Y = OLDY
+*
+*     Now multiply by BETA**EMAX to get RMAX.
+*
+      DO 30 I = 1, EMAX
+         Y = SLAMC3( Y*BETA, ZERO )
+   30 CONTINUE
+*
+      RMAX = Y
+      RETURN
+*
+*     End of SLAMC5
+*
+      END