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

◆ dlamc2()

subroutine dlamc2 ( integer  beta,
integer  t,
logical  rnd,
double precision  eps,
integer  emin,
double precision  rmin,
integer  emax,
double precision  rmax 
)

Definition at line 318 of file dlamch.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 DOUBLE PRECISION EPS, RMAX, RMIN
329* ..
330*
331* Purpose
332* =======
333*
334* DLAMC2 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) DOUBLE PRECISION
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) DOUBLE PRECISION
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) DOUBLE PRECISION
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 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
388 $ SIXTH, SMALL, THIRD, TWO, ZERO
389* ..
390* .. External Functions ..
391 DOUBLE PRECISION DLAMC3
392 EXTERNAL dlamc3
393* ..
394* .. External Subroutines ..
395 EXTERNAL dlamc1, dlamc4, dlamc5
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 DLAMC3 to ensure
419* that relevant values are stored and not held in registers, or
420* are not affected by optimizers.
421*
422* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
423*
424 CALL dlamc1( 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 = dlamc3( b, -half )
437 third = dlamc3( sixth, sixth )
438 b = dlamc3( third, -half )
439 b = dlamc3( 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 = dlamc3( half*leps, ( two**5 )*( leps**2 ) )
451 c = dlamc3( half, -c )
452 b = dlamc3( half, c )
453 c = dlamc3( half, -b )
454 b = dlamc3( 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 = dlamc3( small*rbase, zero )
472 20 CONTINUE
473 a = dlamc3( one, small )
474 CALL dlamc4( ngpmin, one, lbeta )
475 CALL dlamc4( ngnmin, -one, lbeta )
476 CALL dlamc4( gpmin, a, lbeta )
477 CALL dlamc4( 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 DLAMC1. 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 = dlamc3( lrmin*rbase, zero )
546 30 CONTINUE
547*
548* Finally, call DLAMC5 to compute EMAX and RMAX.
549*
550 CALL dlamc5( 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 $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
570*
571* End of DLAMC2
572*
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
double precision function dlamc3(a, b)
Definition tools.f:586
subroutine dlamc4(emin, start, base)
Definition tools.f:624
subroutine dlamc1(beta, t, rnd, ieee1)
Definition tools.f:140
subroutine dlamc5(beta, p, emin, ieee, emax, rmax)
Definition tools.f:708
Here is the call graph for this function: