Bug Summary

File:build/include/omc/c/external_solvers/dlamch.c
Warning:line 159, column 13
Assigned value is garbage or undefined

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/* dlamch.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11*/
12
13#include "f2c.h"
14#include "blaswrap.h"
15#include <stdlib.h>
16#include <math.h>
17
18/* Table of constant values */
19
20static integer c__1 = 1;
21static doublereal c_b32 = 0.;
22
23doublereal dlamch_(char *cmach)
24{
25 /* Initialized data */
26
27 static logical first = TRUE_(1);
28
29 /* System generated locals */
30 integer i__1;
31 doublereal ret_val;
32
33 /* Builtin functions */
34 double pow_di(doublereal *, integer *);
35
36 /* Local variables */
37 static doublereal t;
38 integer it;
39 static doublereal rnd, eps, base;
40 integer beta;
41 static doublereal emin, prec, emax;
42 integer imin, imax;
43 logical lrnd;
44 static doublereal rmin, rmax;
45 doublereal rmach;
1
'rmach' declared without an initial value
46 extern logical lsame_(char *, char *);
47 doublereal small;
48 static doublereal sfmin;
49 extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
50 doublereal *, integer *, doublereal *, integer *, doublereal *);
51
52
53/* -- LAPACK auxiliary routine (version 3.2) -- */
54/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
55/* November 2006 */
56
57/* .. Scalar Arguments .. */
58/* .. */
59
60/* Purpose */
61/* ======= */
62
63/* DLAMCH determines double precision machine parameters. */
64
65/* Arguments */
66/* ========= */
67
68/* CMACH (input) CHARACTER*1 */
69/* Specifies the value to be returned by DLAMCH: */
70/* = 'E' or 'e', DLAMCH := eps */
71/* = 'S' or 's , DLAMCH := sfmin */
72/* = 'B' or 'b', DLAMCH := base */
73/* = 'P' or 'p', DLAMCH := eps*base */
74/* = 'N' or 'n', DLAMCH := t */
75/* = 'R' or 'r', DLAMCH := rnd */
76/* = 'M' or 'm', DLAMCH := emin */
77/* = 'U' or 'u', DLAMCH := rmin */
78/* = 'L' or 'l', DLAMCH := emax */
79/* = 'O' or 'o', DLAMCH := rmax */
80
81/* where */
82
83/* eps = relative machine precision */
84/* sfmin = safe minimum, such that 1/sfmin does not overflow */
85/* base = base of the machine */
86/* prec = eps*base */
87/* t = number of (base) digits in the mantissa */
88/* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
89/* emin = minimum exponent before (gradual) underflow */
90/* rmin = underflow threshold - base**(emin-1) */
91/* emax = largest exponent before overflow */
92/* rmax = overflow threshold - (base**emax)*(1-eps) */
93
94/* ===================================================================== */
95
96/* .. Parameters .. */
97/* .. */
98/* .. Local Scalars .. */
99/* .. */
100/* .. External Functions .. */
101/* .. */
102/* .. External Subroutines .. */
103/* .. */
104/* .. Save statement .. */
105/* .. */
106/* .. Data statements .. */
107/* .. */
108/* .. Executable Statements .. */
109
110 if (first) {
2
Taking true branch
111 dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
112 base = (doublereal) beta;
113 t = (doublereal) it;
114 if (lrnd) {
3
Assuming 'lrnd' is 0
4
Taking false branch
115 rnd = 1.;
116 i__1 = 1 - it;
117 eps = pow_di(&base, &i__1) / 2;
118 } else {
119 rnd = 0.;
120 i__1 = 1 - it;
121 eps = pow_di(&base, &i__1);
122 }
123 prec = eps * base;
124 emin = (doublereal) imin;
125 emax = (doublereal) imax;
126 sfmin = rmin;
127 small = 1. / rmax;
128 if (small >= sfmin) {
5
Taking false branch
129
130/* Use SMALL plus a bit, to avoid the possibility of rounding */
131/* causing overflow when computing 1/sfmin. */
132
133 sfmin = small * (eps + 1.);
134 }
135 }
136
137 if (lsame_(cmach, "E")) {
6
Assuming the condition is false
7
Taking false branch
138 rmach = eps;
139 } else if (lsame_(cmach, "S")) {
8
Assuming the condition is false
9
Taking false branch
140 rmach = sfmin;
141 } else if (lsame_(cmach, "B")) {
10
Assuming the condition is false
11
Taking false branch
142 rmach = base;
143 } else if (lsame_(cmach, "P")) {
12
Assuming the condition is false
13
Taking false branch
144 rmach = prec;
145 } else if (lsame_(cmach, "N")) {
14
Assuming the condition is false
15
Taking false branch
146 rmach = t;
147 } else if (lsame_(cmach, "R")) {
16
Assuming the condition is false
17
Taking false branch
148 rmach = rnd;
149 } else if (lsame_(cmach, "M")) {
18
Assuming the condition is false
19
Taking false branch
150 rmach = emin;
151 } else if (lsame_(cmach, "U")) {
20
Assuming the condition is false
21
Taking false branch
152 rmach = rmin;
153 } else if (lsame_(cmach, "L")) {
22
Assuming the condition is false
23
Taking false branch
154 rmach = emax;
155 } else if (lsame_(cmach, "O")) {
24
Assuming the condition is false
25
Taking false branch
156 rmach = rmax;
157 }
158
159 ret_val = rmach;
26
Assigned value is garbage or undefined
160 first = FALSE_(0);
161 return ret_val;
162
163/* End of DLAMCH */
164
165} /* dlamch_ */
166
167
168/* *********************************************************************** */
169
170/* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical
171 *ieee1)
172{
173 /* Initialized data */
174
175 static logical first = TRUE_(1);
176
177 /* System generated locals */
178 doublereal d__1, d__2;
179
180 /* Local variables */
181 doublereal a, b, c__, f, t1, t2;
182 static integer lt;
183 doublereal one, qtr;
184 static logical lrnd;
185 static integer lbeta;
186 doublereal savec;
187 extern doublereal dlamc3_(doublereal *, doublereal *);
188 static logical lieee1;
189
190
191/* -- LAPACK auxiliary routine (version 3.2) -- */
192/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
193/* November 2006 */
194
195/* .. Scalar Arguments .. */
196/* .. */
197
198/* Purpose */
199/* ======= */
200
201/* DLAMC1 determines the machine parameters given by BETA, T, RND, and */
202/* IEEE1. */
203
204/* Arguments */
205/* ========= */
206
207/* BETA (output) INTEGER */
208/* The base of the machine. */
209
210/* T (output) INTEGER */
211/* The number of ( BETA ) digits in the mantissa. */
212
213/* RND (output) LOGICAL */
214/* Specifies whether proper rounding ( RND = .TRUE. ) or */
215/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
216/* be a reliable guide to the way in which the machine performs */
217/* its arithmetic. */
218
219/* IEEE1 (output) LOGICAL */
220/* Specifies whether rounding appears to be done in the IEEE */
221/* 'round to nearest' style. */
222
223/* Further Details */
224/* =============== */
225
226/* The routine is based on the routine ENVRON by Malcolm and */
227/* incorporates suggestions by Gentleman and Marovich. See */
228
229/* Malcolm M. A. (1972) Algorithms to reveal properties of */
230/* floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
231
232/* Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
233/* that reveal properties of floating point arithmetic units. */
234/* Comms. of the ACM, 17, 276-277. */
235
236/* ===================================================================== */
237
238/* .. Local Scalars .. */
239/* .. */
240/* .. External Functions .. */
241/* .. */
242/* .. Save statement .. */
243/* .. */
244/* .. Data statements .. */
245/* .. */
246/* .. Executable Statements .. */
247
248 if (first) {
249 one = 1.;
250
251/* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
252/* IEEE1, T and RND. */
253
254/* Throughout this routine we use the function DLAMC3 to ensure */
255/* that relevant values are stored and not held in registers, or */
256/* are not affected by optimizers. */
257
258/* Compute a = 2.0**m with the smallest positive integer m such */
259/* that */
260
261/* fl( a + 1.0 ) = a. */
262
263 a = 1.;
264 c__ = 1.;
265
266/* + WHILE( C.EQ.ONE )LOOP */
267L10:
268 if (c__ == one) {
269 a *= 2;
270 c__ = dlamc3_(&a, &one);
271 d__1 = -a;
272 c__ = dlamc3_(&c__, &d__1);
273 goto L10;
274 }
275/* + END WHILE */
276
277/* Now compute b = 2.0**m with the smallest positive integer m */
278/* such that */
279
280/* fl( a + b ) .gt. a. */
281
282 b = 1.;
283 c__ = dlamc3_(&a, &b);
284
285/* + WHILE( C.EQ.A )LOOP */
286L20:
287 if (c__ == a) {
288 b *= 2;
289 c__ = dlamc3_(&a, &b);
290 goto L20;
291 }
292/* + END WHILE */
293
294/* Now compute the base. a and c are neighbouring floating point */
295/* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
296/* their difference is beta. Adding 0.25 to c is to ensure that it */
297/* is truncated to beta and not ( beta - 1 ). */
298
299 qtr = one / 4;
300 savec = c__;
301 d__1 = -a;
302 c__ = dlamc3_(&c__, &d__1);
303 lbeta = (integer) (c__ + qtr);
304
305/* Now determine whether rounding or chopping occurs, by adding a */
306/* bit less than beta/2 and a bit more than beta/2 to a. */
307
308 b = (doublereal) lbeta;
309 d__1 = b / 2;
310 d__2 = -b / 100;
311 f = dlamc3_(&d__1, &d__2);
312 c__ = dlamc3_(&f, &a);
313 if (c__ == a) {
314 lrnd = TRUE_(1);
315 } else {
316 lrnd = FALSE_(0);
317 }
318 d__1 = b / 2;
319 d__2 = b / 100;
320 f = dlamc3_(&d__1, &d__2);
321 c__ = dlamc3_(&f, &a);
322 if (lrnd && c__ == a) {
323 lrnd = FALSE_(0);
324 }
325
326/* Try and decide whether rounding is done in the IEEE 'round to */
327/* nearest' style. B/2 is half a unit in the last place of the two */
328/* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
329/* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
330/* A, but adding B/2 to SAVEC should change SAVEC. */
331
332 d__1 = b / 2;
333 t1 = dlamc3_(&d__1, &a);
334 d__1 = b / 2;
335 t2 = dlamc3_(&d__1, &savec);
336 lieee1 = t1 == a && t2 > savec && lrnd;
337
338/* Now find the mantissa, t. It should be the integer part of */
339/* log to the base beta of a, however it is safer to determine t */
340/* by powering. So we find t as the smallest positive integer for */
341/* which */
342
343/* fl( beta**t + 1.0 ) = 1.0. */
344
345 lt = 0;
346 a = 1.;
347 c__ = 1.;
348
349/* + WHILE( C.EQ.ONE )LOOP */
350L30:
351 if (c__ == one) {
352 ++lt;
353 a *= lbeta;
354 c__ = dlamc3_(&a, &one);
355 d__1 = -a;
356 c__ = dlamc3_(&c__, &d__1);
357 goto L30;
358 }
359/* + END WHILE */
360
361 }
362
363 *beta = lbeta;
364 *t = lt;
365 *rnd = lrnd;
366 *ieee1 = lieee1;
367 first = FALSE_(0);
368 return 0;
369
370/* End of DLAMC1 */
371
372} /* dlamc1_ */
373
374
375/* *********************************************************************** */
376
377/* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd,
378 doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
379 doublereal *rmax)
380{
381 /* Initialized data */
382
383 static logical first = TRUE_(1);
384 static logical iwarn = FALSE_(0);
385
386 /* Format strings */
387 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
388 "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
389 "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
390 " the IF block as marked within the code of routine\002,\002 DLAM"
391 "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
392
393 /* System generated locals */
394 integer i__1;
395 doublereal d__1, d__2, d__3, d__4, d__5;
396
397 /* Builtin functions */
398 double pow_di(doublereal *, integer *);
399 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
400
401 /* Local variables */
402 doublereal a, b, c__;
403 integer i__;
404 static integer lt;
405 doublereal one, two;
406 logical ieee;
407 doublereal half;
408 logical lrnd;
409 static doublereal leps;
410 doublereal zero;
411 static integer lbeta;
412 doublereal rbase;
413 static integer lemin, lemax;
414 integer gnmin;
415 doublereal small;
416 integer gpmin;
417 doublereal third;
418 static doublereal lrmin, lrmax;
419 doublereal sixth;
420 extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *,
421 logical *);
422 extern doublereal dlamc3_(doublereal *, doublereal *);
423 logical lieee1;
424 extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *),
425 dlamc5_(integer *, integer *, integer *, logical *, integer *,
426 doublereal *);
427 integer ngnmin, ngpmin;
428
429 /* Fortran I/O blocks */
430 static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
431
432
433
434/* -- LAPACK auxiliary routine (version 3.2) -- */
435/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
436/* November 2006 */
437
438/* .. Scalar Arguments .. */
439/* .. */
440
441/* Purpose */
442/* ======= */
443
444/* DLAMC2 determines the machine parameters specified in its argument */
445/* list. */
446
447/* Arguments */
448/* ========= */
449
450/* BETA (output) INTEGER */
451/* The base of the machine. */
452
453/* T (output) INTEGER */
454/* The number of ( BETA ) digits in the mantissa. */
455
456/* RND (output) LOGICAL */
457/* Specifies whether proper rounding ( RND = .TRUE. ) or */
458/* chopping ( RND = .FALSE. ) occurs in addition. This may not */
459/* be a reliable guide to the way in which the machine performs */
460/* its arithmetic. */
461
462/* EPS (output) DOUBLE PRECISION */
463/* The smallest positive number such that */
464
465/* fl( 1.0 - EPS ) .LT. 1.0, */
466
467/* where fl denotes the computed value. */
468
469/* EMIN (output) INTEGER */
470/* The minimum exponent before (gradual) underflow occurs. */
471
472/* RMIN (output) DOUBLE PRECISION */
473/* The smallest normalized number for the machine, given by */
474/* BASE**( EMIN - 1 ), where BASE is the floating point value */
475/* of BETA. */
476
477/* EMAX (output) INTEGER */
478/* The maximum exponent before overflow occurs. */
479
480/* RMAX (output) DOUBLE PRECISION */
481/* The largest positive number for the machine, given by */
482/* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
483/* value of BETA. */
484
485/* Further Details */
486/* =============== */
487
488/* The computation of EPS is based on a routine PARANOIA by */
489/* W. Kahan of the University of California at Berkeley. */
490
491/* ===================================================================== */
492
493/* .. Local Scalars .. */
494/* .. */
495/* .. External Functions .. */
496/* .. */
497/* .. External Subroutines .. */
498/* .. */
499/* .. Intrinsic Functions .. */
500/* .. */
501/* .. Save statement .. */
502/* .. */
503/* .. Data statements .. */
504/* .. */
505/* .. Executable Statements .. */
506
507 if (first) {
508 zero = 0.;
509 one = 1.;
510 two = 2.;
511
512/* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
513/* BETA, T, RND, EPS, EMIN and RMIN. */
514
515/* Throughout this routine we use the function DLAMC3 to ensure */
516/* that relevant values are stored and not held in registers, or */
517/* are not affected by optimizers. */
518
519/* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
520
521 dlamc1_(&lbeta, &lt, &lrnd, &lieee1);
522
523/* Start to find EPS. */
524
525 b = (doublereal) lbeta;
526 i__1 = -lt;
527 a = pow_di(&b, &i__1);
528 leps = a;
529
530/* Try some tricks to see whether or not this is the correct EPS. */
531
532 b = two / 3;
533 half = one / 2;
534 d__1 = -half;
535 sixth = dlamc3_(&b, &d__1);
536 third = dlamc3_(&sixth, &sixth);
537 d__1 = -half;
538 b = dlamc3_(&third, &d__1);
539 b = dlamc3_(&b, &sixth);
540 b = fabs(b);
541 if (b < leps) {
542 b = leps;
543 }
544
545 leps = 1.;
546
547/* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
548L10:
549 if (leps > b && b > zero) {
550 leps = b;
551 d__1 = half * leps;
552/* Computing 5th power */
553 d__3 = two, d__4 = d__3, d__3 *= d__3;
554/* Computing 2nd power */
555 d__5 = leps;
556 d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
557 c__ = dlamc3_(&d__1, &d__2);
558 d__1 = -c__;
559 c__ = dlamc3_(&half, &d__1);
560 b = dlamc3_(&half, &c__);
561 d__1 = -b;
562 c__ = dlamc3_(&half, &d__1);
563 b = dlamc3_(&half, &c__);
564 goto L10;
565 }
566/* + END WHILE */
567
568 if (a < leps) {
569 leps = a;
570 }
571
572/* Computation of EPS complete. */
573
574/* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
575/* Keep dividing A by BETA until (gradual) underflow occurs. This */
576/* is detected when we cannot recover the previous A. */
577
578 rbase = one / lbeta;
579 small = one;
580 for (i__ = 1; i__ <= 3; ++i__) {
581 d__1 = small * rbase;
582 small = dlamc3_(&d__1, &zero);
583/* L20: */
584 }
585 a = dlamc3_(&one, &small);
586 dlamc4_(&ngpmin, &one, &lbeta);
587 d__1 = -one;
588 dlamc4_(&ngnmin, &d__1, &lbeta);
589 dlamc4_(&gpmin, &a, &lbeta);
590 d__1 = -a;
591 dlamc4_(&gnmin, &d__1, &lbeta);
592 ieee = FALSE_(0);
593
594 if (ngpmin == ngnmin && gpmin == gnmin) {
595 if (ngpmin == gpmin) {
596 lemin = ngpmin;
597/* ( Non twos-complement machines, no gradual underflow; */
598/* e.g., VAX ) */
599 } else if (gpmin - ngpmin == 3) {
600 lemin = ngpmin - 1 + lt;
601 ieee = TRUE_(1);
602/* ( Non twos-complement machines, with gradual underflow; */
603/* e.g., IEEE standard followers ) */
604 } else {
605 lemin = lmin(ngpmin,gpmin);
606/* ( A guess; no known machine ) */
607 iwarn = TRUE_(1);
608 }
609
610 } else if (ngpmin == gpmin && ngnmin == gnmin) {
611 if ((i__1 = ngpmin - ngnmin, labs(i__1)) == 1) {
612 lemin = lmax(ngpmin,ngnmin);
613/* ( Twos-complement machines, no gradual underflow; */
614/* e.g., CYBER 205 ) */
615 } else {
616 lemin = lmin(ngpmin,ngnmin);
617/* ( A guess; no known machine ) */
618 iwarn = TRUE_(1);
619 }
620
621 } else if ((i__1 = ngpmin - ngnmin, labs(i__1)) == 1 && gpmin == gnmin)
622 {
623 if (gpmin - lmin(ngpmin,ngnmin) == 3) {
624 lemin = lmax(ngpmin,ngnmin) - 1 + lt;
625/* ( Twos-complement machines with gradual underflow; */
626/* no known machine ) */
627 } else {
628 lemin = lmin(ngpmin,ngnmin);
629/* ( A guess; no known machine ) */
630 iwarn = TRUE_(1);
631 }
632
633 } else {
634/* Computing MIN */
635 i__1 = lmin(ngpmin,ngnmin), i__1 = lmin(i__1,gpmin);
636 lemin = lmin(i__1,gnmin);
637/* ( A guess; no known machine ) */
638 iwarn = TRUE_(1);
639 }
640 first = FALSE_(0);
641/* ** */
642/* Comment out this if block if EMIN is ok */
643 if (iwarn) {
644 abort();
645/*
646 first = TRUE_;
647 s_wsfe(&io___58);
648 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
649 e_wsfe();
650*/
651 }
652/* ** */
653
654/* Assume IEEE arithmetic if we found denormalised numbers above, */
655/* or if arithmetic seems to round in the IEEE style, determined */
656/* in routine DLAMC1. A true IEEE machine should have both things */
657/* true; however, faulty machines may have one or the other. */
658
659 ieee = ieee || lieee1;
660
661/* Compute RMIN by successive division by BETA. We could compute */
662/* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
663/* this computation. */
664
665 lrmin = 1.;
666 i__1 = 1 - lemin;
667 for (i__ = 1; i__ <= i__1; ++i__) {
668 d__1 = lrmin * rbase;
669 lrmin = dlamc3_(&d__1, &zero);
670/* L30: */
671 }
672
673/* Finally, call DLAMC5 to compute EMAX and RMAX. */
674
675 dlamc5_(&lbeta, &lt, &lemin, &ieee, &lemax, &lrmax);
676 }
677
678 *beta = lbeta;
679 *t = lt;
680 *rnd = lrnd;
681 *eps = leps;
682 *emin = lemin;
683 *rmin = lrmin;
684 *emax = lemax;
685 *rmax = lrmax;
686
687 return 0;
688
689
690/* End of DLAMC2 */
691
692} /* dlamc2_ */
693
694
695/* *********************************************************************** */
696
697doublereal dlamc3_(doublereal *a, doublereal *b)
698{
699 /* System generated locals */
700 doublereal ret_val;
701
702
703/* -- LAPACK auxiliary routine (version 3.2) -- */
704/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
705/* November 2006 */
706
707/* .. Scalar Arguments .. */
708/* .. */
709
710/* Purpose */
711/* ======= */
712
713/* DLAMC3 is intended to force A and B to be stored prior to doing */
714/* the addition of A and B , for use in situations where optimizers */
715/* might hold one of these in a register. */
716
717/* Arguments */
718/* ========= */
719
720/* A (input) DOUBLE PRECISION */
721/* B (input) DOUBLE PRECISION */
722/* The values A and B. */
723
724/* ===================================================================== */
725
726/* .. Executable Statements .. */
727
728 ret_val = *a + *b;
729
730 return ret_val;
731
732/* End of DLAMC3 */
733
734} /* dlamc3_ */
735
736
737/* *********************************************************************** */
738
739/* Subroutine */ int dlamc4_(integer *emin, doublereal *start, integer *base)
740{
741 /* System generated locals */
742 integer i__1;
743 doublereal d__1;
744
745 /* Local variables */
746 doublereal a;
747 integer i__;
748 doublereal b1, b2, c1, c2, d1, d2, one, zero, rbase;
749 extern doublereal dlamc3_(doublereal *, doublereal *);
750
751
752/* -- LAPACK auxiliary routine (version 3.2) -- */
753/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
754/* November 2006 */
755
756/* .. Scalar Arguments .. */
757/* .. */
758
759/* Purpose */
760/* ======= */
761
762/* DLAMC4 is a service routine for DLAMC2. */
763
764/* Arguments */
765/* ========= */
766
767/* EMIN (output) INTEGER */
768/* The minimum exponent before (gradual) underflow, computed by */
769/* setting A = START and dividing by BASE until the previous A */
770/* can not be recovered. */
771
772/* START (input) DOUBLE PRECISION */
773/* The starting point for determining EMIN. */
774
775/* BASE (input) INTEGER */
776/* The base of the machine. */
777
778/* ===================================================================== */
779
780/* .. Local Scalars .. */
781/* .. */
782/* .. External Functions .. */
783/* .. */
784/* .. Executable Statements .. */
785
786 a = *start;
787 one = 1.;
788 rbase = one / *base;
789 zero = 0.;
790 *emin = 1;
791 d__1 = a * rbase;
792 b1 = dlamc3_(&d__1, &zero);
793 c1 = a;
794 c2 = a;
795 d1 = a;
796 d2 = a;
797/* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
798/* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
799L10:
800 if (c1 == a && c2 == a && d1 == a && d2 == a) {
801 --(*emin);
802 a = b1;
803 d__1 = a / *base;
804 b1 = dlamc3_(&d__1, &zero);
805 d__1 = b1 * *base;
806 c1 = dlamc3_(&d__1, &zero);
807 d1 = zero;
808 i__1 = *base;
809 for (i__ = 1; i__ <= i__1; ++i__) {
810 d1 += b1;
811/* L20: */
812 }
813 d__1 = a * rbase;
814 b2 = dlamc3_(&d__1, &zero);
815 d__1 = b2 / rbase;
816 c2 = dlamc3_(&d__1, &zero);
817 d2 = zero;
818 i__1 = *base;
819 for (i__ = 1; i__ <= i__1; ++i__) {
820 d2 += b2;
821/* L30: */
822 }
823 goto L10;
824 }
825/* + END WHILE */
826
827 return 0;
828
829/* End of DLAMC4 */
830
831} /* dlamc4_ */
832
833
834/* *********************************************************************** */
835
836/* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin,
837 logical *ieee, integer *emax, doublereal *rmax)
838{
839 /* System generated locals */
840 integer i__1;
841 doublereal d__1;
842
843 /* Local variables */
844 integer i__;
845 doublereal y, z__;
846 integer try__, lexp;
847 doublereal oldy;
848 integer uexp, nbits;
849 extern doublereal dlamc3_(doublereal *, doublereal *);
850 doublereal recbas;
851 integer exbits, expsum;
852
853
854/* -- LAPACK auxiliary routine (version 3.2) -- */
855/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
856/* November 2006 */
857
858/* .. Scalar Arguments .. */
859/* .. */
860
861/* Purpose */
862/* ======= */
863
864/* DLAMC5 attempts to compute RMAX, the largest machine floating-point */
865/* number, without overflow. It assumes that EMAX + abs(EMIN) sum */
866/* approximately to a power of 2. It will fail on machines where this */
867/* assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
868/* EMAX = 28718). It will also fail if the value supplied for EMIN is */
869/* too large (i.e. too close to zero), probably with overflow. */
870
871/* Arguments */
872/* ========= */
873
874/* BETA (input) INTEGER */
875/* The base of floating-point arithmetic. */
876
877/* P (input) INTEGER */
878/* The number of base BETA digits in the mantissa of a */
879/* floating-point value. */
880
881/* EMIN (input) INTEGER */
882/* The minimum exponent before (gradual) underflow. */
883
884/* IEEE (input) LOGICAL */
885/* A logical flag specifying whether or not the arithmetic */
886/* system is thought to comply with the IEEE standard. */
887
888/* EMAX (output) INTEGER */
889/* The largest exponent before overflow */
890
891/* RMAX (output) DOUBLE PRECISION */
892/* The largest machine floating-point number. */
893
894/* ===================================================================== */
895
896/* .. Parameters .. */
897/* .. */
898/* .. Local Scalars .. */
899/* .. */
900/* .. External Functions .. */
901/* .. */
902/* .. Intrinsic Functions .. */
903/* .. */
904/* .. Executable Statements .. */
905
906/* First compute LEXP and UEXP, two powers of 2 that bound */
907/* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
908/* approximately to the bound that is closest to abs(EMIN). */
909/* (EMAX is the exponent of the required number RMAX). */
910
911 lexp = 1;
912 exbits = 1;
913L10:
914 try__ = lexp << 1;
915 if (try__ <= -(*emin)) {
916 lexp = try__;
917 ++exbits;
918 goto L10;
919 }
920 if (lexp == -(*emin)) {
921 uexp = lexp;
922 } else {
923 uexp = try__;
924 ++exbits;
925 }
926
927/* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
928/* than or equal to EMIN. EXBITS is the number of bits needed to */
929/* store the exponent. */
930
931 if (uexp + *emin > -lexp - *emin) {
932 expsum = lexp << 1;
933 } else {
934 expsum = uexp << 1;
935 }
936
937/* EXPSUM is the exponent range, approximately equal to */
938/* EMAX - EMIN + 1 . */
939
940 *emax = expsum + *emin - 1;
941 nbits = exbits + 1 + *p;
942
943/* NBITS is the total number of bits needed to store a */
944/* floating-point number. */
945
946 if (nbits % 2 == 1 && *beta == 2) {
947
948/* Either there are an odd number of bits used to store a */
949/* floating-point number, which is unlikely, or some bits are */
950/* not used in the representation of numbers, which is possible, */
951/* (e.g. Cray machines) or the mantissa has an implicit bit, */
952/* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
953/* most likely. We have to assume the last alternative. */
954/* If this is true, then we need to reduce EMAX by one because */
955/* there must be some way of representing zero in an implicit-bit */
956/* system. On machines like Cray, we are reducing EMAX by one */
957/* unnecessarily. */
958
959 --(*emax);
960 }
961
962 if (*ieee) {
963
964/* Assume we are on an IEEE machine which reserves one exponent */
965/* for infinity and NaN. */
966
967 --(*emax);
968 }
969
970/* Now create RMAX, the largest machine number, which should */
971/* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
972
973/* First compute 1.0 - BETA**(-P), being careful that the */
974/* result is less than 1.0 . */
975
976 recbas = 1. / *beta;
977 z__ = *beta - 1.;
978 y = 0.;
979 i__1 = *p;
980 for (i__ = 1; i__ <= i__1; ++i__) {
981 z__ *= recbas;
982 if (y < 1.) {
983 oldy = y;
984 }
985 y = dlamc3_(&y, &z__);
986/* L20: */
987 }
988 if (y >= 1.) {
989 y = oldy;
990 }
991
992/* Now multiply by BETA**EMAX to get RMAX. */
993
994 i__1 = *emax;
995 for (i__ = 1; i__ <= i__1; ++i__) {
996 d__1 = y * *beta;
997 y = dlamc3_(&d__1, &c_b32);
998/* L30: */
999 }
1000
1001 *rmax = y;
1002 return 0;
1003
1004/* End of DLAMC5 */
1005
1006} /* dlamc5_ */