diff -up blacs/BLACS/INSTALL/Makefile.fedora blacs/BLACS/INSTALL/Makefile
--- blacs/BLACS/INSTALL/Makefile.fedora 2006-01-18 14:36:03.000000000 -0700
+++ blacs/BLACS/INSTALL/Makefile 2014-02-05 15:28:39.268960236 -0700
@@ -28,7 +28,7 @@ $(INSTdir)/xintface : Fintface.o Cintfac
$(INSTdir)/xsyserrors : syserrors.o
$(CCLOADER) $(CCLOADFLAGS) -o $@ syserrors.o $(MPILIB)
-$(INSTdir)/xtc_CsameF77 : mpif.h tc_fCsameF77.o tc_cCsameF77.o
+$(INSTdir)/xtc_CsameF77 : tc_fCsameF77.o tc_cCsameF77.o
$(F77LOADER) $(F77LOADFLAGS) -o $@ tc_fCsameF77.o tc_cCsameF77.o $(MPILIB)
$(INSTdir)/xtc_UseMpich : tc_UseMpich.o
@@ -37,17 +37,13 @@ $(INSTdir)/xtc_UseMpich : tc_UseMpich.o
$(INSTdir)/xcmpi_sane : cmpi_sane.o
$(CCLOADER) $(CCLOADFLAGS) -o $@ cmpi_sane.o $(MPILIB)
-$(INSTdir)/xfmpi_sane : mpif.h fmpi_sane.o
+$(INSTdir)/xfmpi_sane : fmpi_sane.o
$(F77LOADER) $(F77LOADFLAGS) -o $@ fmpi_sane.o $(MPILIB)
clean:
rm -f size.o Fintface.o Cintface.o syserrors.o transcomm.o \
mpi_sane.o fmpi_sane.o tc_UseMpich.o tc_fCsameF77.o tc_cCsameF77.o
-mpif.h : $(MPIINCdir)/mpif.h
- rm -f mpif.h
- ln -s $(MPIINCdir)/mpif.h mpif.h
-
.f.o: ; $(F77) -c $(F77FLAGS) $*.f
.c.o:
$(CC) -c $(CCFLAGS) $(BLACSDEFS) $<
diff -up blacs/BLACS/SRC/MPI/Makefile.fedora blacs/BLACS/SRC/MPI/Makefile
--- blacs/BLACS/SRC/MPI/Makefile.fedora 2007-10-04 11:09:14.000000000 -0600
+++ blacs/BLACS/SRC/MPI/Makefile 2014-02-05 15:28:39.267960243 -0700
@@ -88,12 +88,14 @@ all : INTERN $(Fintobj) $(Cintobj)
make $(BLACSCINIT)
$(ARCH) $(ARCHFLAGS) $(BLACSLIB) $(Fintobj) $(Cintobj)
$(RANLIB) $(BLACSLIB)
+ ( mkdir tmp ; cp $(BLACSLIB) tmp ; cd tmp ; ar x $(BLACSLIB) ; mpif77 -shared -o $(BLACSdir)/$(BLACSLIBS).0.0 *.o *.oo -Wl,-soname=$(BLACSLIBS) ; cd .. )
$(BLACSFINIT) :
$(CC) -o Cblacs_pinfo.o -c $(CCFLAGS) $(BLACSDEFS) -DCallFromC -DMainInF77 \
blacs_pinfo_.c
$(CC) -c $(CCFLAGS) $(BLACSDEFS) -DMainInF77 blacs_pinfo_.c
$(ARCH) $(ARCHFLAGS) $(BLACSFINIT) blacs_pinfo_.o Cblacs_pinfo.o
+ $(CC) $(CCFLAGS) -shared -Wl,-soname=$(BLACSFINITS) -o $(BLACSdir)/$(BLACSFINITS).0.0 blacs_pinfo_.o Cblacs_pinfo.o
$(RANLIB) $(BLACSFINIT)
$(BLACSCINIT) :
@@ -101,6 +103,7 @@ $(BLACSCINIT) :
blacs_pinfo_.c
$(CC) -c $(CCFLAGS) $(BLACSDEFS) -DMainInC blacs_pinfo_.c
$(ARCH) $(ARCHFLAGS) $(BLACSCINIT) blacs_pinfo_.o Cblacs_pinfo.o
+ $(CC) $(CCFLAGS) -shared -Wl,-soname=$(BLACSCINITS) -o $(BLACSdir)/$(BLACSCINITS).0.0 blacs_pinfo_.o Cblacs_pinfo.o
$(RANLIB) $(BLACSCINIT)
# ------------------
@@ -129,7 +132,7 @@ clean :
( cd INTERNAL ; rm -f bi_f77_get_constants.o bi_f77_mpi_attr_get.o \
bi_f77_mpi_initialized.o bi_f77_mpi_op_create.o \
bi_f77_mpi_test.o bi_f77_mpi_testall.o \
- bi_f77_init.o Bdef.h Bconfig.h mpif.h )
+ bi_f77_init.o Bdef.h Bconfig.h )
( cd INTERNAL ; rm -f $(internal) )
# -------------------------------------
@@ -142,7 +145,7 @@ killib :
( cd INTERNAL ; rm -f $(internal:.o=.c) Bdef.h Bconfig.h )
( cd INTERNAL ; rm -f bi_f77_get_constants.f bi_f77_mpi_attr_get.f \
bi_f77_mpi_initialized.f bi_f77_mpi_op_create.f \
- bi_f77_mpi_test.f bi_f77_mpi_testall.f mpif.h )
+ bi_f77_mpi_test.f bi_f77_mpi_testall.f )
# -------------------------------------------------------------------------
# Establish how to make logical links to the long-name C interface routines
@@ -180,25 +183,21 @@ Cfree_blacs_system_handle_.oo : free_bla
# -------------------------------------
# Compile the (ouch!) fortran internals
# -------------------------------------
-bi_f77_init.o : mpif.h bi_f77_init.f
+bi_f77_init.o : bi_f77_init.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_get_constants.o : mpif.h bi_f77_get_constants.f
+bi_f77_get_constants.o :bi_f77_get_constants.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_mpi_attr_get.o : mpif.h bi_f77_mpi_attr_get.f
+bi_f77_mpi_attr_get.o : bi_f77_mpi_attr_get.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_mpi_initialized.o : mpif.h bi_f77_mpi_initialized.f
+bi_f77_mpi_initialized.o : bi_f77_mpi_initialized.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_mpi_op_create.o : mpif.h bi_f77_mpi_op_create.f
+bi_f77_mpi_op_create.o : bi_f77_mpi_op_create.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_mpi_test.o : mpif.h bi_f77_mpi_test.f
+bi_f77_mpi_test.o : bi_f77_mpi_test.f
$(F77) -c $(F77FLAGS) $*.f
-bi_f77_mpi_testall.o : mpif.h bi_f77_mpi_testall.f
+bi_f77_mpi_testall.o : bi_f77_mpi_testall.f
$(F77) -c $(F77FLAGS) $*.f
-mpif.h : $(MPIINCdir)/mpif.h
- rm -f mpif.h
- ln -s $(MPIINCdir)/mpif.h mpif.h
-
# ------------------------------------------------------------------------
# We move C .o files to .oo so that we can use the portable suffix rule for
# compilation, and still have them coexist with the fortran interface
diff -up blacs/BLACS/TESTING/Makefile.fedora blacs/BLACS/TESTING/Makefile
--- blacs/BLACS/TESTING/Makefile.fedora 2006-01-18 14:36:03.000000000 -0700
+++ blacs/BLACS/TESTING/Makefile 2014-02-05 15:54:55.718550117 -0700
@@ -18,10 +18,10 @@ obj = blacstest.o btprim_$(COMMLIB).o
$(CTESTexe): $(obj) $(tools)
$(CC) -c $(CCFLAGS) -DBTCINTFACE $(BLACSDEFS) Cbt.c
- $(F77LOADER) $(F77LOADFLAGS) -o $@ $(obj) $(tools) Cbt.o $(BTLIBS)
+ $(F77LOADER) $(F77LOADFLAGS) -o $@ $(obj) $(tools) Cbt.o $(BTLIBS) -llapack
$(FTESTexe): $(obj) $(tools)
- $(F77LOADER) $(F77LOADFLAGS) -o $@ $(obj) $(tools) $(BTLIBS)
+ $(F77LOADER) $(F77LOADFLAGS) -o $@ $(obj) $(tools) $(BTLIBS) -llapack
# --------------------------------------------------------------------
# The files tools.f and blacstest.f are compiled without optimization.
@@ -51,22 +51,17 @@ $(TESTdir)/comb.dat : $(BTOPdir)/TESTING
cp $(BTOPdir)/TESTING/comb.dat $(TESTdir)/
btprim_MPI.o : btprim_MPI.f
- make mpif.h
$(F77) -c $(F77FLAGS) $*.f
btprim_PVM.o : btprim_PVM.f
make fpvm3.h
$(F77) -c $(F77FLAGS) $*.f
-mpif.h : $(MPIINCdir)/mpif.h
- rm -f mpif.h
- ln -s $(MPIINCdir)/mpif.h mpif.h
-
fpvm3.h : $(PVMINCdir)/fpvm3.h
rm -f fpvm3.h
ln -s $(PVMINCdir)/fpvm3.h fpvm3.h
clean :
- rm -f $(obj) tools.o Cbt.o mpif.h fpvm3.h
+ rm -f $(obj) tools.o Cbt.o fpvm3.h
.f.o: ; $(F77) -c $(F77FLAGS) $*.f
diff -up blacs/BLACS/TESTING/tools.f.fedora blacs/BLACS/TESTING/tools.f
--- blacs/BLACS/TESTING/tools.f.fedora 2006-01-18 14:36:03.000000000 -0700
+++ blacs/BLACS/TESTING/tools.f 2014-02-05 15:54:36.031749119 -0700
@@ -6,1807 +6,6 @@
* 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) --