Mercurial > hg > octave-lyh

annotate libcruft/lapack/slamc5.f @ 8489:e76b92c7f779

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[docs] low cost => low-cost
author Brian Gough <bjg@gnu.org>
date Tue, 13 Jan 2009 00:26:46 -0500
parents 82be108cc558
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 * -- LAPACK auxiliary routine (version 3.1) --
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 * November 2006
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 * .. Scalar Arguments ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 LOGICAL IEEE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 INTEGER BETA, EMAX, EMIN, P
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 REAL RMAX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 * ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 * Purpose
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 * =======
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 * SLAMC5 attempts to compute RMAX, the largest machine floating-point
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 * number, without overflow. It assumes that EMAX + abs(EMIN) sum
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
18 * approximately to a power of 2. It will fail on machines where this
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 * EMAX = 28718). It will also fail if the value supplied for EMIN is
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 * too large (i.e. too close to zero), probably with overflow.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 * Arguments
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 * =========
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 * BETA (input) INTEGER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 * The base of floating-point arithmetic.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 * P (input) INTEGER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 * The number of base BETA digits in the mantissa of a
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 * floating-point value.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 * EMIN (input) INTEGER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 * The minimum exponent before (gradual) underflow.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 * IEEE (input) LOGICAL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 * A logical flag specifying whether or not the arithmetic
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 * system is thought to comply with the IEEE standard.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 * EMAX (output) INTEGER
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 * The largest exponent before overflow
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 * RMAX (output) REAL
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 * The largest machine floating-point number.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 * =====================================================================
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 * .. Parameters ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 REAL ZERO, ONE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 * ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 * .. Local Scalars ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 REAL OLDY, RECBAS, Y, Z
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 * ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 * .. External Functions ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 REAL SLAMC3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 EXTERNAL SLAMC3
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 * ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 * .. Intrinsic Functions ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 INTRINSIC MOD
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 * ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 * .. Executable Statements ..
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 * First compute LEXP and UEXP, two powers of 2 that bound
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 * approximately to the bound that is closest to abs(EMIN).
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 * (EMAX is the exponent of the required number RMAX).
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 LEXP = 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 EXBITS = 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 10 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 TRY = LEXP*2
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 IF( TRY.LE.( -EMIN ) ) THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 LEXP = TRY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 EXBITS = EXBITS + 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 GO TO 10
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 END IF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 IF( LEXP.EQ.-EMIN ) THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 UEXP = LEXP
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 ELSE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 UEXP = TRY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 EXBITS = EXBITS + 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 END IF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 * than or equal to EMIN. EXBITS is the number of bits needed to
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 * store the exponent.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 EXPSUM = 2*LEXP
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 ELSE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 EXPSUM = 2*UEXP
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 END IF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 * EXPSUM is the exponent range, approximately equal to
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 * EMAX - EMIN + 1 .
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 EMAX = EXPSUM + EMIN - 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 NBITS = 1 + EXBITS + P
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 * NBITS is the total number of bits needed to store a
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 * floating-point number.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
106 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 * Either there are an odd number of bits used to store a
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 * floating-point number, which is unlikely, or some bits are
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 * not used in the representation of numbers, which is possible,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 * (e.g. Cray machines) or the mantissa has an implicit bit,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 * most likely. We have to assume the last alternative.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
113 * If this is true, then we need to reduce EMAX by one because
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 * there must be some way of representing zero in an implicit-bit
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 * system. On machines like Cray, we are reducing EMAX by one
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 * unnecessarily.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 EMAX = EMAX - 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 END IF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 IF( IEEE ) THEN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 * Assume we are on an IEEE machine which reserves one exponent
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 * for infinity and NaN.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 EMAX = EMAX - 1
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 END IF
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 * Now create RMAX, the largest machine number, which should
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 * First compute 1.0 - BETA**(-P), being careful that the
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 * result is less than 1.0 .
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 RECBAS = ONE / BETA
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 Z = BETA - ONE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
137 Y = ZERO
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 DO 20 I = 1, P
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 Z = Z*RECBAS
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
140 IF( Y.LT.ONE )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 $ OLDY = Y
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 Y = SLAMC3( Y, Z )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 20 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 IF( Y.GE.ONE )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 $ Y = OLDY
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 * Now multiply by BETA**EMAX to get RMAX.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 DO 30 I = 1, EMAX
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 Y = SLAMC3( Y*BETA, ZERO )
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 30 CONTINUE
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 RMAX = Y
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 RETURN
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 * End of SLAMC5
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 *
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
158 END