source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phymar/PHY_CPU_NumPrec.f90 @ 3405

Last change on this file since 3405 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
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  ( beta**t, 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( HALF*LEPS, ( TWO**5 )*( 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.