Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_CPU_NumPrec.f90 @ 4031

Last change on this file since 4031 was 3990, checked in by millour, 8 years ago
An updated LMDZ5 (based on LMDZ rev 2786) to play with. EM
File size: 26.9 KB

Line
1	subroutine PHY_CPU_NumPrec
2
3	!------------------------------------------------------------------------------+
4	! Sat 30-Mar-2013 MAR \|
5	! MAR PHY_CPU_NumPrec \|
6	! \|
7	! SubRoutine relrea Computes Machine Precision \|
8	! using a LAPACK auxiliary routine \|
9	! (Originating from CONVEX -- Ruud VAN DER PAS) \|
10	! \|
11	! version 3.p.4.1 adapted by H. Gallee, Sat 30-Mar-2013 \|
12	! Last Modification by H. Gallee, Sat 30-Mar-2013 \|
13	!------------------------------------------------------------------------------+
14
15	use Mod_Real
16	use Mod_PHY____dat
17
18
19	IMPLICIT NONE
20
21	character(len=1) :: cc
22	real(kind=real8) :: eps
23	real(kind=real8) :: sf_MIN
24	real(kind=real8) :: rr_MIN
25	real(kind=real8) :: rr_MAX
26	integer :: iargum,i__MIN,i__MAX
27
28	REAL :: SLAMCH
29
30
31
32	! Machine Precision
33	! =================
34
35	cc ='E'
36	eps = SLAMCH(cc)
37	cc ='S'
38	sf_MIN = SLAMCH(cc)
39	cc ='U'
40	rr_MIN = SLAMCH(cc)
41	cc ='O'
42	rr_MAX = SLAMCH(cc)
43
44
45	write(6,600) eps,sf_MIN,rr_MIN,rr_MAX
46	600 format(' Precision relative : ',e12.4, &
47	& /,' Plus petit nombre : ',e12.4, &
48	& /,' Underflow Threshold = ',e12.4, &
49	& /,' Overflow Threshold = ',e12.4)
50
51
52
53
54	! Min and Max Arguments of Function exp(x)
55	! ========================================
56
57	ea_MIN = log(rr_MIN)
58	iargum = ea_MIN
59	i__MIN = iargum + 7
60	ea_MAX = log(rr_MAX)
61	iargum = ea_MAX
62	i__MAX = iargum - 8
63
64	write(6,601) ea_MIN,i__MIN,ea_MAX,i__MAX
65	601 format(/,' Function exp(x) : Arguments:', &
66	& /,' Minimum Value : ',e12.4, 5x,'==> (',i3,')', &
67	& /,' Maximum Value : ',e12.4, 5x,'==> (',i3,')')
68
69	ea_MIN = i__MIN
70	ea_MAX = i__MAX
71
72
73
74	return
75	end subroutine PHY_CPU_NumPrec
76
77
78	REAL FUNCTION SLAMCH( CMACH )
79
80	! -- LAPACK auxiliary routine (version 1.0) --
81	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
82	! Courant Institute, Argonne National Lab, and Rice University
83	! February 29, 1992
84
85	! .. Scalar Arguments ..
86	CHARACTER CMACH
87	! ..
88
89	! Purpose
90	! =======
91
92	! SLAMCH determines single precision machine parameters.
93
94	! Arguments
95	! =========
96
97	! CMACH (input) CHARACTER*1
98	! Specifies the value to be returned by SLAMCH:
99	! = 'E' or 'e', SLAMCH := eps
100	! = 'S' or 's , SLAMCH := sfmin
101	! = 'B' or 'b', SLAMCH := base
102	! = 'P' or 'p', SLAMCH := eps*base
103	! = 'N' or 'n', SLAMCH := t
104	! = 'R' or 'r', SLAMCH := rnd
105	! = 'M' or 'm', SLAMCH := emin
106	! = 'U' or 'u', SLAMCH := rmin
107	! = 'L' or 'l', SLAMCH := emax
108	! = 'O' or 'o', SLAMCH := rmax
109
110	! where
111
112	! eps = relative machine precision
113	! sfmin = safe minimum, such that 1/sfmin does not overflow
114	! base = base of the machine
115	! prec = eps*base
116	! t = number of (base) digits in the mantissa
117	! rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
118	! emin = minimum exponent before (gradual) underflow
119	! rmin = underflow threshold - base**(emin-1)
120	! emax = largest exponent before overflow
121	! rmax = overflow threshold - (base*emax)(1-eps)
122
123
124	! .. Parameters ..
125	REAL ONE, ZERO
126	PARAMETER ( ONE = 1.0d+0, ZERO = 0.0d+0 )
127	! ..
128	! .. Local Scalars ..
129	LOGICAL FIRST, LRND
130	INTEGER BETA, IMAX, IMIN, IT
131	REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN
132	REAL RND, SFMIN, SMALL, T
133	! ..
134	! .. External Functions ..
135	LOGICAL LSAME
136	EXTERNAL LSAME
137	! ..
138	! .. External Subroutines ..
139	EXTERNAL SLAMC2
140	! ..
141	! .. Save statement ..
142	SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN
143	SAVE EMAX, RMAX, PREC
144	! ..
145	! .. Data statements ..
146	DATA FIRST / .TRUE. /
147	! ..
148	! .. Executable Statements ..
149
150	IF( FIRST ) THEN
151	FIRST = .FALSE.
152	CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
153	BASE = BETA
154	T = IT
155	IF( LRND ) THEN
156	RND = ONE
157	EPS = ( BASE**( 1-IT ) ) / 2
158	ELSE
159	RND = ZERO
160	EPS = BASE**( 1-IT )
161	END IF
162	PREC = EPS*BASE
163	EMIN = IMIN
164	EMAX = IMAX
165	SFMIN = RMIN
166	SMALL = ONE / RMAX
167	IF( SMALL.GE.SFMIN ) THEN
168
169	! Use SMALL plus a bit, to avoid the possibility of rounding
170	! causing overflow when computing 1/sfmin.
171
172	SFMIN = SMALL*( ONE+EPS )
173	END IF
174	END IF
175
176	IF( LSAME( CMACH, 'E' ) ) THEN
177	RMACH = EPS
178	ELSE IF( LSAME( CMACH, 'S' ) ) THEN
179	RMACH = SFMIN
180	ELSE IF( LSAME( CMACH, 'B' ) ) THEN
181	RMACH = BASE
182	ELSE IF( LSAME( CMACH, 'P' ) ) THEN
183	RMACH = PREC
184	ELSE IF( LSAME( CMACH, 'N' ) ) THEN
185	RMACH = T
186	ELSE IF( LSAME( CMACH, 'R' ) ) THEN
187	RMACH = RND
188	ELSE IF( LSAME( CMACH, 'M' ) ) THEN
189	RMACH = EMIN
190	ELSE IF( LSAME( CMACH, 'U' ) ) THEN
191	RMACH = RMIN
192	ELSE IF( LSAME( CMACH, 'L' ) ) THEN
193	RMACH = EMAX
194	ELSE IF( LSAME( CMACH, 'O' ) ) THEN
195	RMACH = RMAX
196	END IF
197
198	SLAMCH = RMACH
199	RETURN
200
201	! End of SLAMCH
202
203	END
204
205	! ***********************************************************************
206
207	SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
208
209	! -- LAPACK auxiliary routine (version 1.0) --
210	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
211	! Courant Institute, Argonne National Lab, and Rice University
212	! February 29, 1992
213
214	! .. Scalar Arguments ..
215	LOGICAL IEEE1, RND
216	INTEGER BETA, T
217	! ..
218
219	! Purpose
220	! =======
221
222	! SLAMC1 determines the machine parameters given by BETA, T, RND, and
223	! IEEE1.
224
225	! Arguments
226	! =========
227
228	! BETA (output) INTEGER
229	! The base of the machine.
230
231	! T (output) INTEGER
232	! The number of ( BETA ) digits in the mantissa.
233
234	! RND (output) LOGICAL
235	! Specifies whether proper rounding ( RND = .TRUE. ) or
236	! chopping ( RND = .FALSE. ) occurs in addition. This may not
237	! be a reliable guide to the way in which the machine performs
238	! its arithmetic.
239
240	! IEEE1 (output) LOGICAL
241	! Specifies whether rounding appears to be done in the IEEE
242	! 'round to nearest' style.
243
244	! Further Details
245	! ===============
246
247	! The routine is based on the routine ENVRON by Malcolm and
248	! incorporates suggestions by Gentleman and Marovich. See
249
250	! Malcolm M. A. (1972) Algorithms to reveal properties of
251	! floating-point arithmetic. Comms. of the ACM, 15, 949-951.
252
253	! Gentleman W. M. and Marovich S. B. (1974) More on algorithms
254	! that reveal properties of floating point arithmetic units.
255	! Comms. of the ACM, 17, 276-277.
256
257
258	! .. Local Scalars ..
259	LOGICAL FIRST, LIEEE1, LRND
260	INTEGER LBETA, LT
261	REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
262	! ..
263	! .. External Functions ..
264	REAL SLAMC3
265	EXTERNAL SLAMC3
266	! ..
267	! .. Save statement ..
268	SAVE FIRST, LIEEE1, LBETA, LRND, LT
269	! ..
270	! .. Data statements ..
271	DATA FIRST / .TRUE. /
272	! ..
273	! .. Executable Statements ..
274
275	IF( FIRST ) THEN
276	FIRST = .FALSE.
277	ONE = 1
278
279	! LBETA, LIEEE1, LT and LRND are the local values of BETA,
280	! IEEE1, T and RND.
281
282	! Throughout this routine we use the Function SLAMC3 to ensure
283	! that relevant values are stored and not held in registers, or
284	! are not affected by optimizers.
285
286	! Compute a = 2.0**m with the smallest positive integer m such
287	! that
288
289	! fl( a + 1.0 ) = a.
290
291	A = 1
292	C = 1
293
294	! WHILE( C.EQ.ONE )LOOP
295	10 CONTINUE
296	IF( C.EQ.ONE ) THEN
297	A = 2*A
298	C = SLAMC3( A, ONE )
299	C = SLAMC3( C, -A )
300	GO TO 10
301	END IF
302	! END WHILE
303
304	! Now compute b = 2.0**m with the smallest positive integer m
305	! such that
306
307	! fl( a + b ) .gt. a.
308
309	B = 1
310	C = SLAMC3( A, B )
311
312	! WHILE( C.EQ.A )LOOP
313	20 CONTINUE
314	IF( C.EQ.A ) THEN
315	B = 2*B
316	C = SLAMC3( A, B )
317	GO TO 20
318	END IF
319	! END WHILE
320
321	! Now compute the base. a and c are neighbouring floating point
322	! numbers in the interval ( betat, beta( t + 1 ) ) and so
323	! their difference is beta. Adding 0.25 to c is to ensure that it
324	! is truncated to beta and not ( beta - 1 ).
325
326	QTR = ONE / 4
327	SAVEC = C
328	C = SLAMC3( C, -A )
329	LBETA = C + QTR
330
331	! Now determine whether rounding or chopping occurs, by adding a
332	! bit less than beta/2 and a bit more than beta/2 to a.
333
334	B = LBETA
335	F = SLAMC3( B / 2, -B / 100 )
336	C = SLAMC3( F, A )
337	IF( C.EQ.A ) THEN
338	LRND = .TRUE.
339	ELSE
340	LRND = .FALSE.
341	END IF
342	F = SLAMC3( B / 2, B / 100 )
343	C = SLAMC3( F, A )
344	IF( ( LRND ) .AND. ( C.EQ.A ) ) &
345	& LRND = .FALSE.
346
347	! Try and decide whether rounding is done in the IEEE 'round to
348	! nearest' style. B/2 is half a unit in the last place of the two
349	! numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
350	! zero, and SAVEC is odd. Thus adding B/2 to A should not change
351	! A, but adding B/2 to SAVEC should change SAVEC.
352
353	T1 = SLAMC3( B / 2, A )
354	T2 = SLAMC3( B / 2, SAVEC )
355	LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
356
357	! Now find the mantissa, t. It should be the integer part of
358	! log to the base beta of a, however it is safer to determine t
359	! by powering. So we find t as the smallest positive integer for
360	! which
361
362	! fl( beta**t + 1.0 ) = 1.0.
363
364	LT = 0
365	A = 1
366	C = 1
367
368	! WHILE( C.EQ.ONE )LOOP
369	30 CONTINUE
370	IF( C.EQ.ONE ) THEN
371	LT = LT + 1
372	A = A*LBETA
373	C = SLAMC3( A, ONE )
374	C = SLAMC3( C, -A )
375	GO TO 30
376	END IF
377	! END WHILE
378
379	END IF
380
381	BETA = LBETA
382	T = LT
383	RND = LRND
384	IEEE1 = LIEEE1
385	RETURN
386
387	! End of SLAMC1
388
389	END
390
391	! ***********************************************************************
392
393	SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
394
395	! -- LAPACK auxiliary routine (version 1.0) --
396	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
397	! Courant Institute, Argonne National Lab, and Rice University
398	! February 29, 1992
399
400	! .. Scalar Arguments ..
401	LOGICAL RND
402	INTEGER BETA, EMAX, EMIN, T
403	REAL EPS, RMAX, RMIN
404	! ..
405
406	! Purpose
407	! =======
408
409	! SLAMC2 determines the machine parameters specified in its argument
410	! list.
411
412	! Arguments
413	! =========
414
415	! BETA (output) INTEGER
416	! The base of the machine.
417
418	! T (output) INTEGER
419	! The number of ( BETA ) digits in the mantissa.
420
421	! RND (output) LOGICAL
422	! Specifies whether proper rounding ( RND = .TRUE. ) or
423	! chopping ( RND = .FALSE. ) occurs in addition. This may not
424	! be a reliable guide to the way in which the machine performs
425	! its arithmetic.
426
427	! EPS (output) REAL
428	! The smallest positive number such that
429
430	! fl( 1.0 - EPS ) .LT. 1.0,
431
432	! where fl denotes the computed value.
433
434	! EMIN (output) INTEGER
435	! The minimum exponent before (gradual) underflow occurs.
436
437	! RMIN (output) REAL
438	! The smallest normalized number for the machine, given by
439	! BASE**( EMIN - 1 ), where BASE is the floating point value
440	! of BETA.
441
442	! EMAX (output) INTEGER
443	! The maximum exponent before overflow occurs.
444
445	! RMAX (output) REAL
446	! The largest positive number for the machine, given by
447	! BASE*EMAX ( 1 - EPS ), where BASE is the floating point
448	! value of BETA.
449
450	! Further Details
451	! ===============
452
453	! The computation of EPS is based on a routine PARANOIA by
454	! W. Kahan of the University of California at Berkeley.
455
456
457	! .. Local Scalars ..
458	LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
459	INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT
460	INTEGER NGNMIN, NGPMIN
461	REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE
462	REAL SIXTH, SMALL, THIRD, TWO, ZERO
463	! ..
464	! .. External Functions ..
465	REAL SLAMC3
466	EXTERNAL SLAMC3
467	! ..
468	! .. External Subroutines ..
469	EXTERNAL SLAMC1, SLAMC4, SLAMC5
470	! ..
471	! .. Intrinsic Functions ..
472	INTRINSIC ABS, MAX, MIN
473	! ..
474	! .. Save statement ..
475	SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX
476	SAVE LRMIN, LT
477	! ..
478	! .. Data statements ..
479	DATA FIRST / .TRUE. / , IWARN / .FALSE. /
480	! ..
481	! .. Executable Statements ..
482
483	IF( FIRST ) THEN
484	FIRST = .FALSE.
485	ZERO = 0
486	ONE = 1
487	TWO = 2
488
489	! LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
490	! BETA, T, RND, EPS, EMIN and RMIN.
491
492	! Throughout this routine we use the Function SLAMC3 to ensure
493	! that relevant values are stored and not held in registers, or
494	! are not affected by optimizers.
495
496	! SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
497
498	CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
499
500	! Start to find EPS.
501
502	B = LBETA
503	A = B**( -LT )
504	LEPS = A
505
506	! Try some tricks to see whether or not this is the correct EPS.
507
508	B = TWO / 3
509	HALF = ONE / 2
510	SIXTH = SLAMC3( B, -HALF )
511	THIRD = SLAMC3( SIXTH, SIXTH )
512	B = SLAMC3( THIRD, -HALF )
513	B = SLAMC3( B, SIXTH )
514	B = ABS( B )
515	IF( B.LT.LEPS ) &
516	& B = LEPS
517
518	LEPS = 1
519
520	! WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
521	10 CONTINUE
522	IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
523	LEPS = B
524	C = SLAMC3( HALFLEPS, ( TWO5 )( LEPS**2 ) )
525	C = SLAMC3( HALF, -C )
526	B = SLAMC3( HALF, C )
527	C = SLAMC3( HALF, -B )
528	B = SLAMC3( HALF, C )
529	GO TO 10
530	END IF
531	! END WHILE
532
533	IF( A.LT.LEPS ) &
534	& LEPS = A
535
536	! Computation of EPS complete.
537
538	! Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
539	! Keep dividing A by BETA until (gradual) underflow occurs. This
540	! is detected when we cannot recover the previous A.
541
542	RBASE = ONE / LBETA
543	SMALL = ONE
544	DO 20 I = 1, 3
545	SMALL = SLAMC3( SMALL*RBASE, ZERO )
546	20 CONTINUE
547	A = SLAMC3( ONE, SMALL )
548	CALL SLAMC4( NGPMIN, ONE, LBETA )
549	CALL SLAMC4( NGNMIN, -ONE, LBETA )
550	CALL SLAMC4( GPMIN, A, LBETA )
551	CALL SLAMC4( GNMIN, -A, LBETA )
552	IEEE = .FALSE.
553
554	IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
555	IF( NGPMIN.EQ.GPMIN ) THEN
556	LEMIN = NGPMIN
557	! ( Non twos-complement machines, no gradual underflow;
558	! e.g., VAX )
559	ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
560	LEMIN = NGPMIN - 1 + LT
561	IEEE = .TRUE.
562	! ( Non twos-complement machines, with gradual underflow;
563	! e.g., IEEE standard followers )
564	ELSE
565	LEMIN = MIN( NGPMIN, GPMIN )
566	! ( A guess; no known machine )
567	IWARN = .TRUE.
568	END IF
569
570	ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
571	IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
572	LEMIN = MAX( NGPMIN, NGNMIN )
573	! ( Twos-complement machines, no gradual underflow;
574	! e.g., CYBER 205 )
575	ELSE
576	LEMIN = MIN( NGPMIN, NGNMIN )
577	! ( A guess; no known machine )
578	IWARN = .TRUE.
579	END IF
580
581	ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. &
582	& ( GPMIN.EQ.GNMIN ) ) THEN
583	IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
584	LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
585	! ( Twos-complement machines with gradual underflow;
586	! no known machine )
587	ELSE
588	LEMIN = MIN( NGPMIN, NGNMIN )
589	! ( A guess; no known machine )
590	IWARN = .TRUE.
591	END IF
592
593	ELSE
594	LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
595	! ( A guess; no known machine )
596	IWARN = .TRUE.
597	END IF
598	! **
599	! Comment out this if block if EMIN is ok
600	IF( IWARN ) THEN
601	FIRST = .TRUE.
602	WRITE( 6, FMT = 9999 )LEMIN
603	END IF
604	! **
605
606	! Assume IEEE arithmetic if we found denormalised numbers above,
607	! or if arithmetic seems to round in the IEEE style, determined
608	! in routine SLAMC1. A true IEEE machine should have both things
609	! true; however, faulty machines may have one or the other.
610
611	IEEE = IEEE .OR. LIEEE1
612
613	! Compute RMIN by successive division by BETA. We could compute
614	! RMIN as BASE**( EMIN - 1 ), but some machines underflow during
615	! this computation.
616
617	LRMIN = 1
618	DO 30 I = 1, 1 - LEMIN
619	LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
620	30 CONTINUE
621
622	! Finally, call SLAMC5 to compute EMAX and RMAX.
623
624	CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
625	END IF
626
627	BETA = LBETA
628	T = LT
629	RND = LRND
630	EPS = LEPS
631	EMIN = LEMIN
632	RMIN = LRMIN
633	EMAX = LEMAX
634	RMAX = LRMAX
635
636	RETURN
637
638	9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', &
639	& ' EMIN = ', I8, / &
640	& ' If, after inspection, the value EMIN looks', &
641	& ' acceptable please comment out ', &
642	& / ' the IF block as marked within the code of routine', &
643	& ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
644
645	! End of SLAMC2
646
647	END
648
649	! ***********************************************************************
650
651	REAL FUNCTION SLAMC3( A, B )
652
653	! -- LAPACK auxiliary routine (version 1.0) --
654	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
655	! Courant Institute, Argonne National Lab, and Rice University
656	! February 29, 1992
657
658	! .. Scalar Arguments ..
659	REAL A, B
660	! ..
661
662	! Purpose
663	! =======
664
665	! SLAMC3 is intended to force A and B to be stored prior to doing
666	! the addition of A and B , for use in situations where optimizers
667	! might hold one of these in a register.
668
669	! Arguments
670	! =========
671
672	! A, B (input) REAL
673	! The values A and B.
674
675
676	! .. Executable Statements ..
677
678	SLAMC3 = A + B
679
680	RETURN
681
682	! End of SLAMC3
683
684	END
685
686	! ***********************************************************************
687
688	SUBROUTINE SLAMC4( EMIN, START, BASE )
689
690	! -- LAPACK auxiliary routine (version 1.0) --
691	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
692	! Courant Institute, Argonne National Lab, and Rice University
693	! February 29, 1992
694
695	! .. Scalar Arguments ..
696	INTEGER BASE, EMIN
697	REAL START
698	! ..
699
700	! Purpose
701	! =======
702
703	! SLAMC4 is a service routine for SLAMC2.
704
705	! Arguments
706	! =========
707
708	! EMIN (output) EMIN
709	! The minimum exponent before (gradual) underflow, computed by
710	! setting A = START and dividing by BASE until the previous A
711	! can not be recovered.
712
713	! START (input) REAL
714	! The starting point for determining EMIN.
715
716	! BASE (input) INTEGER
717	! The base of the machine.
718
719
720	! .. Local Scalars ..
721	INTEGER I
722	REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
723	! ..
724	! .. External Functions ..
725	REAL SLAMC3
726	EXTERNAL SLAMC3
727	! ..
728	! .. Executable Statements ..
729
730	A = START
731	ONE = 1
732	RBASE = ONE / BASE
733	ZERO = 0
734	EMIN = 1
735	B1 = SLAMC3( A*RBASE, ZERO )
736	C1 = A
737	C2 = A
738	D1 = A
739	D2 = A
740	! WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
741	! $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
742	10 CONTINUE
743	IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. &
744	& ( D2.EQ.A ) ) THEN
745	EMIN = EMIN - 1
746	A = B1
747	B1 = SLAMC3( A / BASE, ZERO )
748	C1 = SLAMC3( B1*BASE, ZERO )
749	D1 = ZERO
750	DO 20 I = 1, BASE
751	D1 = D1 + B1
752	20 CONTINUE
753	B2 = SLAMC3( A*RBASE, ZERO )
754	C2 = SLAMC3( B2 / RBASE, ZERO )
755	D2 = ZERO
756	DO 30 I = 1, BASE
757	D2 = D2 + B2
758	30 CONTINUE
759	GO TO 10
760	END IF
761	! END WHILE
762
763	RETURN
764
765	! End of SLAMC4
766
767	END
768
769	! ***********************************************************************
770
771	SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
772
773	! -- LAPACK auxiliary routine (version 1.0) --
774	! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
775	! Courant Institute, Argonne National Lab, and Rice University
776	! February 29, 1992
777
778	! .. Scalar Arguments ..
779	LOGICAL IEEE
780	INTEGER BETA, EMAX, EMIN, P
781	REAL RMAX
782	! ..
783
784	! Purpose
785	! =======
786
787	! SLAMC5 attempts to compute RMAX, the largest machine floating-point
788	! number, without overflow. It assumes that EMAX + abs(EMIN) sum
789	! approximately to a power of 2. It will fail on machines where this
790	! assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
791	! EMAX = 28718). It will also fail if the value supplied for EMIN is
792	! too large (i.e. too close to zero), probably with overflow.
793
794	! Arguments
795	! =========
796
797	! BETA (input) INTEGER
798	! The base of floating-point arithmetic.
799
800	! P (input) INTEGER
801	! The number of base BETA digits in the mantissa of a
802	! floating-point value.
803
804	! EMIN (input) INTEGER
805	! The minimum exponent before (gradual) underflow.
806
807	! IEEE (input) LOGICAL
808	! A logical flag specifying whether or not the arithmetic
809	! system is thought to comply with the IEEE standard.
810
811	! EMAX (output) INTEGER
812	! The largest exponent before overflow
813
814	! RMAX (output) REAL
815	! The largest machine floating-point number.
816
817
818	! .. Parameters ..
819	REAL ZERO, ONE
820	PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
821	! ..
822	! .. Local Scalars ..
823	INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
824	REAL OLDY, RECBAS, Y, Z
825	! ..
826	! .. External Functions ..
827	REAL SLAMC3
828	EXTERNAL SLAMC3
829	! ..
830	! .. Intrinsic Functions ..
831	INTRINSIC MOD
832	! ..
833	! .. Executable Statements ..
834
835	! First compute LEXP and UEXP, two powers of 2 that bound
836	! abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
837	! approximately to the bound that is closest to abs(EMIN).
838	! (EMAX is the exponent of the required number RMAX).
839
840	LEXP = 1
841	EXBITS = 1
842	10 CONTINUE
843	TRY = LEXP*2
844	IF( TRY.LE.( -EMIN ) ) THEN
845	LEXP = TRY
846	EXBITS = EXBITS + 1
847	GO TO 10
848	END IF
849	IF( LEXP.EQ.-EMIN ) THEN
850	UEXP = LEXP
851	ELSE
852	UEXP = TRY
853	EXBITS = EXBITS + 1
854	END IF
855
856	! Now -LEXP is less than or equal to EMIN, and -UEXP is greater
857	! than or equal to EMIN. EXBITS is the number of bits needed to
858	! store the exponent.
859
860	IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
861	EXPSUM = 2*LEXP
862	ELSE
863	EXPSUM = 2*UEXP
864	END IF
865
866	! EXPSUM is the exponent range, approximately equal to
867	! EMAX - EMIN + 1 .
868
869	EMAX = EXPSUM + EMIN - 1
870	NBITS = 1 + EXBITS + P
871
872	! NBITS is the total number of bits needed to store a
873	! floating-point number.
874
875	IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
876
877	! Either there are an odd number of bits used to store a
878	! floating-point number, which is unlikely, or some bits are
879	! not used in the representation of numbers, which is possible,
880	! (e.g. Cray machines) or the mantissa has an implicit bit,
881	! (e.g. IEEE machines, Dec Vax machines), which is perhaps the
882	! most likely. We have to assume the last alternative.
883	! If this is true, then we need to reduce EMAX by one because
884	! there must be some way of representing zero in an implicit-bit
885	! system. On machines like Cray, we are reducing EMAX by one
886	! unnecessarily.
887
888	EMAX = EMAX - 1
889	END IF
890
891	IF( IEEE ) THEN
892
893	! Assume we are on an IEEE machine which reserves one exponent
894	! for infinity and NaN.
895
896	EMAX = EMAX - 1
897	END IF
898
899	! Now create RMAX, the largest machine number, which should
900	! be equal to (1.0 - BETA*(-P)) BETA**EMAX .
901
902	! First compute 1.0 - BETA**(-P), being careful that the
903	! result is less than 1.0 .
904
905	RECBAS = ONE / BETA
906	Z = BETA - ONE
907	Y = ZERO
908	DO 20 I = 1, P
909	Z = Z*RECBAS
910	IF( Y.LT.ONE ) &
911	& OLDY = Y
912	Y = SLAMC3( Y, Z )
913	20 CONTINUE
914	IF( Y.GE.ONE ) &
915	& Y = OLDY
916
917	! Now multiply by BETA**EMAX to get RMAX.
918
919	DO 30 I = 1, EMAX
920	Y = SLAMC3( Y*BETA, ZERO )
921	30 CONTINUE
922
923	RMAX = Y
924	RETURN
925
926	! End of SLAMC5
927
928	END
929
930
931	LOGICAL FUNCTION LSAME(CH1,CH2)
932
933	!------------------------------------------+
934	! cfr. Hubert Gallee, 30 octobre 92 \|
935	!------------------------------------------+
936
937	CHARACTER CH1,CH2
938	IF (CH1.EQ.CH2) THEN
939	LSAME = .TRUE.
940	ELSE
941	LSAME = .FALSE.
942	END IF
943	RETURN
944	END
945

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format