tools.f

Go to the documentation of this file.

*  ================================================================
*  This file contains the following LAPACK routines, for use by the
*  BLACS tester:  LSAME, SLAMCH, DLAMCH, DLARND, ZLARND, DLARAN,
*  and ZLARAN. If you have ScaLAPACK or LAPACK, all of these files
*  are present in your library, and you may discard this file and 
*  point to the appropriate archive instead.
*  ================================================================
 
      DOUBLE PRECISION FUNCTION dlamch( CMACH )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
      RETURN
*
*     End of DLAMCH
*
      END
*
************************************************************************
*
      SUBROUTINE dlamc1( BETA, T, RND, IEEE1 )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
      RETURN
*
*     End of DLAMC1
*
      END
*
************************************************************************
*
      SUBROUTINE dlamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
***
* 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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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, 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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      INTEGER            BASE, EMIN
      DOUBLE PRECISION   START
*     ..
*
*  Purpose
*  =======
*
*  DLAMC4 is a service routine for DLAMC2.
*
*  Arguments
*  =========
*
*  EMIN    (output) EMIN
*          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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
      REAL             function slamch( cmach )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
      RETURN
*
*     End of SLAMCH
*
      END
*
************************************************************************
*
      SUBROUTINE slamc1( BETA, T, RND, IEEE1 )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
      RETURN
*
*     End of SLAMC1
*
      END
*
************************************************************************
*
      SUBROUTINE slamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
         first = .false.
         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
***
* 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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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, 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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. Scalar Arguments ..
      INTEGER            BASE, EMIN
      REAL               START
*     ..
*
*  Purpose
*  =======
*
*  SLAMC4 is a service routine for SLAMC2.
*
*  Arguments
*  =========
*
*  EMIN    (output) EMIN
*          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 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*     .. 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
      LOGICAL          FUNCTION lsame( CA, CB )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1994
*
*     .. Scalar Arguments ..
      CHARACTER          ca, cb
*     ..
*
*  Purpose
*  =======
*
*  LSAME returns .TRUE. if CA is the same letter as CB regardless of
*  case.
*
*  Arguments
*  =========
*
*  CA      (input) CHARACTER*1
*  CB      (input) CHARACTER*1
*          CA and CB specify the single characters to be compared.
*
* =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          ichar
*     ..
*     .. Local Scalars ..
      INTEGER            inta, intb, zcode
*     ..
*     .. Executable Statements ..
*
*     Test if the characters are equal
*
      lsame = ca.EQ.cb
      IF( lsame )
     $   RETURN
*
*     Now test for equivalence if both characters are alphabetic.
*
      zcode = ichar( 'Z' )
*
*     Use 'Z' rather than 'A' so that ASCII can be detected on Prime
*     machines, on which ICHAR returns a value with bit 8 set.
*     ICHAR('A') on Prime machines returns 193 which is the same as
*     ICHAR('A') on an EBCDIC machine.
*
      inta = ichar( ca )
      intb = ichar( cb )
*
      IF( zcode.EQ.90 .OR. zcode.EQ.122 ) THEN
*
*        ASCII is assumed - ZCODE is the ASCII code of either lower or
*        upper case 'Z'.
*
         IF( inta.GE.97 .AND. inta.LE.122 ) inta = inta - 32
         IF( intb.GE.97 .AND. intb.LE.122 ) intb = intb - 32
*
      ELSE IF( zcode.EQ.233 .OR. zcode.EQ.169 ) THEN
*
*        EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
*        upper case 'Z'.
*
         IF( inta.GE.129 .AND. inta.LE.137 .OR.
     $       inta.GE.145 .AND. inta.LE.153 .OR.
     $       inta.GE.162 .AND. inta.LE.169 ) inta = inta + 64
         IF( intb.GE.129 .AND. intb.LE.137 .OR.
     $       intb.GE.145 .AND. intb.LE.153 .OR.
     $       intb.GE.162 .AND. intb.LE.169 ) intb = intb + 64
