LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ slamc5()

subroutine slamc5 ( integer  BETA,
integer  P,
integer  EMIN,
logical  IEEE,
integer  EMAX,
real  RMAX 
)

SLAMC5

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.
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 802 of file slamchf77.f.

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