The cost of underflows: some empirical data

[David Hough]

Here's a little short program I've used for testing IEEE underflow/subnormal
performance.   Be
careful, not all Fortran compilers handle 2.0e0**(-149) and 2.0d0**(-1074)
correctly.   _PLAIN_UNIX is for systems other than SunOS that don't have
a way of enabling abrupt underflow. 
The timing and random number routines may require some attention, as is often
the case when porting Fortran 77 programs.

#ifdef __hpux
#define STDERR 7
#else
#define STDERR 0
#endif
#define STDOUT 6

#ifdef SP
#define REAL real
#define MINNORM (2.0e0**(-126))
#define MINSUB (2.0e0**(-149))
#endif
#ifdef DP
#define REAL doubleprecision
#define MINNORM (2.0d0**(-1022))
#define MINSUB (2.0d0**(-1074))
#endif
#ifdef QP
#define REAL real*16
#define MINNORM (2.0q0**(-16382))
#define MINSUB (2.0q0**(-16494))
#endif
#define	N 20000	/* size of vectors */
#define M 10  /* repetition count */

	program fastmode
	REAL x(N,3), y(N,3), z(N)
	logical gradual
	integer i
	do 1 i = 1,3
	call setup(i,x(1,i),y(1,i),N)
 1	continue
#ifdef _PLAIN_UNIX
	if (gradual()) then
	write (STDOUT,*) 
     1	' results for gradual underflow/subnormal/exp cases '
	write (STDERR,*) 
     1	' results for gradual underflow/subnormal/exp cases '
	call dotests(M,0,x,y,z,N)
	else
	write (STDOUT,*) 
     1	' results for abrupt underflow/subnormal/exp cases '
	write (STDERR,*) 
     1	' results for abrupt underflow/subnormal/exp cases '
	call dotests(M,1,x,y,z,N)
	endif
#else
	write (STDOUT,*) 
     1	' results for abrupt underflow/subnormal/exp cases '
	write (STDERR,*) 
     1	' results for abrupt underflow/subnormal/exp cases '
	call dotests(M,1,x,y,z,N)
	write (STDOUT,*) 
     1	' results for gradual underflow/subnormal/exp cases '
	write (STDERR,*) 
     1	' results for gradual underflow/subnormal/exp cases '
	call dotests(M,0,x,y,z,N)
#endif
	end

	subroutine dotests(reps,abrupt,x,y,z,n)
	REAL x(n,3), y(n,3), z(n)
	REAL abrupttimes(3)
	integer n, reps, abrupt
	doubleprecision check
	doubleprecision t1, r, u, s
	integer i,j
	real*4 dtime,tarray(2)
	
	do 1 i = 1, 3
#ifndef _PLAIN_UNIX 
	if (abrupt .eq. 1) then
	call abrupt_underflow 
	endif
#endif 
#ifdef __NEWVALID
	call startvalidtime
#else
	t1=dtime(tarray)
#endif
	do 3 j = 1, reps
	call test(i,x(1,i),y(1,i),z,n)
 3	continue
#ifdef __NEWVALID
	call validtime(r,u,s)
	t1=u+s
#else
	t1=dtime(tarray)
#endif
	if (abrupt .eq. 1) then
	abrupttimes(i)=t1
#ifndef _PLAIN_UNIX 
	call gradual_underflow 
#endif
	endif
	check = 0
	do 2 j = 1,N
	check = check + z(j)
c	print *,j,z(j)
 2	continue
	write(STDOUT,*) ' test ',i,' checksum ',check
	if (i .eq. 3) then
		write(STDERR,*) ' test ',i,' KFLOPS ',
     1			nint(10*reps*(N/(1000*(t1))))
c			exp deserves about 10 FLOPS
	else
		write(STDERR,*) ' test ',i,' KFLOPS ',
     1			nint(reps*(N/(1000*(t1))))
	endif
#ifndef _PLAIN_UNIX
	if (abrupt.eq.0) then
		write(STDERR,*) 
     1	' gradual/abrupt ratio ',nint((t1)/abrupttimes(i)),' X'
	endif
#endif _PLAIN_UNIX
 1	continue
	return
	end

	subroutine test(itest,x, y, z, n)
	REAL x(n),y(n),z(n)
	integer itest,n,i
	goto(1,2,3),itest
 1	continue
	do 10 i = 1, n
 10	z(i) = x(i) * y(i)
	goto 9
 2	continue
	do 20 i = 1, n
 20	z(i) = x(i) * y(i)
	goto 9
 3	continue
	do 30 i = 1, n
 30	z(i) = exp(x(i))
	goto 9
 9	continue
	end	
	
	subroutine setup(itest,x,y,n)
	REAL x(n),y(n)
	integer itest,n,i
	doubleprecision u
#ifdef __NEWVALID
	doubleprecision doublevalidrandom
#define dorand(X)	doublevalidrandom()
#else
	doubleprecision dorand
#endif

	goto(1,2,3),itest
 1	continue
	do 10 i = 1, n
	u = dorand(0)
	x(i) = sqrt(MINSUB) * 
     1	(MINNORM/MINSUB)**(u/2)
	u = dorand(0)
	y(i) = sqrt(MINSUB) * 
     1	(MINNORM/MINSUB)**(u/2)
 10	continue
	goto 9
 2	continue
	do 20 i = 1, n
	u = dorand(0)
	x(i) = MINNORM * 
     1	(MINSUB/MINNORM)**u
	u = dorand(0)
	y(i) = (MINNORM/MINSUB) * 
     1	((MINSUB/MINNORM)/MINNORM) ** u
 20	continue
	goto 9
 3	continue
	do 30 i = 1, n
 30	x(i) = -i
	goto 9
 9	continue
	end	

#ifndef __NEWVALID
	double precision function dorand(n)
	integer n
	real rand
#ifdef _AIX
	dorand = rand()
#else
	dorand = rand(n)
#endif
	end
#endif

	REAL function dummy(t)
	REAL t
	dummy = t
	end

	logical function gradual()
	REAL t,t1,t2,delta,dummy
	integer i
	t = 1
	t1=1
	t2=1
	do 1 i = 1,1000000
	t2=t1
	t1=t
	t = dummy(0.5*t)
	if (t .ge. t1) goto 2
 1	continue
 2	continue
	delta = (dummy(1.01 * t2) - t2)/t2
	if ((delta .gt. .005).and.(delta .lt. .02)) then
	gradual = .false.
	else
	gradual = .true.
	endif
	end