LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlamc5 ( integer  BETA,
integer  P,
integer  EMIN,
logical  IEEE,
integer  EMAX,
double precision  RMAX 
)

DLAMC5

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.
Parameters
[in]BETA
          The base of floating-point arithmetic.
[in]P
          The number of base BETA digits in the mantissa of a
          floating-point value.
[in]EMIN
          The minimum exponent before (gradual) underflow.
[in]IEEE
          A logical flag specifying whether or not the arithmetic
          system is thought to comply with the IEEE standard.
[out]EMAX
          The largest exponent before overflow
[out]RMAX
          The largest machine floating-point number.

Definition at line 797 of file dlamchf77.f.

797 *
798 * -- LAPACK auxiliary routine (version 3.4.1) --
799 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
800 * November 2010
801 *
802 * .. Scalar Arguments ..
803  LOGICAL ieee
804  INTEGER beta, emax, emin, p
805  DOUBLE PRECISION rmax
806 * ..
807 * =====================================================================
808 *
809 * .. Parameters ..
810  DOUBLE PRECISION zero, one
811  parameter ( zero = 0.0d0, one = 1.0d0 )
812 * ..
813 * .. Local Scalars ..
814  INTEGER exbits, expsum, i, lexp, nbits, try, uexp
815  DOUBLE PRECISION oldy, recbas, y, z
816 * ..
817 * .. External Functions ..
818  DOUBLE PRECISION dlamc3
819  EXTERNAL dlamc3
820 * ..
821 * .. Intrinsic Functions ..
822  INTRINSIC mod
823 * ..
824 * .. Executable Statements ..
825 *
826 * First compute LEXP and UEXP, two powers of 2 that bound
827 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
828 * approximately to the bound that is closest to abs(EMIN).
829 * (EMAX is the exponent of the required number RMAX).
830 *
831  lexp = 1
832  exbits = 1
833  10 CONTINUE
834  try = lexp*2
835  IF( try.LE.( -emin ) ) THEN
836  lexp = try
837  exbits = exbits + 1
838  GO TO 10
839  END IF
840  IF( lexp.EQ.-emin ) THEN
841  uexp = lexp
842  ELSE
843  uexp = try
844  exbits = exbits + 1
845  END IF
846 *
847 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
848 * than or equal to EMIN. EXBITS is the number of bits needed to
849 * store the exponent.
850 *
851  IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
852  expsum = 2*lexp
853  ELSE
854  expsum = 2*uexp
855  END IF
856 *
857 * EXPSUM is the exponent range, approximately equal to
858 * EMAX - EMIN + 1 .
859 *
860  emax = expsum + emin - 1
861  nbits = 1 + exbits + p
862 *
863 * NBITS is the total number of bits needed to store a
864 * floating-point number.
865 *
866  IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
867 *
868 * Either there are an odd number of bits used to store a
869 * floating-point number, which is unlikely, or some bits are
870 * not used in the representation of numbers, which is possible,
871 * (e.g. Cray machines) or the mantissa has an implicit bit,
872 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
873 * most likely. We have to assume the last alternative.
874 * If this is true, then we need to reduce EMAX by one because
875 * there must be some way of representing zero in an implicit-bit
876 * system. On machines like Cray, we are reducing EMAX by one
877 * unnecessarily.
878 *
879  emax = emax - 1
880  END IF
881 *
882  IF( ieee ) THEN
883 *
884 * Assume we are on an IEEE machine which reserves one exponent
885 * for infinity and NaN.
886 *
887  emax = emax - 1
888  END IF
889 *
890 * Now create RMAX, the largest machine number, which should
891 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
892 *
893 * First compute 1.0 - BETA**(-P), being careful that the
894 * result is less than 1.0 .
895 *
896  recbas = one / beta
897  z = beta - one
898  y = zero
899  DO 20 i = 1, p
900  z = z*recbas
901  IF( y.LT.one )
902  $ oldy = y
903  y = dlamc3( y, z )
904  20 CONTINUE
905  IF( y.GE.one )
906  $ y = oldy
907 *
908 * Now multiply by BETA**EMAX to get RMAX.
909 *
910  DO 30 i = 1, emax
911  y = dlamc3( y*beta, zero )
912  30 CONTINUE
913 *
914  rmax = y
915  RETURN
916 *
917 * End of DLAMC5
918 *
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169

Here is the caller graph for this function: