SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slamc2()

subroutine slamc2 ( integer  beta,
integer  t,
logical  rnd,
real  eps,
integer  emin,
real  rmin,
integer  emax,
real  rmax 
)

Definition at line 318 of file slamch.f.

319*
320* -- LAPACK auxiliary routine (version 2.1) --
321* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
322* Courant Institute, Argonne National Lab, and Rice University
323* October 31, 1992
324*
325* .. Scalar Arguments ..
326 LOGICAL RND
327 INTEGER BETA, EMAX, EMIN, T
328 REAL EPS, RMAX, RMIN
329* ..
330*
331* Purpose
332* =======
333*
334* SLAMC2 determines the machine parameters specified in its argument
335* list.
336*
337* Arguments
338* =========
339*
340* BETA (output) INTEGER
341* The base of the machine.
342*
343* T (output) INTEGER
344* The number of ( BETA ) digits in the mantissa.
345*
346* RND (output) LOGICAL
347* Specifies whether proper rounding ( RND = .TRUE. ) or
348* chopping ( RND = .FALSE. ) occurs in addition. This may not
349* be a reliable guide to the way in which the machine performs
350* its arithmetic.
351*
352* EPS (output) REAL
353* The smallest positive number such that
354*
355* fl( 1.0 - EPS ) .LT. 1.0,
356*
357* where fl denotes the computed value.
358*
359* EMIN (output) INTEGER
360* The minimum exponent before (gradual) underflow occurs.
361*
362* RMIN (output) REAL
363* The smallest normalized number for the machine, given by
364* BASE**( EMIN - 1 ), where BASE is the floating point value
365* of BETA.
366*
367* EMAX (output) INTEGER
368* The maximum exponent before overflow occurs.
369*
370* RMAX (output) REAL
371* The largest positive number for the machine, given by
372* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
373* value of BETA.
374*
375* Further Details
376* ===============
377*
378* The computation of EPS is based on a routine PARANOIA by
379* W. Kahan of the University of California at Berkeley.
380*
381* =====================================================================
382*
383* .. Local Scalars ..
384 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
385 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
386 $ NGNMIN, NGPMIN
387 REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
388 $ SIXTH, SMALL, THIRD, TWO, ZERO
389* ..
390* .. External Functions ..
391 REAL SLAMC3
392 EXTERNAL slamc3
393* ..
394* .. External Subroutines ..
395 EXTERNAL slamc1, slamc4, slamc5
396* ..
397* .. Intrinsic Functions ..
398 INTRINSIC abs, max, min
399* ..
400* .. Save statement ..
401 SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
402 $ lrmin, lt
403* ..
404* .. Data statements ..
405 DATA first / .true. / , iwarn / .false. /
406* ..
407* .. Executable Statements ..
408*
409 IF( first ) THEN
410 first = .false.
411 zero = 0
412 one = 1
413 two = 2
414*
415* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
416* BETA, T, RND, EPS, EMIN and RMIN.
417*
418* Throughout this routine we use the function SLAMC3 to ensure
419* that relevant values are stored and not held in registers, or
420* are not affected by optimizers.
421*
422* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
423*
424 CALL slamc1( lbeta, lt, lrnd, lieee1 )
425*
426* Start to find EPS.
427*
428 b = lbeta
429 a = b**( -lt )
430 leps = a
431*
432* Try some tricks to see whether or not this is the correct EPS.
433*
434 b = two / 3
435 half = one / 2
436 sixth = slamc3( b, -half )
437 third = slamc3( sixth, sixth )
438 b = slamc3( third, -half )
439 b = slamc3( b, sixth )
440 b = abs( b )
441 IF( b.LT.leps )
442 $ b = leps
443*
444 leps = 1
445*
446*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
447 10 CONTINUE
448 IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
449 leps = b
450 c = slamc3( half*leps, ( two**5 )*( leps**2 ) )
451 c = slamc3( half, -c )
452 b = slamc3( half, c )
453 c = slamc3( half, -b )
454 b = slamc3( half, c )
455 GO TO 10
456 END IF
457*+ END WHILE
458*
459 IF( a.LT.leps )
460 $ leps = a
461*
462* Computation of EPS complete.
463*
464* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
465* Keep dividing A by BETA until (gradual) underflow occurs. This
466* is detected when we cannot recover the previous A.
467*
468 rbase = one / lbeta
469 small = one
470 DO 20 i = 1, 3
471 small = slamc3( small*rbase, zero )
472 20 CONTINUE
473 a = slamc3( one, small )
474 CALL slamc4( ngpmin, one, lbeta )
475 CALL slamc4( ngnmin, -one, lbeta )
476 CALL slamc4( gpmin, a, lbeta )
477 CALL slamc4( gnmin, -a, lbeta )
478 ieee = .false.
479*
480 IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
481 IF( ngpmin.EQ.gpmin ) THEN
482 lemin = ngpmin
483* ( Non twos-complement machines, no gradual underflow;
484* e.g., VAX )
485 ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
486 lemin = ngpmin - 1 + lt
487 ieee = .true.
488* ( Non twos-complement machines, with gradual underflow;
489* e.g., IEEE standard followers )
490 ELSE
491 lemin = min( ngpmin, gpmin )
492* ( A guess; no known machine )
493 iwarn = .true.
494 END IF
495*
496 ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
497 IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
498 lemin = max( ngpmin, ngnmin )
499* ( Twos-complement machines, no gradual underflow;
500* e.g., CYBER 205 )
501 ELSE
502 lemin = min( ngpmin, ngnmin )
503* ( A guess; no known machine )
504 iwarn = .true.
505 END IF
506*
507 ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
508 $ ( gpmin.EQ.gnmin ) ) THEN
509 IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
510 lemin = max( ngpmin, ngnmin ) - 1 + lt
511* ( Twos-complement machines with gradual underflow;
512* no known machine )
513 ELSE
514 lemin = min( ngpmin, ngnmin )
515* ( A guess; no known machine )
516 iwarn = .true.
517 END IF
518*
519 ELSE
520 lemin = min( ngpmin, ngnmin, gpmin, gnmin )
521* ( A guess; no known machine )
522 iwarn = .true.
523 END IF
524***
525* Comment out this if block if EMIN is ok
526 IF( iwarn ) THEN
527 first = .true.
528 WRITE( 6, fmt = 9999 )lemin
529 END IF
530***
531*
532* Assume IEEE arithmetic if we found denormalised numbers above,
533* or if arithmetic seems to round in the IEEE style, determined
534* in routine SLAMC1. A true IEEE machine should have both things
535* true; however, faulty machines may have one or the other.
536*
537 ieee = ieee .OR. lieee1
538*
539* Compute RMIN by successive division by BETA. We could compute
540* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
541* this computation.
542*
543 lrmin = 1
544 DO 30 i = 1, 1 - lemin
545 lrmin = slamc3( lrmin*rbase, zero )
546 30 CONTINUE
547*
548* Finally, call SLAMC5 to compute EMAX and RMAX.
549*
550 CALL slamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
551 END IF
552*
553 beta = lbeta
554 t = lt
555 rnd = lrnd
556 eps = leps
557 emin = lemin
558 rmin = lrmin
559 emax = lemax
560 rmax = lrmax
561*
562 RETURN
563*
564 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
565 $ ' EMIN = ', i8, /
566 $ ' If, after inspection, the value EMIN looks',
567 $ ' acceptable please comment out ',
568 $ / ' the IF block as marked within the code of routine',
569 $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
570*
571* End of SLAMC2
572*
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
subroutine slamc5(beta, p, emin, ieee, emax, rmax)
Definition tools.f:1565
subroutine slamc1(beta, t, rnd, ieee1)
Definition tools.f:997
subroutine slamc4(emin, start, base)
Definition tools.f:1481
real function slamc3(a, b)
Definition tools.f:1443
Here is the call graph for this function: