OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
tools.f File Reference

Go to the source code of this file.

Functions/Subroutines

double precision function dlamch (cmach)
subroutine dlamc1 (beta, t, rnd, ieee1)
subroutine dlamc2 (beta, t, rnd, eps, emin, rmin, emax, rmax)
double precision function dlamc3 (a, b)
subroutine dlamc4 (emin, start, base)
subroutine dlamc5 (beta, p, emin, ieee, emax, rmax)
real function slamch (cmach)
subroutine slamc1 (beta, t, rnd, ieee1)
subroutine slamc2 (beta, t, rnd, eps, emin, rmin, emax, rmax)
real function slamc3 (a, b)
subroutine slamc4 (emin, start, base)
subroutine slamc5 (beta, p, emin, ieee, emax, rmax)
logical function lsame (ca, cb)
double precision function dlarnd (idist, iseed)
double complex function zlarnd (idist, iseed)
double precision function dlaran (iseed)

Function/Subroutine Documentation

◆ dlamc1()

subroutine dlamc1 ( integer beta,
integer t,
logical rnd,
logical ieee1 )

Definition at line 139 of file tools.f.

140*
141* -- LAPACK auxiliary routine (version 2.0) --
142* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
143* Courant Institute, Argonne National Lab, and Rice University
144* October 31, 1992
145*
146* .. Scalar Arguments ..
147 LOGICAL IEEE1, RND
148 INTEGER BETA, T
149* ..
150*
151* Purpose
152* =======
153*
154* DLAMC1 determines the machine parameters given by BETA, T, RND, and
155* IEEE1.
156*
157* Arguments
158* =========
159*
160* BETA (output) INTEGER
161* The base of the machine.
162*
163* T (output) INTEGER
164* The number of ( BETA ) digits in the mantissa.
165*
166* RND (output) LOGICAL
167* Specifies whether proper rounding ( RND = .TRUE. ) or
168* chopping ( RND = .FALSE. ) occurs in addition. This may not
169* be a reliable guide to the way in which the machine performs
170* its arithmetic.
171*
172* IEEE1 (output) LOGICAL
173* Specifies whether rounding appears to be done in the IEEE
174* 'round to nearest' style.
175*
176* Further Details
177* ===============
178*
179* The routine is based on the routine ENVRON by Malcolm and
180* incorporates suggestions by Gentleman and Marovich. See
181*
182* Malcolm M. A. (1972) Algorithms to reveal properties of
183* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
184*
185* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
186* that reveal properties of floating point arithmetic units.
187* Comms. of the ACM, 17, 276-277.
188*
189* =====================================================================
190*
191* .. Local Scalars ..
192 LOGICAL FIRST, LIEEE1, LRND
193 INTEGER LBETA, LT
194 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
195* ..
196* .. External Functions ..
197 DOUBLE PRECISION DLAMC3
198 EXTERNAL dlamc3
199* ..
200* .. Save statement ..
201 SAVE first, lieee1, lbeta, lrnd, lt
202* ..
203* .. Data statements ..
204 DATA first / .true. /
205* ..
206* .. Executable Statements ..
207*
208 IF( first ) THEN
209 first = .false.
210 one = 1
211*
212* LBETA, LIEEE1, LT and LRND are the local values of BETA,
213* IEEE1, T and RND.
214*
215* Throughout this routine we use the function DLAMC3 to ensure
216* that relevant values are stored and not held in registers, or
217* are not affected by optimizers.
218*
219* Compute a = 2.0**m with the smallest positive integer m such
220* that
221*
222* fl( a + 1.0 ) = a.
223*
224 a = 1
225 c = 1
226*
227*+ WHILE( C.EQ.ONE )LOOP
228 10 CONTINUE
229 IF( c.EQ.one ) THEN
230 a = 2*a
231 c = dlamc3( a, one )
232 c = dlamc3( c, -a )
233 GO TO 10
234 END IF
235*+ END WHILE
236*
237* Now compute b = 2.0**m with the smallest positive integer m
238* such that
239*
240* fl( a + b ) .gt. a.
241*
242 b = 1
243 c = dlamc3( a, b )
244*
245*+ WHILE( C.EQ.A )LOOP
246 20 CONTINUE
247 IF( c.EQ.a ) THEN
248 b = 2*b
249 c = dlamc3( a, b )
250 GO TO 20
251 END IF
252*+ END WHILE
253*
254* Now compute the base. a and c are neighbouring floating point
255* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
256* their difference is beta. Adding 0.25 to c is to ensure that it
257* is truncated to beta and not ( beta - 1 ).
258*
259 qtr = one / 4
260 savec = c
261 c = dlamc3( c, -a )
262 lbeta = c + qtr
263*
264* Now determine whether rounding or chopping occurs, by adding a
265* bit less than beta/2 and a bit more than beta/2 to a.
266*
267 b = lbeta
268 f = dlamc3( b / 2, -b / 100 )
269 c = dlamc3( f, a )
270 IF( c.EQ.a ) THEN
271 lrnd = .true.
272 ELSE
273 lrnd = .false.
274 END IF
275 f = dlamc3( b / 2, b / 100 )
276 c = dlamc3( f, a )
277 IF( ( lrnd ) .AND. ( c.EQ.a ) )
278 $ lrnd = .false.
279*
280* Try and decide whether rounding is done in the IEEE 'round to
281* nearest' style. B/2 is half a unit in the last place of the two
282* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
283* zero, and SAVEC is odd. Thus adding B/2 to A should not change
284* A, but adding B/2 to SAVEC should change SAVEC.
285*
286 t1 = dlamc3( b / 2, a )
287 t2 = dlamc3( b / 2, savec )
288 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
289*
290* Now find the mantissa, t. It should be the integer part of
291* log to the base beta of a, however it is safer to determine t
292* by powering. So we find t as the smallest positive integer for
293* which
294*
295* fl( beta**t + 1.0 ) = 1.0.
296*
297 lt = 0
298 a = 1
299 c = 1
300*
301*+ WHILE( C.EQ.ONE )LOOP
302 30 CONTINUE
303 IF( c.EQ.one ) THEN
304 lt = lt + 1
305 a = a*lbeta
306 c = dlamc3( a, one )
307 c = dlamc3( c, -a )
308 GO TO 30
309 END IF
310*+ END WHILE
311*
312 END IF
313*
314 beta = lbeta
315 t = lt
316 rnd = lrnd
317 ieee1 = lieee1
318 RETURN
319*
320* End of DLAMC1
321*
dlamc3
double precision function dlamc3(a, b)
Definition tools.f:586

◆ 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 326 of file tools.f.

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

◆ dlamc3()

double precision function dlamc3 ( double precision a,
double precision b )

Definition at line 585 of file tools.f.

586*
587* -- LAPACK auxiliary routine (version 2.0) --
588* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
589* Courant Institute, Argonne National Lab, and Rice University
590* October 31, 1992
591*
592* .. Scalar Arguments ..
593 DOUBLE PRECISION A, B
594* ..
595*
596* Purpose
597* =======
598*
599* DLAMC3 is intended to force A and B to be stored prior to doing
600* the addition of A and B , for use in situations where optimizers
601* might hold one of these in a register.
602*
603* Arguments
604* =========
605*
606* A, B (input) DOUBLE PRECISION
607* The values A and B.
608*
609* =====================================================================
610*
611* .. Executable Statements ..
612*
613 dlamc3 = a + b
614*
615 RETURN
616*
617* End of DLAMC3
618*

◆ dlamc4()

subroutine dlamc4 ( integer emin,
double precision start,
integer base )

Definition at line 623 of file tools.f.

624*
625* -- LAPACK auxiliary routine (version 2.0) --
626* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
627* Courant Institute, Argonne National Lab, and Rice University
628* October 31, 1992
629*
630* .. Scalar Arguments ..
631 INTEGER BASE, EMIN
632 DOUBLE PRECISION START
633* ..
634*
635* Purpose
636* =======
637*
638* DLAMC4 is a service routine for DLAMC2.
639*
640* Arguments
641* =========
642*
643* EMIN (output) EMIN
644* The minimum exponent before (gradual) underflow, computed by
645* setting A = START and dividing by BASE until the previous A
646* can not be recovered.
647*
648* START (input) DOUBLE PRECISION
649* The starting point for determining EMIN.
650*
651* BASE (input) INTEGER
652* The base of the machine.
653*
654* =====================================================================
655*
656* .. Local Scalars ..
657 INTEGER I
658 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
659* ..
660* .. External Functions ..
661 DOUBLE PRECISION DLAMC3
662 EXTERNAL dlamc3
663* ..
664* .. Executable Statements ..
665*
666 a = start
667 one = 1
668 rbase = one / base
669 zero = 0
670 emin = 1
671 b1 = dlamc3( a*rbase, zero )
672 c1 = a
673 c2 = a
674 d1 = a
675 d2 = a
676*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
677* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
678 10 CONTINUE
679 IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
680 $ ( d2.EQ.a ) ) THEN
681 emin = emin - 1
682 a = b1
683 b1 = dlamc3( a / base, zero )
684 c1 = dlamc3( b1*base, zero )
685 d1 = zero
686 DO 20 i = 1, base
687 d1 = d1 + b1
688 20 CONTINUE
689 b2 = dlamc3( a*rbase, zero )
690 c2 = dlamc3( b2 / rbase, zero )
691 d2 = zero
692 DO 30 i = 1, base
693 d2 = d2 + b2
694 30 CONTINUE
695 GO TO 10
696 END IF
697*+ END WHILE
698*
699 RETURN
700*
701* End of DLAMC4
702*

◆ dlamc5()

subroutine dlamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
double precision rmax )

Definition at line 707 of file tools.f.

708*
709* -- LAPACK auxiliary routine (version 2.0) --
710* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
711* Courant Institute, Argonne National Lab, and Rice University
712* October 31, 1992
713*
714* .. Scalar Arguments ..
715 LOGICAL IEEE
716 INTEGER BETA, EMAX, EMIN, P
717 DOUBLE PRECISION RMAX
718* ..
719*
720* Purpose
721* =======
722*
723* DLAMC5 attempts to compute RMAX, the largest machine floating-point
724* number, without overflow. It assumes that EMAX + abs(EMIN) sum
725* approximately to a power of 2. It will fail on machines where this
726* assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
727* EMAX = 28718). It will also fail if the value supplied for EMIN is
728* too large (i.e. too close to zero), probably with overflow.
729*
730* Arguments
731* =========
732*
733* BETA (input) INTEGER
734* The base of floating-point arithmetic.
735*
736* P (input) INTEGER
737* The number of base BETA digits in the mantissa of a
738* floating-point value.
739*
740* EMIN (input) INTEGER
741* The minimum exponent before (gradual) underflow.
742*
743* IEEE (input) LOGICAL
744* A logical flag specifying whether or not the arithmetic
745* system is thought to comply with the IEEE standard.
746*
747* EMAX (output) INTEGER
748* The largest exponent before overflow
749*
750* RMAX (output) DOUBLE PRECISION
751* The largest machine floating-point number.
752*
753* =====================================================================
754*
755* .. Parameters ..
756 DOUBLE PRECISION ZERO, ONE
757 parameter( zero = 0.0d0, one = 1.0d0 )
758* ..
759* .. Local Scalars ..
760 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
761 DOUBLE PRECISION OLDY, RECBAS, Y, Z
762* ..
763* .. External Functions ..
764 DOUBLE PRECISION DLAMC3
765 EXTERNAL dlamc3
766* ..
767* .. Intrinsic Functions ..
768 INTRINSIC mod
769* ..
770* .. Executable Statements ..
771*
772* First compute LEXP and UEXP, two powers of 2 that bound
773* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
774* approximately to the bound that is closest to abs(EMIN).
775* (EMAX is the exponent of the required number RMAX).
776*
777 lexp = 1
778 exbits = 1
779 10 CONTINUE
780 try = lexp*2
781 IF( try.LE.( -emin ) ) THEN
782 lexp = try
783 exbits = exbits + 1
784 GO TO 10
785 END IF
786 IF( lexp.EQ.-emin ) THEN
787 uexp = lexp
788 ELSE
789 uexp = try
790 exbits = exbits + 1
791 END IF
792*
793* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
794* than or equal to EMIN. EXBITS is the number of bits needed to
795* store the exponent.
796*
797 IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
798 expsum = 2*lexp
799 ELSE
800 expsum = 2*uexp
801 END IF
802*
803* EXPSUM is the exponent range, approximately equal to
804* EMAX - EMIN + 1 .
805*
806 emax = expsum + emin - 1
807 nbits = 1 + exbits + p
808*
809* NBITS is the total number of bits needed to store a
810* floating-point number.
811*
812 IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
813*
814* Either there are an odd number of bits used to store a
815* floating-point number, which is unlikely, or some bits are
816* not used in the representation of numbers, which is possible,
817* (e.g. Cray machines) or the mantissa has an implicit bit,
818* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
819* most likely. We have to assume the last alternative.
820* If this is true, then we need to reduce EMAX by one because
821* there must be some way of representing zero in an implicit-bit
822* system. On machines like Cray, we are reducing EMAX by one
823* unnecessarily.
824*
825 emax = emax - 1
826 END IF
827*
828 IF( ieee ) THEN
829*
830* Assume we are on an IEEE machine which reserves one exponent
831* for infinity and NaN.
832*
833 emax = emax - 1
834 END IF
835*
836* Now create RMAX, the largest machine number, which should
837* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
838*
839* First compute 1.0 - BETA**(-P), being careful that the
840* result is less than 1.0 .
841*
842 recbas = one / beta
843 z = beta - one
844 y = zero
845 DO 20 i = 1, p
846 z = z*recbas
847 IF( y.LT.one )
848 $ oldy = y
849 y = dlamc3( y, z )
850 20 CONTINUE
851 IF( y.GE.one )
852 $ y = oldy
853*
854* Now multiply by BETA**EMAX to get RMAX.
855*
856 DO 30 i = 1, emax
857 y = dlamc3( y*beta, zero )
858 30 CONTINUE
859*
860 rmax = y
861 RETURN
862*
863* End of DLAMC5
864*

◆ dlamch()

double precision function dlamch ( character cmach)

Definition at line 9 of file tools.f.

10*
11* -- LAPACK auxiliary routine (version 2.0) --
12* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
13* Courant Institute, Argonne National Lab, and Rice University
14* October 31, 1992
15*
16* .. Scalar Arguments ..
17 CHARACTER CMACH
18* ..
19*
20* Purpose
21* =======
22*
23* DLAMCH determines double precision machine parameters.
24*
25* Arguments
26* =========
27*
28* CMACH (input) CHARACTER*1
29* Specifies the value to be returned by DLAMCH:
30* = 'E' or 'e', DLAMCH := eps
31* = 'S' or 's , DLAMCH := sfmin
32* = 'B' or 'b', DLAMCH := base
33* = 'P' or 'p', DLAMCH := eps*base
34* = 'N' or 'n', DLAMCH := t
35* = 'R' or 'r', DLAMCH := rnd
36* = 'M' or 'm', DLAMCH := emin
37* = 'U' or 'u', DLAMCH := rmin
38* = 'L' or 'l', DLAMCH := emax
39* = 'O' or 'o', DLAMCH := rmax
40*
41* where
42*
43* eps = relative machine precision
44* sfmin = safe minimum, such that 1/sfmin does not overflow
45* base = base of the machine
46* prec = eps*base
47* t = number of (base) digits in the mantissa
48* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
49* emin = minimum exponent before (gradual) underflow
50* rmin = underflow threshold - base**(emin-1)
51* emax = largest exponent before overflow
52* rmax = overflow threshold - (base**emax)*(1-eps)
53*
54* =====================================================================
55*
56* .. Parameters ..
57 DOUBLE PRECISION ONE, ZERO
58 parameter( one = 1.0d+0, zero = 0.0d+0 )
59* ..
60* .. Local Scalars ..
61 LOGICAL FIRST, LRND
62 INTEGER BETA, IMAX, IMIN, IT
63 DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
64 $ RND, SFMIN, SMALL, T
65* ..
66* .. External Functions ..
67 LOGICAL LSAME
68 EXTERNAL lsame
69* ..
70* .. External Subroutines ..
71 EXTERNAL dlamc2
72* ..
73* .. Save statement ..
74 SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
75 $ emax, rmax, prec
76* ..
77* .. Data statements ..
78 DATA first / .true. /
79* ..
80* .. Executable Statements ..
81*
82 IF( first ) THEN
83 first = .false.
84 CALL dlamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
85 base = beta
86 t = it
87 IF( lrnd ) THEN
88 rnd = one
89 eps = ( base**( 1-it ) ) / 2
90 ELSE
91 rnd = zero
92 eps = base**( 1-it )
93 END IF
94 prec = eps*base
95 emin = imin
96 emax = imax
97 sfmin = rmin
98 small = one / rmax
99 IF( small.GE.sfmin ) THEN
100*
101* Use SMALL plus a bit, to avoid the possibility of rounding
102* causing overflow when computing 1/sfmin.
103*
104 sfmin = small*( one+eps )
105 END IF
106 END IF
107*
108 IF( lsame( cmach, 'E' ) ) THEN
109 rmach = eps
110 ELSE IF( lsame( cmach, 'S' ) ) THEN
111 rmach = sfmin
112 ELSE IF( lsame( cmach, 'B' ) ) THEN
113 rmach = base
114 ELSE IF( lsame( cmach, 'P' ) ) THEN
115 rmach = prec
116 ELSE IF( lsame( cmach, 'N' ) ) THEN
117 rmach = t
118 ELSE IF( lsame( cmach, 'R' ) ) THEN
119 rmach = rnd
120 ELSE IF( lsame( cmach, 'M' ) ) THEN
121 rmach = emin
122 ELSE IF( lsame( cmach, 'U' ) ) THEN
123 rmach = rmin
124 ELSE IF( lsame( cmach, 'L' ) ) THEN
125 rmach = emax
126 ELSE IF( lsame( cmach, 'O' ) ) THEN
127 rmach = rmax
128 END IF
129*
130 dlamch = rmach
131 RETURN
132*
133* End of DLAMCH
134*
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
dlamch
double precision function dlamch(cmach)
Definition tools.f:10
dlamc2
subroutine dlamc2(beta, t, rnd, eps, emin, rmin, emax, rmax)
Definition tools.f:327

◆ dlaran()

double precision function dlaran ( integer, dimension( 4 ) iseed)

Definition at line 1999 of file tools.f.

2000*
2001* -- LAPACK auxiliary routine (version 2.0) --
2002* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
2003* Courant Institute, Argonne National Lab, and Rice University
2004* February 29, 1992
2005*
2006* .. Array Arguments ..
2007 INTEGER ISEED( 4 )
2008* ..
2009*
2010* Purpose
2011* =======
2012*
2013* DLARAN returns a random real number from a uniform (0,1)
2014* distribution.
2015*
2016* Arguments
2017* =========
2018*
2019* ISEED (input/output) INTEGER array, dimension (4)
2020* On entry, the seed of the random number generator; the array
2021* elements must be between 0 and 4095, and ISEED(4) must be
2022* odd.
2023* On exit, the seed is updated.
2024*
2025* Further Details
2026* ===============
2027*
2028* This routine uses a multiplicative congruential method with modulus
2029* 2**48 and multiplier 33952834046453 (see G.S.Fishman,
2030* 'Multiplicative congruential random number generators with modulus
2031* 2**b: an exhaustive analysis for b = 32 and a partial analysis for
2032* b = 48', Math. Comp. 189, pp 331-344, 1990).
2033*
2034* 48-bit integers are stored in 4 integer array elements with 12 bits
2035* per element. Hence the routine is portable across machines with
2036* integers of 32 bits or more.
2037*
2038* =====================================================================
2039*
2040* .. Parameters ..
2041 INTEGER M1, M2, M3, M4
2042 parameter( m1 = 494, m2 = 322, m3 = 2508, m4 = 2549 )
2043 DOUBLE PRECISION ONE
2044 parameter( one = 1.0d+0 )
2045 INTEGER IPW2
2046 DOUBLE PRECISION R
2047 parameter( ipw2 = 4096, r = one / ipw2 )
2048* ..
2049* .. Local Scalars ..
2050 INTEGER IT1, IT2, IT3, IT4
2051* ..
2052* .. Intrinsic Functions ..
2053 INTRINSIC dble, mod
2054* ..
2055* .. Executable Statements ..
2056*
2057* multiply the seed by the multiplier modulo 2**48
2058*
2059 it4 = iseed( 4 )*m4
2060 it3 = it4 / ipw2
2061 it4 = it4 - ipw2*it3
2062 it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
2063 it2 = it3 / ipw2
2064 it3 = it3 - ipw2*it2
2065 it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
2066 it1 = it2 / ipw2
2067 it2 = it2 - ipw2*it1
2068 it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 +
2069 $ iseed( 4 )*m1
2070 it1 = mod( it1, ipw2 )
2071*
2072* return updated seed
2073*
2074 iseed( 1 ) = it1
2075 iseed( 2 ) = it2
2076 iseed( 3 ) = it3
2077 iseed( 4 ) = it4
2078*
2079* convert 48-bit integer to a real number in the interval (0,1)
2080*
2081 dlaran = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
2082 $ ( dble( it4 ) ) ) ) )
2083 RETURN
2084*
2085* End of DLARAN
2086*
dlaran
double precision function dlaran(iseed)
Definition tools.f:2000

◆ dlarnd()

double precision function dlarnd ( integer idist,
integer, dimension( 4 ) iseed )

Definition at line 1810 of file tools.f.

1811*
1812* -- LAPACK auxiliary routine (version 2.0) --
1813* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1814* Courant Institute, Argonne National Lab, and Rice University
1815* June 30, 1994
1816*
1817* .. Scalar Arguments ..
1818 INTEGER IDIST
1819* ..
1820* .. Array Arguments ..
1821 INTEGER ISEED( 4 )
1822* ..
1823*
1824* Purpose
1825* =======
1826*
1827* DLARND returns a random real number from a uniform or normal
1828* distribution.
1829*
1830* Arguments
1831* =========
1832*
1833* IDIST (input) INTEGER
1834* Specifies the distribution of the random numbers:
1835* = 1: uniform (0,1)
1836* = 2: uniform (-1,1)
1837* = 3: normal (0,1)
1838*
1839* ISEED (input/output) INTEGER array, dimension (4)
1840* On entry, the seed of the random number generator; the array
1841* elements must be between 0 and 4095, and ISEED(4) must be
1842* odd.
1843* On exit, the seed is updated.
1844*
1845* Further Details
1846* ===============
1847*
1848* This routine calls the auxiliary routine DLARAN to generate a random
1849* real number from a uniform (0,1) distribution. The Box-Muller method
1850* is used to transform numbers from a uniform to a normal distribution.
1851*
1852* =====================================================================
1853*
1854* .. Parameters ..
1855 DOUBLE PRECISION ONE, TWO
1856 parameter( one = 1.0d+0, two = 2.0d+0 )
1857 DOUBLE PRECISION TWOPI
1858 parameter( twopi = 6.2831853071795864769252867663d+0 )
1859* ..
1860* .. Local Scalars ..
1861 DOUBLE PRECISION T1, T2
1862* ..
1863* .. External Functions ..
1864 DOUBLE PRECISION DLARAN
1865 EXTERNAL dlaran
1866* ..
1867* .. Intrinsic Functions ..
1868 INTRINSIC cos, log, sqrt
1869* ..
1870* .. Executable Statements ..
1871*
1872* Generate a real random number from a uniform (0,1) distribution
1873*
1874 t1 = dlaran( iseed )
1875*
1876 IF( idist.EQ.1 ) THEN
1877*
1878* uniform (0,1)
1879*
1880 dlarnd = t1
1881 ELSE IF( idist.EQ.2 ) THEN
1882*
1883* uniform (-1,1)
1884*
1885 dlarnd = two*t1 - one
1886 ELSE IF( idist.EQ.3 ) THEN
1887*
1888* normal (0,1)
1889*
1890 t2 = dlaran( iseed )
1891 dlarnd = sqrt( -two*log( t1 ) )*cos( twopi*t2 )
1892 END IF
1893 RETURN
1894*
1895* End of DLARND
1896*
dlarnd
double precision function dlarnd(idist, iseed)
Definition tools.f:1811

◆ lsame()

logical function lsame ( character ca,
character cb )

Definition at line 1723 of file tools.f.

1724*
1725* -- LAPACK auxiliary routine (version 2.0) --
1726* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1727* Courant Institute, Argonne National Lab, and Rice University
1728* June 30, 1994
1729*
1730* .. Scalar Arguments ..
1731 CHARACTER CA, CB
1732* ..
1733*
1734* Purpose
1735* =======
1736*
1737* LSAME returns .TRUE. if CA is the same letter as CB regardless of
1738* case.
1739*
1740* Arguments
1741* =========
1742*
1743* CA (input) CHARACTER*1
1744* CB (input) CHARACTER*1
1745* CA and CB specify the single characters to be compared.
1746*
1747* =====================================================================
1748*
1749* .. Intrinsic Functions ..
1750 INTRINSIC ichar
1751* ..
1752* .. Local Scalars ..
1753 INTEGER INTA, INTB, ZCODE
1754* ..
1755* .. Executable Statements ..
1756*
1757* Test if the characters are equal
1758*
1759 lsame = ca.EQ.cb
1760 IF( lsame )
1761 $ RETURN
1762*
1763* Now test for equivalence if both characters are alphabetic.
1764*
1765 zcode = ichar( 'Z' )
1766*
1767* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
1768* machines, on which ICHAR returns a value with bit 8 set.
1769* ICHAR('A') on Prime machines returns 193 which is the same as
1770* ICHAR('A') on an EBCDIC machine.
1771*
1772 inta = ichar( ca )
1773 intb = ichar( cb )
1774*
1775 IF( zcode.EQ.90 .OR. zcode.EQ.122 ) THEN
1776*
1777* ASCII is assumed - ZCODE is the ASCII code of either lower or
1778* upper case 'Z'.
1779*
1780 IF( inta.GE.97 .AND. inta.LE.122 ) inta = inta - 32
1781 IF( intb.GE.97 .AND. intb.LE.122 ) intb = intb - 32
1782*
1783 ELSE IF( zcode.EQ.233 .OR. zcode.EQ.169 ) THEN
1784*
1785* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
1786* upper case 'Z'.
1787*
1788 IF( inta.GE.129 .AND. inta.LE.137 .OR.
1789 $ inta.GE.145 .AND. inta.LE.153 .OR.
1790 $ inta.GE.162 .AND. inta.LE.169 ) inta = inta + 64
1791 IF( intb.GE.129 .AND. intb.LE.137 .OR.
1792 $ intb.GE.145 .AND. intb.LE.153 .OR.
1793 $ intb.GE.162 .AND. intb.LE.169 ) intb = intb + 64
1794*
1795 ELSE IF( zcode.EQ.218 .OR. zcode.EQ.250 ) THEN
1796*
1797* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
1798* plus 128 of either lower or upper case 'Z'.
1799*
1800 IF( inta.GE.225 .AND. inta.LE.250 ) inta = inta - 32
1801 IF( intb.GE.225 .AND. intb.LE.250 ) intb = intb - 32
1802 END IF
1803 lsame = inta.EQ.intb
1804*
1805* RETURN
1806*
1807* End of LSAME
1808*

◆ slamc1()

subroutine slamc1 ( integer beta,
integer t,
logical rnd,
logical ieee1 )

Definition at line 996 of file tools.f.

997*
998* -- LAPACK auxiliary routine (version 2.0) --
999* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1000* Courant Institute, Argonne National Lab, and Rice University
1001* October 31, 1992
1002*
1003* .. Scalar Arguments ..
1004 LOGICAL IEEE1, RND
1005 INTEGER BETA, T
1006* ..
1007*
1008* Purpose
1009* =======
1010*
1011* SLAMC1 determines the machine parameters given by BETA, T, RND, and
1012* IEEE1.
1013*
1014* Arguments
1015* =========
1016*
1017* BETA (output) INTEGER
1018* The base of the machine.
1019*
1020* T (output) INTEGER
1021* The number of ( BETA ) digits in the mantissa.
1022*
1023* RND (output) LOGICAL
1024* Specifies whether proper rounding ( RND = .TRUE. ) or
1025* chopping ( RND = .FALSE. ) occurs in addition. This may not
1026* be a reliable guide to the way in which the machine performs
1027* its arithmetic.
1028*
1029* IEEE1 (output) LOGICAL
1030* Specifies whether rounding appears to be done in the IEEE
1031* 'round to nearest' style.
1032*
1033* Further Details
1034* ===============
1035*
1036* The routine is based on the routine ENVRON by Malcolm and
1037* incorporates suggestions by Gentleman and Marovich. See
1038*
1039* Malcolm M. A. (1972) Algorithms to reveal properties of
1040* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
1041*
1042* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
1043* that reveal properties of floating point arithmetic units.
1044* Comms. of the ACM, 17, 276-277.
1045*
1046* =====================================================================
1047*
1048* .. Local Scalars ..
1049 LOGICAL FIRST, LIEEE1, LRND
1050 INTEGER LBETA, LT
1051 REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
1052* ..
1053* .. External Functions ..
1054 REAL SLAMC3
1055 EXTERNAL slamc3
1056* ..
1057* .. Save statement ..
1058 SAVE first, lieee1, lbeta, lrnd, lt
1059* ..
1060* .. Data statements ..
1061 DATA first / .true. /
1062* ..
1063* .. Executable Statements ..
1064*
1065 IF( first ) THEN
1066 first = .false.
1067 one = 1
1068*
1069* LBETA, LIEEE1, LT and LRND are the local values of BETA,
1070* IEEE1, T and RND.
1071*
1072* Throughout this routine we use the function SLAMC3 to ensure
1073* that relevant values are stored and not held in registers, or
1074* are not affected by optimizers.
1075*
1076* Compute a = 2.0**m with the smallest positive integer m such
1077* that
1078*
1079* fl( a + 1.0 ) = a.
1080*
1081 a = 1
1082 c = 1
1083*
1084*+ WHILE( C.EQ.ONE )LOOP
1085 10 CONTINUE
1086 IF( c.EQ.one ) THEN
1087 a = 2*a
1088 c = slamc3( a, one )
1089 c = slamc3( c, -a )
1090 GO TO 10
1091 END IF
1092*+ END WHILE
1093*
1094* Now compute b = 2.0**m with the smallest positive integer m
1095* such that
1096*
1097* fl( a + b ) .gt. a.
1098*
1099 b = 1
1100 c = slamc3( a, b )
1101*
1102*+ WHILE( C.EQ.A )LOOP
1103 20 CONTINUE
1104 IF( c.EQ.a ) THEN
1105 b = 2*b
1106 c = slamc3( a, b )
1107 GO TO 20
1108 END IF
1109*+ END WHILE
1110*
1111* Now compute the base. a and c are neighbouring floating point
1112* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
1113* their difference is beta. Adding 0.25 to c is to ensure that it
1114* is truncated to beta and not ( beta - 1 ).
1115*
1116 qtr = one / 4
1117 savec = c
1118 c = slamc3( c, -a )
1119 lbeta = c + qtr
1120*
1121* Now determine whether rounding or chopping occurs, by adding a
1122* bit less than beta/2 and a bit more than beta/2 to a.
1123*
1124 b = lbeta
1125 f = slamc3( b / 2, -b / 100 )
1126 c = slamc3( f, a )
1127 IF( c.EQ.a ) THEN
1128 lrnd = .true.
1129 ELSE
1130 lrnd = .false.
1131 END IF
1132 f = slamc3( b / 2, b / 100 )
1133 c = slamc3( f, a )
1134 IF( ( lrnd ) .AND. ( c.EQ.a ) )
1135 $ lrnd = .false.
1136*
1137* Try and decide whether rounding is done in the IEEE 'round to
1138* nearest' style. B/2 is half a unit in the last place of the two
1139* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
1140* zero, and SAVEC is odd. Thus adding B/2 to A should not change
1141* A, but adding B/2 to SAVEC should change SAVEC.
1142*
1143 t1 = slamc3( b / 2, a )
1144 t2 = slamc3( b / 2, savec )
1145 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
1146*
1147* Now find the mantissa, t. It should be the integer part of
1148* log to the base beta of a, however it is safer to determine t
1149* by powering. So we find t as the smallest positive integer for
1150* which
1151*
1152* fl( beta**t + 1.0 ) = 1.0.
1153*
1154 lt = 0
1155 a = 1
1156 c = 1
1157*
1158*+ WHILE( C.EQ.ONE )LOOP
1159 30 CONTINUE
1160 IF( c.EQ.one ) THEN
1161 lt = lt + 1
1162 a = a*lbeta
1163 c = slamc3( a, one )
1164 c = slamc3( c, -a )
1165 GO TO 30
1166 END IF
1167*+ END WHILE
1168*
1169 END IF
1170*
1171 beta = lbeta
1172 t = lt
1173 rnd = lrnd
1174 ieee1 = lieee1
1175 RETURN
1176*
1177* End of SLAMC1
1178*
slamc3
real function slamc3(a, b)
Definition tools.f:1443

◆ slamc2()

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

Definition at line 1183 of file tools.f.

1184*
1185* -- LAPACK auxiliary routine (version 2.0) --
1186* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1187* Courant Institute, Argonne National Lab, and Rice University
1188* October 31, 1992
1189*
1190* .. Scalar Arguments ..
1191 LOGICAL RND
1192 INTEGER BETA, EMAX, EMIN, T
1193 REAL EPS, RMAX, RMIN
1194* ..
1195*
1196* Purpose
1197* =======
1198*
1199* SLAMC2 determines the machine parameters specified in its argument
1200* list.
1201*
1202* Arguments
1203* =========
1204*
1205* BETA (output) INTEGER
1206* The base of the machine.
1207*
1208* T (output) INTEGER
1209* The number of ( BETA ) digits in the mantissa.
1210*
1211* RND (output) LOGICAL
1212* Specifies whether proper rounding ( RND = .TRUE. ) or
1213* chopping ( RND = .FALSE. ) occurs in addition. This may not
1214* be a reliable guide to the way in which the machine performs
1215* its arithmetic.
1216*
1217* EPS (output) REAL
1218* The smallest positive number such that
1219*
1220* fl( 1.0 - EPS ) .LT. 1.0,
1221*
1222* where fl denotes the computed value.
1223*
1224* EMIN (output) INTEGER
1225* The minimum exponent before (gradual) underflow occurs.
1226*
1227* RMIN (output) REAL
1228* The smallest normalized number for the machine, given by
1229* BASE**( EMIN - 1 ), where BASE is the floating point value
1230* of BETA.
1231*
1232* EMAX (output) INTEGER
1233* The maximum exponent before overflow occurs.
1234*
1235* RMAX (output) REAL
1236* The largest positive number for the machine, given by
1237* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
1238* value of BETA.
1239*
1240* Further Details
1241* ===============
1242*
1243* The computation of EPS is based on a routine PARANOIA by
1244* W. Kahan of the University of California at Berkeley.
1245*
1246* =====================================================================
1247*
1248* .. Local Scalars ..
1249 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
1250 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
1251 $ NGNMIN, NGPMIN
1252 REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
1253 $ SIXTH, SMALL, THIRD, TWO, ZERO
1254* ..
1255* .. External Functions ..
1256 REAL SLAMC3
1257 EXTERNAL slamc3
1258* ..
1259* .. External Subroutines ..
1260 EXTERNAL slamc1, slamc4, slamc5
1261* ..
1262* .. Intrinsic Functions ..
1263 INTRINSIC abs, max, min
1264* ..
1265* .. Save statement ..
1266 SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
1267 $ lrmin, lt
1268* ..
1269* .. Data statements ..
1270 DATA first / .true. / , iwarn / .false. /
1271* ..
1272* .. Executable Statements ..
1273*
1274 IF( first ) THEN
1275 first = .false.
1276 zero = 0
1277 one = 1
1278 two = 2
1279*
1280* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
1281* BETA, T, RND, EPS, EMIN and RMIN.
1282*
1283* Throughout this routine we use the function SLAMC3 to ensure
1284* that relevant values are stored and not held in registers, or
1285* are not affected by optimizers.
1286*
1287* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
1288*
1289 CALL slamc1( lbeta, lt, lrnd, lieee1 )
1290*
1291* Start to find EPS.
1292*
1293 b = lbeta
1294 a = b**( -lt )
1295 leps = a
1296*
1297* Try some tricks to see whether or not this is the correct EPS.
1298*
1299 b = two / 3
1300 half = one / 2
1301 sixth = slamc3( b, -half )
1302 third = slamc3( sixth, sixth )
1303 b = slamc3( third, -half )
1304 b = slamc3( b, sixth )
1305 b = abs( b )
1306 IF( b.LT.leps )
1307 $ b = leps
1308*
1309 leps = 1
1310*
1311*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
1312 10 CONTINUE
1313 IF( ( leps.GT.b ) .AND. ( b.GT.zero ) ) THEN
1314 leps = b
1315 c = slamc3( half*leps, ( two**5 )*( leps**2 ) )
1316 c = slamc3( half, -c )
1317 b = slamc3( half, c )
1318 c = slamc3( half, -b )
1319 b = slamc3( half, c )
1320 GO TO 10
1321 END IF
1322*+ END WHILE
1323*
1324 IF( a.LT.leps )
1325 $ leps = a
1326*
1327* Computation of EPS complete.
1328*
1329* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
1330* Keep dividing A by BETA until (gradual) underflow occurs. This
1331* is detected when we cannot recover the previous A.
1332*
1333 rbase = one / lbeta
1334 small = one
1335 DO 20 i = 1, 3
1336 small = slamc3( small*rbase, zero )
1337 20 CONTINUE
1338 a = slamc3( one, small )
1339 CALL slamc4( ngpmin, one, lbeta )
1340 CALL slamc4( ngnmin, -one, lbeta )
1341 CALL slamc4( gpmin, a, lbeta )
1342 CALL slamc4( gnmin, -a, lbeta )
1343 ieee = .false.
1344*
1345 IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) ) THEN
1346 IF( ngpmin.EQ.gpmin ) THEN
1347 lemin = ngpmin
1348* ( Non twos-complement machines, no gradual underflow;
1349* e.g., VAX )
1350 ELSE IF( ( gpmin-ngpmin ).EQ.3 ) THEN
1351 lemin = ngpmin - 1 + lt
1352 ieee = .true.
1353* ( Non twos-complement machines, with gradual underflow;
1354* e.g., IEEE standard followers )
1355 ELSE
1356 lemin = min( ngpmin, gpmin )
1357* ( A guess; no known machine )
1358 iwarn = .true.
1359 END IF
1360*
1361 ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) ) THEN
1362 IF( abs( ngpmin-ngnmin ).EQ.1 ) THEN
1363 lemin = max( ngpmin, ngnmin )
1364* ( Twos-complement machines, no gradual underflow;
1365* e.g., CYBER 205 )
1366 ELSE
1367 lemin = min( ngpmin, ngnmin )
1368* ( A guess; no known machine )
1369 iwarn = .true.
1370 END IF
1371*
1372 ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
1373 $ ( gpmin.EQ.gnmin ) ) THEN
1374 IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 ) THEN
1375 lemin = max( ngpmin, ngnmin ) - 1 + lt
1376* ( Twos-complement machines with gradual underflow;
1377* no known machine )
1378 ELSE
1379 lemin = min( ngpmin, ngnmin )
1380* ( A guess; no known machine )
1381 iwarn = .true.
1382 END IF
1383*
1384 ELSE
1385 lemin = min( ngpmin, ngnmin, gpmin, gnmin )
1386* ( A guess; no known machine )
1387 iwarn = .true.
1388 END IF
1389***
1390* Comment out this if block if EMIN is ok
1391 IF( iwarn ) THEN
1392 first = .true.
1393 WRITE( 6, fmt = 9999 )lemin
1394 END IF
1395***
1396*
1397* Assume IEEE arithmetic if we found denormalised numbers above,
1398* or if arithmetic seems to round in the IEEE style, determined
1399* in routine SLAMC1. A true IEEE machine should have both things
1400* true; however, faulty machines may have one or the other.
1401*
1402 ieee = ieee .OR. lieee1
1403*
1404* Compute RMIN by successive division by BETA. We could compute
1405* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
1406* this computation.
1407*
1408 lrmin = 1
1409 DO 30 i = 1, 1 - lemin
1410 lrmin = slamc3( lrmin*rbase, zero )
1411 30 CONTINUE
1412*
1413* Finally, call SLAMC5 to compute EMAX and RMAX.
1414*
1415 CALL slamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
1416 END IF
1417*
1418 beta = lbeta
1419 t = lt
1420 rnd = lrnd
1421 eps = leps
1422 emin = lemin
1423 rmin = lrmin
1424 emax = lemax
1425 rmax = lrmax
1426*
1427 RETURN
1428*
1429 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
1430 $ ' EMIN = ', i8, /
1431 $ ' If, after inspection, the value EMIN looks',
1432 $ ' acceptable please comment out ',
1433 $ / ' the IF block as marked within the code of routine',
1434 $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
1435*
1436* End of SLAMC2
1437*
slamc5
subroutine slamc5(beta, p, emin, ieee, emax, rmax)
Definition tools.f:1565
slamc1
subroutine slamc1(beta, t, rnd, ieee1)
Definition tools.f:997
slamc4
subroutine slamc4(emin, start, base)
Definition tools.f:1481

◆ slamc3()

real function slamc3 ( real a,
real b )

Definition at line 1442 of file tools.f.

1443*
1444* -- LAPACK auxiliary routine (version 2.0) --
1445* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1446* Courant Institute, Argonne National Lab, and Rice University
1447* October 31, 1992
1448*
1449* .. Scalar Arguments ..
1450 REAL A, B
1451* ..
1452*
1453* Purpose
1454* =======
1455*
1456* SLAMC3 is intended to force A and B to be stored prior to doing
1457* the addition of A and B , for use in situations where optimizers
1458* might hold one of these in a register.
1459*
1460* Arguments
1461* =========
1462*
1463* A, B (input) REAL
1464* The values A and B.
1465*
1466* =====================================================================
1467*
1468* .. Executable Statements ..
1469*
1470 slamc3 = a + b
1471*
1472 RETURN
1473*
1474* End of SLAMC3
1475*

◆ slamc4()

subroutine slamc4 ( integer emin,
real start,
integer base )

Definition at line 1480 of file tools.f.

1481*
1482* -- LAPACK auxiliary routine (version 2.0) --
1483* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1484* Courant Institute, Argonne National Lab, and Rice University
1485* October 31, 1992
1486*
1487* .. Scalar Arguments ..
1488 INTEGER BASE, EMIN
1489 REAL START
1490* ..
1491*
1492* Purpose
1493* =======
1494*
1495* SLAMC4 is a service routine for SLAMC2.
1496*
1497* Arguments
1498* =========
1499*
1500* EMIN (output) EMIN
1501* The minimum exponent before (gradual) underflow, computed by
1502* setting A = START and dividing by BASE until the previous A
1503* can not be recovered.
1504*
1505* START (input) REAL
1506* The starting point for determining EMIN.
1507*
1508* BASE (input) INTEGER
1509* The base of the machine.
1510*
1511* =====================================================================
1512*
1513* .. Local Scalars ..
1514 INTEGER I
1515 REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
1516* ..
1517* .. External Functions ..
1518 REAL SLAMC3
1519 EXTERNAL slamc3
1520* ..
1521* .. Executable Statements ..
1522*
1523 a = start
1524 one = 1
1525 rbase = one / base
1526 zero = 0
1527 emin = 1
1528 b1 = slamc3( a*rbase, zero )
1529 c1 = a
1530 c2 = a
1531 d1 = a
1532 d2 = a
1533*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
1534* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
1535 10 CONTINUE
1536 IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
1537 $ ( d2.EQ.a ) ) THEN
1538 emin = emin - 1
1539 a = b1
1540 b1 = slamc3( a / base, zero )
1541 c1 = slamc3( b1*base, zero )
1542 d1 = zero
1543 DO 20 i = 1, base
1544 d1 = d1 + b1
1545 20 CONTINUE
1546 b2 = slamc3( a*rbase, zero )
1547 c2 = slamc3( b2 / rbase, zero )
1548 d2 = zero
1549 DO 30 i = 1, base
1550 d2 = d2 + b2
1551 30 CONTINUE
1552 GO TO 10
1553 END IF
1554*+ END WHILE
1555*
1556 RETURN
1557*
1558* End of SLAMC4
1559*

◆ slamc5()

subroutine slamc5 ( integer beta,
integer p,
integer emin,
logical ieee,
integer emax,
real rmax )

Definition at line 1564 of file tools.f.

1565*
1566* -- LAPACK auxiliary routine (version 2.0) --
1567* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1568* Courant Institute, Argonne National Lab, and Rice University
1569* October 31, 1992
1570*
1571* .. Scalar Arguments ..
1572 LOGICAL IEEE
1573 INTEGER BETA, EMAX, EMIN, P
1574 REAL RMAX
1575* ..
1576*
1577* Purpose
1578* =======
1579*
1580* SLAMC5 attempts to compute RMAX, the largest machine floating-point
1581* number, without overflow. It assumes that EMAX + abs(EMIN) sum
1582* approximately to a power of 2. It will fail on machines where this
1583* assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
1584* EMAX = 28718). It will also fail if the value supplied for EMIN is
1585* too large (i.e. too close to zero), probably with overflow.
1586*
1587* Arguments
1588* =========
1589*
1590* BETA (input) INTEGER
1591* The base of floating-point arithmetic.
1592*
1593* P (input) INTEGER
1594* The number of base BETA digits in the mantissa of a
1595* floating-point value.
1596*
1597* EMIN (input) INTEGER
1598* The minimum exponent before (gradual) underflow.
1599*
1600* IEEE (input) LOGICAL
1601* A logical flag specifying whether or not the arithmetic
1602* system is thought to comply with the IEEE standard.
1603*
1604* EMAX (output) INTEGER
1605* The largest exponent before overflow
1606*
1607* RMAX (output) REAL
1608* The largest machine floating-point number.
1609*
1610* =====================================================================
1611*
1612* .. Parameters ..
1613 REAL ZERO, ONE
1614 parameter( zero = 0.0e0, one = 1.0e0 )
1615* ..
1616* .. Local Scalars ..
1617 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
1618 REAL OLDY, RECBAS, Y, Z
1619* ..
1620* .. External Functions ..
1621 REAL SLAMC3
1622 EXTERNAL slamc3
1623* ..
1624* .. Intrinsic Functions ..
1625 INTRINSIC mod
1626* ..
1627* .. Executable Statements ..
1628*
1629* First compute LEXP and UEXP, two powers of 2 that bound
1630* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
1631* approximately to the bound that is closest to abs(EMIN).
1632* (EMAX is the exponent of the required number RMAX).
1633*
1634 lexp = 1
1635 exbits = 1
1636 10 CONTINUE
1637 try = lexp*2
1638 IF( try.LE.( -emin ) ) THEN
1639 lexp = try
1640 exbits = exbits + 1
1641 GO TO 10
1642 END IF
1643 IF( lexp.EQ.-emin ) THEN
1644 uexp = lexp
1645 ELSE
1646 uexp = try
1647 exbits = exbits + 1
1648 END IF
1649*
1650* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
1651* than or equal to EMIN. EXBITS is the number of bits needed to
1652* store the exponent.
1653*
1654 IF( ( uexp+emin ).GT.( -lexp-emin ) ) THEN
1655 expsum = 2*lexp
1656 ELSE
1657 expsum = 2*uexp
1658 END IF
1659*
1660* EXPSUM is the exponent range, approximately equal to
1661* EMAX - EMIN + 1 .
1662*
1663 emax = expsum + emin - 1
1664 nbits = 1 + exbits + p
1665*
1666* NBITS is the total number of bits needed to store a
1667* floating-point number.
1668*
1669 IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) ) THEN
1670*
1671* Either there are an odd number of bits used to store a
1672* floating-point number, which is unlikely, or some bits are
1673* not used in the representation of numbers, which is possible,
1674* (e.g. Cray machines) or the mantissa has an implicit bit,
1675* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
1676* most likely. We have to assume the last alternative.
1677* If this is true, then we need to reduce EMAX by one because
1678* there must be some way of representing zero in an implicit-bit
1679* system. On machines like Cray, we are reducing EMAX by one
1680* unnecessarily.
1681*
1682 emax = emax - 1
1683 END IF
1684*
1685 IF( ieee ) THEN
1686*
1687* Assume we are on an IEEE machine which reserves one exponent
1688* for infinity and NaN.
1689*
1690 emax = emax - 1
1691 END IF
1692*
1693* Now create RMAX, the largest machine number, which should
1694* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
1695*
1696* First compute 1.0 - BETA**(-P), being careful that the
1697* result is less than 1.0 .
1698*
1699 recbas = one / beta
1700 z = beta - one
1701 y = zero
1702 DO 20 i = 1, p
1703 z = z*recbas
1704 IF( y.LT.one )
1705 $ oldy = y
1706 y = slamc3( y, z )
1707 20 CONTINUE
1708 IF( y.GE.one )
1709 $ y = oldy
1710*
1711* Now multiply by BETA**EMAX to get RMAX.
1712*
1713 DO 30 i = 1, emax
1714 y = slamc3( y*beta, zero )
1715 30 CONTINUE
1716*
1717 rmax = y
1718 RETURN
1719*
1720* End of SLAMC5
1721*

◆ slamch()

real function slamch ( character cmach)

Definition at line 866 of file tools.f.

867*
868* -- LAPACK auxiliary routine (version 2.0) --
869* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
870* Courant Institute, Argonne National Lab, and Rice University
871* October 31, 1992
872*
873* .. Scalar Arguments ..
874 CHARACTER CMACH
875* ..
876*
877* Purpose
878* =======
879*
880* SLAMCH determines single precision machine parameters.
881*
882* Arguments
883* =========
884*
885* CMACH (input) CHARACTER*1
886* Specifies the value to be returned by SLAMCH:
887* = 'E' or 'e', SLAMCH := eps
888* = 'S' or 's , SLAMCH := sfmin
889* = 'B' or 'b', SLAMCH := base
890* = 'P' or 'p', SLAMCH := eps*base
891* = 'N' or 'n', SLAMCH := t
892* = 'R' or 'r', SLAMCH := rnd
893* = 'M' or 'm', SLAMCH := emin
894* = 'U' or 'u', SLAMCH := rmin
895* = 'L' or 'l', SLAMCH := emax
896* = 'O' or 'o', SLAMCH := rmax
897*
898* where
899*
900* eps = relative machine precision
901* sfmin = safe minimum, such that 1/sfmin does not overflow
902* base = base of the machine
903* prec = eps*base
904* t = number of (base) digits in the mantissa
905* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
906* emin = minimum exponent before (gradual) underflow
907* rmin = underflow threshold - base**(emin-1)
908* emax = largest exponent before overflow
909* rmax = overflow threshold - (base**emax)*(1-eps)
910*
911* =====================================================================
912*
913* .. Parameters ..
914 REAL ONE, ZERO
915 parameter( one = 1.0e+0, zero = 0.0e+0 )
916* ..
917* .. Local Scalars ..
918 LOGICAL FIRST, LRND
919 INTEGER BETA, IMAX, IMIN, IT
920 REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
921 $ RND, SFMIN, SMALL, T
922* ..
923* .. External Functions ..
924 LOGICAL LSAME
925 EXTERNAL lsame
926* ..
927* .. External Subroutines ..
928 EXTERNAL slamc2
929* ..
930* .. Save statement ..
931 SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
932 $ emax, rmax, prec
933* ..
934* .. Data statements ..
935 DATA first / .true. /
936* ..
937* .. Executable Statements ..
938*
939 IF( first ) THEN
940 first = .false.
941 CALL slamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
942 base = beta
943 t = it
944 IF( lrnd ) THEN
945 rnd = one
946 eps = ( base**( 1-it ) ) / 2
947 ELSE
948 rnd = zero
949 eps = base**( 1-it )
950 END IF
951 prec = eps*base
952 emin = imin
953 emax = imax
954 sfmin = rmin
955 small = one / rmax
956 IF( small.GE.sfmin ) THEN
957*
958* Use SMALL plus a bit, to avoid the possibility of rounding
959* causing overflow when computing 1/sfmin.
960*
961 sfmin = small*( one+eps )
962 END IF
963 END IF
964*
965 IF( lsame( cmach, 'E' ) ) THEN
966 rmach = eps
967 ELSE IF( lsame( cmach, 'S' ) ) THEN
968 rmach = sfmin
969 ELSE IF( lsame( cmach, 'B' ) ) THEN
970 rmach = base
971 ELSE IF( lsame( cmach, 'P' ) ) THEN
972 rmach = prec
973 ELSE IF( lsame( cmach, 'N' ) ) THEN
974 rmach = t
975 ELSE IF( lsame( cmach, 'R' ) ) THEN
976 rmach = rnd
977 ELSE IF( lsame( cmach, 'M' ) ) THEN
978 rmach = emin
979 ELSE IF( lsame( cmach, 'U' ) ) THEN
980 rmach = rmin
981 ELSE IF( lsame( cmach, 'L' ) ) THEN
982 rmach = emax
983 ELSE IF( lsame( cmach, 'O' ) ) THEN
984 rmach = rmax
985 END IF
986*
987 slamch = rmach
988 RETURN
989*
990* End of SLAMCH
991*
slamch
real function slamch(cmach)
Definition tools.f:867
slamc2
subroutine slamc2(beta, t, rnd, eps, emin, rmin, emax, rmax)
Definition tools.f:1184

◆ zlarnd()

double complex function zlarnd ( integer idist,
integer, dimension( 4 ) iseed )

Definition at line 1898 of file tools.f.

1899*
1900* -- LAPACK auxiliary routine (version 2.0) --
1901* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1902* Courant Institute, Argonne National Lab, and Rice University
1903* June 30, 1994
1904*
1905* .. Scalar Arguments ..
1906 INTEGER IDIST
1907* ..
1908* .. Array Arguments ..
1909 INTEGER ISEED( 4 )
1910* ..
1911*
1912* Purpose
1913* =======
1914*
1915* ZLARND returns a random complex number from a uniform or normal
1916* distribution.
1917*
1918* Arguments
1919* =========
1920*
1921* IDIST (input) INTEGER
1922* Specifies the distribution of the random numbers:
1923* = 1: real and imaginary parts each uniform (0,1)
1924* = 2: real and imaginary parts each uniform (-1,1)
1925* = 3: real and imaginary parts each normal (0,1)
1926* = 4: uniformly distributed on the disc abs(z) <= 1
1927* = 5: uniformly distributed on the circle abs(z) = 1
1928*
1929* ISEED (input/output) INTEGER array, dimension (4)
1930* On entry, the seed of the random number generator; the array
1931* elements must be between 0 and 4095, and ISEED(4) must be
1932* odd.
1933* On exit, the seed is updated.
1934*
1935* Further Details
1936* ===============
1937*
1938* This routine calls the auxiliary routine DLARAN to generate a random
1939* real number from a uniform (0,1) distribution. The Box-Muller method
1940* is used to transform numbers from a uniform to a normal distribution.
1941*
1942* =====================================================================
1943*
1944* .. Parameters ..
1945 DOUBLE PRECISION ZERO, ONE, TWO
1946 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
1947 DOUBLE PRECISION TWOPI
1948 parameter( twopi = 6.2831853071795864769252867663d+0 )
1949* ..
1950* .. Local Scalars ..
1951 DOUBLE PRECISION T1, T2
1952* ..
1953* .. External Functions ..
1954 DOUBLE PRECISION DLARAN
1955 EXTERNAL dlaran
1956* ..
1957* .. Intrinsic Functions ..
1958 INTRINSIC dcmplx, exp, log, sqrt
1959* ..
1960* .. Executable Statements ..
1961*
1962* Generate a pair of real random numbers from a uniform (0,1)
1963* distribution
1964*
1965 t1 = dlaran( iseed )
1966 t2 = dlaran( iseed )
1967*
1968 IF( idist.EQ.1 ) THEN
1969*
1970* real and imaginary parts each uniform (0,1)
1971*
1972 zlarnd = dcmplx( t1, t2 )
1973 ELSE IF( idist.EQ.2 ) THEN
1974*
1975* real and imaginary parts each uniform (-1,1)
1976*
1977 zlarnd = dcmplx( two*t1-one, two*t2-one )
1978 ELSE IF( idist.EQ.3 ) THEN
1979*
1980* real and imaginary parts each normal (0,1)
1981*
1982 zlarnd = sqrt( -two*log( t1 ) )*exp( dcmplx( zero, twopi*t2 ) )
1983 ELSE IF( idist.EQ.4 ) THEN
1984*
1985* uniform distribution on the unit disc abs(z) <= 1
1986*
1987 zlarnd = sqrt( t1 )*exp( dcmplx( zero, twopi*t2 ) )
1988 ELSE IF( idist.EQ.5 ) THEN
1989*
1990* uniform distribution on the unit circle abs(z) = 1
1991*
1992 zlarnd = exp( dcmplx( zero, twopi*t2 ) )
1993 END IF
1994 RETURN
1995*
1996* End of ZLARND
1997*
zlarnd
double complex function zlarnd(idist, iseed)
Definition tools.f:1899