Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

PHY_CPU_NumPrec.f90 @ 5451

Last change on this file since 5451 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 26.9 KB

Rev	Line
[2089]	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