*
      ELSE IF( zcode.EQ.218 .OR. zcode.EQ.250 ) THEN
*
*        ASCII is assumed, on Prime machines - ZCODE is the ASCII code
*        plus 128 of either lower or upper case 'Z'.
*
         IF( inta.GE.225 .AND. inta.LE.250 ) inta = inta - 32
         IF( intb.GE.225 .AND. intb.LE.250 ) intb = intb - 32
      END IF
      lsame = inta.EQ.intb
*
*     RETURN
*
*     End of LSAME
*
      END
      DOUBLE PRECISION FUNCTION dlarnd( IDIST, ISEED )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1994
*
*     .. Scalar Arguments ..
      INTEGER            idist
*     ..
*     .. Array Arguments ..
      INTEGER            iseed( 4 )
*     ..
*
*  Purpose
*  =======
*
*  DLARND returns a random real number from a uniform or normal
*  distribution.
*
*  Arguments
*  =========
*
*  IDIST   (input) INTEGER
*          Specifies the distribution of the random numbers:
*          = 1:  uniform (0,1)
*          = 2:  uniform (-1,1)
*          = 3:  normal (0,1)
*
*  ISEED   (input/output) INTEGER array, dimension (4)
*          On entry, the seed of the random number generator; the array
*          elements must be between 0 and 4095, and ISEED(4) must be
*          odd.
*          On exit, the seed is updated.
*
*  Further Details
*  ===============
*
*  This routine calls the auxiliary routine DLARAN to generate a random
*  real number from a uniform (0,1) distribution. The Box-Muller method
*  is used to transform numbers from a uniform to a normal distribution.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   one, two
      parameter( one = 1.0d+0, two = 2.0d+0 )
      DOUBLE PRECISION   twopi
      parameter( twopi = 6.2831853071795864769252867663d+0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   t1, t2
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlaran
      EXTERNAL           dlaran
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cos, log, sqrt
*     ..
*     .. Executable Statements ..
*
*     Generate a real random number from a uniform (0,1) distribution
*
      t1 = dlaran( iseed )
*
      IF( idist.EQ.1 ) THEN
*
*        uniform (0,1)
*
         dlarnd = t1
      ELSE IF( idist.EQ.2 ) THEN
*
*        uniform (-1,1)
*
         dlarnd = two*t1 - one
      ELSE IF( idist.EQ.3 ) THEN
*
*        normal (0,1)
*
         t2 = dlaran( iseed )
         dlarnd = sqrt( -two*log( t1 ) )*cos( twopi*t2 )
      END IF
      RETURN
*
*     End of DLARND
*
      END
      DOUBLE COMPLEX   FUNCTION zlarnd( IDIST, ISEED )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1994
*
*     .. Scalar Arguments ..
      INTEGER            idist
*     ..
*     .. Array Arguments ..
      INTEGER            iseed( 4 )
*     ..
*
*  Purpose
*  =======
*
*  ZLARND returns a random complex number from a uniform or normal
*  distribution.
*
*  Arguments
*  =========
*
*  IDIST   (input) INTEGER
*          Specifies the distribution of the random numbers:
*          = 1:  real and imaginary parts each uniform (0,1)
*          = 2:  real and imaginary parts each uniform (-1,1)
*          = 3:  real and imaginary parts each normal (0,1)
*          = 4:  uniformly distributed on the disc abs(z) <= 1
*          = 5:  uniformly distributed on the circle abs(z) = 1
*
*  ISEED   (input/output) INTEGER array, dimension (4)
*          On entry, the seed of the random number generator; the array
*          elements must be between 0 and 4095, and ISEED(4) must be
*          odd.
*          On exit, the seed is updated.
*
*  Further Details
*  ===============
*
*  This routine calls the auxiliary routine DLARAN to generate a random
*  real number from a uniform (0,1) distribution. The Box-Muller method
*  is used to transform numbers from a uniform to a normal distribution.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, one, two
      parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
      DOUBLE PRECISION   twopi
      parameter( twopi = 6.2831853071795864769252867663d+0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   t1, t2
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlaran
      EXTERNAL           dlaran
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dcmplx, exp, log, sqrt
*     ..
*     .. Executable Statements ..
*
*     Generate a pair of real random numbers from a uniform (0,1)
*     distribution
*
      t1 = dlaran( iseed )
      t2 = dlaran( iseed )
*
      IF( idist.EQ.1 ) THEN
*
*        real and imaginary parts each uniform (0,1)
*
         zlarnd = dcmplx( t1, t2 )
      ELSE IF( idist.EQ.2 ) THEN
*
*        real and imaginary parts each uniform (-1,1)
*
         zlarnd = dcmplx( two*t1-one, two*t2-one )
      ELSE IF( idist.EQ.3 ) THEN
*
*        real and imaginary parts each normal (0,1)
*
         zlarnd = sqrt( -two*log( t1 ) )*exp( dcmplx( zero, twopi*t2 ) )
      ELSE IF( idist.EQ.4 ) THEN
*
*        uniform distribution on the unit disc abs(z) <= 1
*
         zlarnd = sqrt( t1 )*exp( dcmplx( zero, twopi*t2 ) )
      ELSE IF( idist.EQ.5 ) THEN
*
*        uniform distribution on the unit circle abs(z) = 1
*
         zlarnd = exp( dcmplx( zero, twopi*t2 ) )
      END IF
      RETURN
*
*     End of ZLARND
*
      END
      DOUBLE PRECISION FUNCTION dlaran( ISEED )
*
*  -- LAPACK auxiliary routine (version 2.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     February 29, 1992
*
*     .. Array Arguments ..
      INTEGER            iseed( 4 )
*     ..
*
*  Purpose
*  =======
*
*  DLARAN returns a random real number from a uniform (0,1)
*  distribution.
*
*  Arguments
*  =========
*
*  ISEED   (input/output) INTEGER array, dimension (4)
*          On entry, the seed of the random number generator; the array
*          elements must be between 0 and 4095, and ISEED(4) must be
*          odd.
*          On exit, the seed is updated.
*
*  Further Details
*  ===============
*
*  This routine uses a multiplicative congruential method with modulus
*  2**48 and multiplier 33952834046453 (see G.S.Fishman,
*  'Multiplicative congruential random number generators with modulus
*  2**b: an exhaustive analysis for b = 32 and a partial analysis for
*  b = 48', Math. Comp. 189, pp 331-344, 1990).
*
*  48-bit integers are stored in 4 integer array elements with 12 bits
*  per element. Hence the routine is portable across machines with
*  integers of 32 bits or more.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            m1, m2, m3, m4
      parameter( m1 = 494, m2 = 322, m3 = 2508, m4 = 2549 )
      DOUBLE PRECISION   one
      parameter( one = 1.0d+0 )
      INTEGER            ipw2
      DOUBLE PRECISION   r
      parameter( ipw2 = 4096, r = one / ipw2 )
*     ..
*     .. Local Scalars ..
      INTEGER            it1, it2, it3, it4
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, mod
*     ..
*     .. Executable Statements ..
*
*     multiply the seed by the multiplier modulo 2**48
*
      it4 = iseed( 4 )*m4
      it3 = it4 / ipw2
      it4 = it4 - ipw2*it3
      it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
      it2 = it3 / ipw2
      it3 = it3 - ipw2*it2
      it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
      it1 = it2 / ipw2
      it2 = it2 - ipw2*it1
      it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 +
     $      iseed( 4 )*m1
      it1 = mod( it1, ipw2 )
*
*     return updated seed
*
      iseed( 1 ) = it1
      iseed( 2 ) = it2
      iseed( 3 ) = it3
      iseed( 4 ) = it4
*
*     convert 48-bit integer to a real number in the interval (0,1)
*
      dlaran = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
     $         ( dble( it4 ) ) ) ) )
      RETURN
*
*     End of DLARAN
*
      END