001:       REAL             FUNCTION SLAMCH( CMACH )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          CMACH
009: *     ..
010: *
011: *  Purpose
012: *  =======
013: *
014: *  SLAMCH determines single precision machine parameters.
015: *
016: *  Arguments
017: *  =========
018: *
019: *  CMACH   (input) CHARACTER*1
020: *          Specifies the value to be returned by SLAMCH:
021: *          = 'E' or 'e',   SLAMCH := eps
022: *          = 'S' or 's ,   SLAMCH := sfmin
023: *          = 'B' or 'b',   SLAMCH := base
024: *          = 'P' or 'p',   SLAMCH := eps*base
025: *          = 'N' or 'n',   SLAMCH := t
026: *          = 'R' or 'r',   SLAMCH := rnd
027: *          = 'M' or 'm',   SLAMCH := emin
028: *          = 'U' or 'u',   SLAMCH := rmin
029: *          = 'L' or 'l',   SLAMCH := emax
030: *          = 'O' or 'o',   SLAMCH := rmax
031: *
032: *          where
033: *
034: *          eps   = relative machine precision
035: *          sfmin = safe minimum, such that 1/sfmin does not overflow
036: *          base  = base of the machine
037: *          prec  = eps*base
038: *          t     = number of (base) digits in the mantissa
039: *          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
040: *          emin  = minimum exponent before (gradual) underflow
041: *          rmin  = underflow threshold - base**(emin-1)
042: *          emax  = largest exponent before overflow
043: *          rmax  = overflow threshold  - (base**emax)*(1-eps)
044: *
045: * =====================================================================
046: *
047: *     .. Parameters ..
048:       REAL               ONE, ZERO
049:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
050: *     ..
051: *     .. Local Scalars ..
052:       LOGICAL            FIRST, LRND
053:       INTEGER            BETA, IMAX, IMIN, IT
054:       REAL               BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
055:      $                   RND, SFMIN, SMALL, T
056: *     ..
057: *     .. External Functions ..
058:       LOGICAL            LSAME
059:       EXTERNAL           LSAME
060: *     ..
061: *     .. External Subroutines ..
062:       EXTERNAL           SLAMC2
063: *     ..
064: *     .. Save statement ..
065:       SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
066:      $                   EMAX, RMAX, PREC
067: *     ..
068: *     .. Data statements ..
069:       DATA               FIRST / .TRUE. /
070: *     ..
071: *     .. Executable Statements ..
072: *
073:       IF( FIRST ) THEN
074:          CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
075:          BASE = BETA
076:          T = IT
077:          IF( LRND ) THEN
078:             RND = ONE
079:             EPS = ( BASE**( 1-IT ) ) / 2
080:          ELSE
081:             RND = ZERO
082:             EPS = BASE**( 1-IT )
083:          END IF
084:          PREC = EPS*BASE
085:          EMIN = IMIN
086:          EMAX = IMAX
087:          SFMIN = RMIN
088:          SMALL = ONE / RMAX
089:          IF( SMALL.GE.SFMIN ) THEN
090: *
091: *           Use SMALL plus a bit, to avoid the possibility of rounding
092: *           causing overflow when computing  1/sfmin.
093: *
094:             SFMIN = SMALL*( ONE+EPS )
095:          END IF
096:       END IF
097: *
098:       IF( LSAME( CMACH, 'E' ) ) THEN
099:          RMACH = EPS
100:       ELSE IF( LSAME( CMACH, 'S' ) ) THEN
101:          RMACH = SFMIN
102:       ELSE IF( LSAME( CMACH, 'B' ) ) THEN
103:          RMACH = BASE
104:       ELSE IF( LSAME( CMACH, 'P' ) ) THEN
105:          RMACH = PREC
106:       ELSE IF( LSAME( CMACH, 'N' ) ) THEN
107:          RMACH = T
108:       ELSE IF( LSAME( CMACH, 'R' ) ) THEN
109:          RMACH = RND
110:       ELSE IF( LSAME( CMACH, 'M' ) ) THEN
111:          RMACH = EMIN
112:       ELSE IF( LSAME( CMACH, 'U' ) ) THEN
113:          RMACH = RMIN
114:       ELSE IF( LSAME( CMACH, 'L' ) ) THEN
115:          RMACH = EMAX
116:       ELSE IF( LSAME( CMACH, 'O' ) ) THEN
117:          RMACH = RMAX
118:       END IF
119: *
120:       SLAMCH = RMACH
121:       FIRST  = .FALSE.
122:       RETURN
123: *
124: *     End of SLAMCH
125: *
126:       END
127: *
128: ************************************************************************
129: *
130:       SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
131: *
132: *  -- LAPACK auxiliary routine (version 3.2) --
133: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
134: *     November 2006
135: *
136: *     .. Scalar Arguments ..
137:       LOGICAL            IEEE1, RND
138:       INTEGER            BETA, T
139: *     ..
140: *
141: *  Purpose
142: *  =======
143: *
144: *  SLAMC1 determines the machine parameters given by BETA, T, RND, and
145: *  IEEE1.
146: *
147: *  Arguments
148: *  =========
149: *
150: *  BETA    (output) INTEGER
151: *          The base of the machine.
152: *
153: *  T       (output) INTEGER
154: *          The number of ( BETA ) digits in the mantissa.
155: *
156: *  RND     (output) LOGICAL
157: *          Specifies whether proper rounding  ( RND = .TRUE. )  or
158: *          chopping  ( RND = .FALSE. )  occurs in addition. This may not
159: *          be a reliable guide to the way in which the machine performs
160: *          its arithmetic.
161: *
162: *  IEEE1   (output) LOGICAL
163: *          Specifies whether rounding appears to be done in the IEEE
164: *          'round to nearest' style.
165: *
166: *  Further Details
167: *  ===============
168: *
169: *  The routine is based on the routine  ENVRON  by Malcolm and
170: *  incorporates suggestions by Gentleman and Marovich. See
171: *
172: *     Malcolm M. A. (1972) Algorithms to reveal properties of
173: *        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
174: *
175: *     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
176: *        that reveal properties of floating point arithmetic units.
177: *        Comms. of the ACM, 17, 276-277.
178: *
179: * =====================================================================
180: *
181: *     .. Local Scalars ..
182:       LOGICAL            FIRST, LIEEE1, LRND
183:       INTEGER            LBETA, LT
184:       REAL               A, B, C, F, ONE, QTR, SAVEC, T1, T2
185: *     ..
186: *     .. External Functions ..
187:       REAL               SLAMC3
188:       EXTERNAL           SLAMC3
189: *     ..
190: *     .. Save statement ..
191:       SAVE               FIRST, LIEEE1, LBETA, LRND, LT
192: *     ..
193: *     .. Data statements ..
194:       DATA               FIRST / .TRUE. /
195: *     ..
196: *     .. Executable Statements ..
197: *
198:       IF( FIRST ) THEN
199:          ONE = 1
200: *
201: *        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
202: *        IEEE1, T and RND.
203: *
204: *        Throughout this routine  we use the function  SLAMC3  to ensure
205: *        that relevant values are  stored and not held in registers,  or
206: *        are not affected by optimizers.
207: *
208: *        Compute  a = 2.0**m  with the  smallest positive integer m such
209: *        that
210: *
211: *           fl( a + 1.0 ) = a.
212: *
213:          A = 1
214:          C = 1
215: *
216: *+       WHILE( C.EQ.ONE )LOOP
217:    10    CONTINUE
218:          IF( C.EQ.ONE ) THEN
219:             A = 2*A
220:             C = SLAMC3( A, ONE )
221:             C = SLAMC3( C, -A )
222:             GO TO 10
223:          END IF
224: *+       END WHILE
225: *
226: *        Now compute  b = 2.0**m  with the smallest positive integer m
227: *        such that
228: *
229: *           fl( a + b ) .gt. a.
230: *
231:          B = 1
232:          C = SLAMC3( A, B )
233: *
234: *+       WHILE( C.EQ.A )LOOP
235:    20    CONTINUE
236:          IF( C.EQ.A ) THEN
237:             B = 2*B
238:             C = SLAMC3( A, B )
239:             GO TO 20
240:          END IF
241: *+       END WHILE
242: *
243: *        Now compute the base.  a and c  are neighbouring floating point
244: *        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
245: *        their difference is beta. Adding 0.25 to c is to ensure that it
246: *        is truncated to beta and not ( beta - 1 ).
247: *
248:          QTR = ONE / 4
249:          SAVEC = C
250:          C = SLAMC3( C, -A )
251:          LBETA = C + QTR
252: *
253: *        Now determine whether rounding or chopping occurs,  by adding a
254: *        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
255: *
256:          B = LBETA
257:          F = SLAMC3( B / 2, -B / 100 )
258:          C = SLAMC3( F, A )
259:          IF( C.EQ.A ) THEN
260:             LRND = .TRUE.
261:          ELSE
262:             LRND = .FALSE.
263:          END IF
264:          F = SLAMC3( B / 2, B / 100 )
265:          C = SLAMC3( F, A )
266:          IF( ( LRND ) .AND. ( C.EQ.A ) )
267:      $      LRND = .FALSE.
268: *
269: *        Try and decide whether rounding is done in the  IEEE  'round to
270: *        nearest' style. B/2 is half a unit in the last place of the two
271: *        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
272: *        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
273: *        A, but adding B/2 to SAVEC should change SAVEC.
274: *
275:          T1 = SLAMC3( B / 2, A )
276:          T2 = SLAMC3( B / 2, SAVEC )
277:          LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
278: *
279: *        Now find  the  mantissa, t.  It should  be the  integer part of
280: *        log to the base beta of a,  however it is safer to determine  t
281: *        by powering.  So we find t as the smallest positive integer for
282: *        which
283: *
284: *           fl( beta**t + 1.0 ) = 1.0.
285: *
286:          LT = 0
287:          A = 1
288:          C = 1
289: *
290: *+       WHILE( C.EQ.ONE )LOOP
291:    30    CONTINUE
292:          IF( C.EQ.ONE ) THEN
293:             LT = LT + 1
294:             A = A*LBETA
295:             C = SLAMC3( A, ONE )
296:             C = SLAMC3( C, -A )
297:             GO TO 30
298:          END IF
299: *+       END WHILE
300: *
301:       END IF
302: *
303:       BETA = LBETA
304:       T = LT
305:       RND = LRND
306:       IEEE1 = LIEEE1
307:       FIRST = .FALSE.
308:       RETURN
309: *
310: *     End of SLAMC1
311: *
312:       END
313: *
314: ************************************************************************
315: *
316:       SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
317: *
318: *  -- LAPACK auxiliary routine (version 3.2) --
319: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
320: *     November 2006
321: *
322: *     .. Scalar Arguments ..
323:       LOGICAL            RND
324:       INTEGER            BETA, EMAX, EMIN, T
325:       REAL               EPS, RMAX, RMIN
326: *     ..
327: *
328: *  Purpose
329: *  =======
330: *
331: *  SLAMC2 determines the machine parameters specified in its argument
332: *  list.
333: *
334: *  Arguments
335: *  =========
336: *
337: *  BETA    (output) INTEGER
338: *          The base of the machine.
339: *
340: *  T       (output) INTEGER
341: *          The number of ( BETA ) digits in the mantissa.
342: *
343: *  RND     (output) LOGICAL
344: *          Specifies whether proper rounding  ( RND = .TRUE. )  or
345: *          chopping  ( RND = .FALSE. )  occurs in addition. This may not
346: *          be a reliable guide to the way in which the machine performs
347: *          its arithmetic.
348: *
349: *  EPS     (output) REAL
350: *          The smallest positive number such that
351: *
352: *             fl( 1.0 - EPS ) .LT. 1.0,
353: *
354: *          where fl denotes the computed value.
355: *
356: *  EMIN    (output) INTEGER
357: *          The minimum exponent before (gradual) underflow occurs.
358: *
359: *  RMIN    (output) REAL
360: *          The smallest normalized number for the machine, given by
361: *          BASE**( EMIN - 1 ), where  BASE  is the floating point value
362: *          of BETA.
363: *
364: *  EMAX    (output) INTEGER
365: *          The maximum exponent before overflow occurs.
366: *
367: *  RMAX    (output) REAL
368: *          The largest positive number for the machine, given by
369: *          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
370: *          value of BETA.
371: *
372: *  Further Details
373: *  ===============
374: *
375: *  The computation of  EPS  is based on a routine PARANOIA by
376: *  W. Kahan of the University of California at Berkeley.
377: *
378: * =====================================================================
379: *
380: *     .. Local Scalars ..
381:       LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
382:       INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
383:      $                   NGNMIN, NGPMIN
384:       REAL               A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
385:      $                   SIXTH, SMALL, THIRD, TWO, ZERO
386: *     ..
387: *     .. External Functions ..
388:       REAL               SLAMC3
389:       EXTERNAL           SLAMC3
390: *     ..
391: *     .. External Subroutines ..
392:       EXTERNAL           SLAMC1, SLAMC4, SLAMC5
393: *     ..
394: *     .. Intrinsic Functions ..
395:       INTRINSIC          ABS, MAX, MIN
396: *     ..
397: *     .. Save statement ..
398:       SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
399:      $                   LRMIN, LT
400: *     ..
401: *     .. Data statements ..
402:       DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
403: *     ..
404: *     .. Executable Statements ..
405: *
406:       IF( FIRST ) THEN
407:          ZERO = 0
408:          ONE = 1
409:          TWO = 2
410: *
411: *        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
412: *        BETA, T, RND, EPS, EMIN and RMIN.
413: *
414: *        Throughout this routine  we use the function  SLAMC3  to ensure
415: *        that relevant values are stored  and not held in registers,  or
416: *        are not affected by optimizers.
417: *
418: *        SLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
419: *
420:          CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
421: *
422: *        Start to find EPS.
423: *
424:          B = LBETA
425:          A = B**( -LT )
426:          LEPS = A
427: *
428: *        Try some tricks to see whether or not this is the correct  EPS.
429: *
430:          B = TWO / 3
431:          HALF = ONE / 2
432:          SIXTH = SLAMC3( B, -HALF )
433:          THIRD = SLAMC3( SIXTH, SIXTH )
434:          B = SLAMC3( THIRD, -HALF )
435:          B = SLAMC3( B, SIXTH )
436:          B = ABS( B )
437:          IF( B.LT.LEPS )
438:      $      B = LEPS
439: *
440:          LEPS = 1
441: *
442: *+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
443:    10    CONTINUE
444:          IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
445:             LEPS = B
446:             C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
447:             C = SLAMC3( HALF, -C )
448:             B = SLAMC3( HALF, C )
449:             C = SLAMC3( HALF, -B )
450:             B = SLAMC3( HALF, C )
451:             GO TO 10
452:          END IF
453: *+       END WHILE
454: *
455:          IF( A.LT.LEPS )
456:      $      LEPS = A
457: *
458: *        Computation of EPS complete.
459: *
460: *        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
461: *        Keep dividing  A by BETA until (gradual) underflow occurs. This
462: *        is detected when we cannot recover the previous A.
463: *
464:          RBASE = ONE / LBETA
465:          SMALL = ONE
466:          DO 20 I = 1, 3
467:             SMALL = SLAMC3( SMALL*RBASE, ZERO )
468:    20    CONTINUE
469:          A = SLAMC3( ONE, SMALL )
470:          CALL SLAMC4( NGPMIN, ONE, LBETA )
471:          CALL SLAMC4( NGNMIN, -ONE, LBETA )
472:          CALL SLAMC4( GPMIN, A, LBETA )
473:          CALL SLAMC4( GNMIN, -A, LBETA )
474:          IEEE = .FALSE.
475: *
476:          IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
477:             IF( NGPMIN.EQ.GPMIN ) THEN
478:                LEMIN = NGPMIN
479: *            ( Non twos-complement machines, no gradual underflow;
480: *              e.g.,  VAX )
481:             ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
482:                LEMIN = NGPMIN - 1 + LT
483:                IEEE = .TRUE.
484: *            ( Non twos-complement machines, with gradual underflow;
485: *              e.g., IEEE standard followers )
486:             ELSE
487:                LEMIN = MIN( NGPMIN, GPMIN )
488: *            ( A guess; no known machine )
489:                IWARN = .TRUE.
490:             END IF
491: *
492:          ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
493:             IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
494:                LEMIN = MAX( NGPMIN, NGNMIN )
495: *            ( Twos-complement machines, no gradual underflow;
496: *              e.g., CYBER 205 )
497:             ELSE
498:                LEMIN = MIN( NGPMIN, NGNMIN )
499: *            ( A guess; no known machine )
500:                IWARN = .TRUE.
501:             END IF
502: *
503:          ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
504:      $            ( GPMIN.EQ.GNMIN ) ) THEN
505:             IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
506:                LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
507: *            ( Twos-complement machines with gradual underflow;
508: *              no known machine )
509:             ELSE
510:                LEMIN = MIN( NGPMIN, NGNMIN )
511: *            ( A guess; no known machine )
512:                IWARN = .TRUE.
513:             END IF
514: *
515:          ELSE
516:             LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
517: *         ( A guess; no known machine )
518:             IWARN = .TRUE.
519:          END IF
520:          FIRST = .FALSE.
521: ***
522: * Comment out this if block if EMIN is ok
523:          IF( IWARN ) THEN
524:             FIRST = .TRUE.
525:             WRITE( 6, FMT = 9999 )LEMIN
526:          END IF
527: ***
528: *
529: *        Assume IEEE arithmetic if we found denormalised  numbers above,
530: *        or if arithmetic seems to round in the  IEEE style,  determined
531: *        in routine SLAMC1. A true IEEE machine should have both  things
532: *        true; however, faulty machines may have one or the other.
533: *
534:          IEEE = IEEE .OR. LIEEE1
535: *
536: *        Compute  RMIN by successive division by  BETA. We could compute
537: *        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
538: *        this computation.
539: *
540:          LRMIN = 1
541:          DO 30 I = 1, 1 - LEMIN
542:             LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
543:    30    CONTINUE
544: *
545: *        Finally, call SLAMC5 to compute EMAX and RMAX.
546: *
547:          CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
548:       END IF
549: *
550:       BETA = LBETA
551:       T = LT
552:       RND = LRND
553:       EPS = LEPS
554:       EMIN = LEMIN
555:       RMIN = LRMIN
556:       EMAX = LEMAX
557:       RMAX = LRMAX
558: *
559:       RETURN
560: *
561:  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
562:      $      '  EMIN = ', I8, /
563:      $      ' If, after inspection, the value EMIN looks',
564:      $      ' acceptable please comment out ',
565:      $      / ' the IF block as marked within the code of routine',
566:      $      ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
567: *
568: *     End of SLAMC2
569: *
570:       END
571: *
572: ************************************************************************
573: *
574:       REAL             FUNCTION SLAMC3( A, B )
575: *
576: *  -- LAPACK auxiliary routine (version 3.2) --
577: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
578: *     November 2006
579: *
580: *     .. Scalar Arguments ..
581:       REAL               A, B
582: *     ..
583: *
584: *  Purpose
585: *  =======
586: *
587: *  SLAMC3  is intended to force  A  and  B  to be stored prior to doing
588: *  the addition of  A  and  B ,  for use in situations where optimizers
589: *  might hold one of these in a register.
590: *
591: *  Arguments
592: *  =========
593: *
594: *  A       (input) REAL
595: *  B       (input) REAL
596: *          The values A and B.
597: *
598: * =====================================================================
599: *
600: *     .. Executable Statements ..
601: *
602:       SLAMC3 = A + B
603: *
604:       RETURN
605: *
606: *     End of SLAMC3
607: *
608:       END
609: *
610: ************************************************************************
611: *
612:       SUBROUTINE SLAMC4( EMIN, START, BASE )
613: *
614: *  -- LAPACK auxiliary routine (version 3.2) --
615: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
616: *     November 2006
617: *
618: *     .. Scalar Arguments ..
619:       INTEGER            BASE
620:       INTEGER            EMIN
621:       REAL               START
622: *     ..
623: *
624: *  Purpose
625: *  =======
626: *
627: *  SLAMC4 is a service routine for SLAMC2.
628: *
629: *  Arguments
630: *  =========
631: *
632: *  EMIN    (output) INTEGER 
633: *          The minimum exponent before (gradual) underflow, computed by
634: *          setting A = START and dividing by BASE until the previous A
635: *          can not be recovered.
636: *
637: *  START   (input) REAL
638: *          The starting point for determining EMIN.
639: *
640: *  BASE    (input) INTEGER
641: *          The base of the machine.
642: *
643: * =====================================================================
644: *
645: *     .. Local Scalars ..
646:       INTEGER            I
647:       REAL               A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
648: *     ..
649: *     .. External Functions ..
650:       REAL               SLAMC3
651:       EXTERNAL           SLAMC3
652: *     ..
653: *     .. Executable Statements ..
654: *
655:       A = START
656:       ONE = 1
657:       RBASE = ONE / BASE
658:       ZERO = 0
659:       EMIN = 1
660:       B1 = SLAMC3( A*RBASE, ZERO )
661:       C1 = A
662:       C2 = A
663:       D1 = A
664:       D2 = A
665: *+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
666: *    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP
667:    10 CONTINUE
668:       IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
669:      $    ( D2.EQ.A ) ) THEN
670:          EMIN = EMIN - 1
671:          A = B1
672:          B1 = SLAMC3( A / BASE, ZERO )
673:          C1 = SLAMC3( B1*BASE, ZERO )
674:          D1 = ZERO
675:          DO 20 I = 1, BASE
676:             D1 = D1 + B1
677:    20    CONTINUE
678:          B2 = SLAMC3( A*RBASE, ZERO )
679:          C2 = SLAMC3( B2 / RBASE, ZERO )
680:          D2 = ZERO
681:          DO 30 I = 1, BASE
682:             D2 = D2 + B2
683:    30    CONTINUE
684:          GO TO 10
685:       END IF
686: *+    END WHILE
687: *
688:       RETURN
689: *
690: *     End of SLAMC4
691: *
692:       END
693: *
694: ************************************************************************
695: *
696:       SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
697: *
698: *  -- LAPACK auxiliary routine (version 3.2) --
699: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
700: *     November 2006
701: *
702: *     .. Scalar Arguments ..
703:       LOGICAL            IEEE
704:       INTEGER            BETA, EMAX, EMIN, P
705:       REAL               RMAX
706: *     ..
707: *
708: *  Purpose
709: *  =======
710: *
711: *  SLAMC5 attempts to compute RMAX, the largest machine floating-point
712: *  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
713: *  approximately to a power of 2.  It will fail on machines where this
714: *  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
715: *  EMAX = 28718).  It will also fail if the value supplied for EMIN is
716: *  too large (i.e. too close to zero), probably with overflow.
717: *
718: *  Arguments
719: *  =========
720: *
721: *  BETA    (input) INTEGER
722: *          The base of floating-point arithmetic.
723: *
724: *  P       (input) INTEGER
725: *          The number of base BETA digits in the mantissa of a
726: *          floating-point value.
727: *
728: *  EMIN    (input) INTEGER
729: *          The minimum exponent before (gradual) underflow.
730: *
731: *  IEEE    (input) LOGICAL
732: *          A logical flag specifying whether or not the arithmetic
733: *          system is thought to comply with the IEEE standard.
734: *
735: *  EMAX    (output) INTEGER
736: *          The largest exponent before overflow
737: *
738: *  RMAX    (output) REAL
739: *          The largest machine floating-point number.
740: *
741: * =====================================================================
742: *
743: *     .. Parameters ..
744:       REAL               ZERO, ONE
745:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
746: *     ..
747: *     .. Local Scalars ..
748:       INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
749:       REAL               OLDY, RECBAS, Y, Z
750: *     ..
751: *     .. External Functions ..
752:       REAL               SLAMC3
753:       EXTERNAL           SLAMC3
754: *     ..
755: *     .. Intrinsic Functions ..
756:       INTRINSIC          MOD
757: *     ..
758: *     .. Executable Statements ..
759: *
760: *     First compute LEXP and UEXP, two powers of 2 that bound
761: *     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
762: *     approximately to the bound that is closest to abs(EMIN).
763: *     (EMAX is the exponent of the required number RMAX).
764: *
765:       LEXP = 1
766:       EXBITS = 1
767:    10 CONTINUE
768:       TRY = LEXP*2
769:       IF( TRY.LE.( -EMIN ) ) THEN
770:          LEXP = TRY
771:          EXBITS = EXBITS + 1
772:          GO TO 10
773:       END IF
774:       IF( LEXP.EQ.-EMIN ) THEN
775:          UEXP = LEXP
776:       ELSE
777:          UEXP = TRY
778:          EXBITS = EXBITS + 1
779:       END IF
780: *
781: *     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
782: *     than or equal to EMIN. EXBITS is the number of bits needed to
783: *     store the exponent.
784: *
785:       IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
786:          EXPSUM = 2*LEXP
787:       ELSE
788:          EXPSUM = 2*UEXP
789:       END IF
790: *
791: *     EXPSUM is the exponent range, approximately equal to
792: *     EMAX - EMIN + 1 .
793: *
794:       EMAX = EXPSUM + EMIN - 1
795:       NBITS = 1 + EXBITS + P
796: *
797: *     NBITS is the total number of bits needed to store a
798: *     floating-point number.
799: *
800:       IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
801: *
802: *        Either there are an odd number of bits used to store a
803: *        floating-point number, which is unlikely, or some bits are
804: *        not used in the representation of numbers, which is possible,
805: *        (e.g. Cray machines) or the mantissa has an implicit bit,
806: *        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
807: *        most likely. We have to assume the last alternative.
808: *        If this is true, then we need to reduce EMAX by one because
809: *        there must be some way of representing zero in an implicit-bit
810: *        system. On machines like Cray, we are reducing EMAX by one
811: *        unnecessarily.
812: *
813:          EMAX = EMAX - 1
814:       END IF
815: *
816:       IF( IEEE ) THEN
817: *
818: *        Assume we are on an IEEE machine which reserves one exponent
819: *        for infinity and NaN.
820: *
821:          EMAX = EMAX - 1
822:       END IF
823: *
824: *     Now create RMAX, the largest machine number, which should
825: *     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
826: *
827: *     First compute 1.0 - BETA**(-P), being careful that the
828: *     result is less than 1.0 .
829: *
830:       RECBAS = ONE / BETA
831:       Z = BETA - ONE
832:       Y = ZERO
833:       DO 20 I = 1, P
834:          Z = Z*RECBAS
835:          IF( Y.LT.ONE )
836:      $      OLDY = Y
837:          Y = SLAMC3( Y, Z )
838:    20 CONTINUE
839:       IF( Y.GE.ONE )
840:      $   Y = OLDY
841: *
842: *     Now multiply by BETA**EMAX to get RMAX.
843: *
844:       DO 30 I = 1, EMAX
845:          Y = SLAMC3( Y*BETA, ZERO )
846:    30 CONTINUE
847: *
848:       RMAX = Y
849:       RETURN
850: *
851: *     End of SLAMC5
852: *
853:       END
854: