LAPACK  3.10.1
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 792 of file slamchf77.f.

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