rts-lapack.c
1 //!@file rts-lapack.c
2 //!@author J. Marcel van der Veer
3 //
4 // !@section Copyright
5 //
6 // This file is part of VIF - vintage FORTRAN compiler.
7 // Copyright 2020-2025 J. Marcel van der Veer <algol68g@xs4all.nl>.
8 //
9 //! @section License
10 //
11 // This program is free software; you can redistribute it and/or modify it
12 // under the terms of the GNU General Public License as published by the
13 // Free Software Foundation; either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful, but
17 // WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 // more details. You should have received a copy of the GNU General Public
20 // License along with this program. If not, see <http://www.gnu.org/licenses/>.
21
22 //! @section Synopsis
23 //!
24 //! Runtime support for LAPACK subprograms.
25
26 // Auxilliary routines for LAPACK.
27
28 #include "vif.h"
29
30 logical_4 _lsame (char _p_, char _p_);
31 logical_4 _lsamen (int_4 *, char _p_, char _p_);
32 real_8 _dlamch (char _p_ cmach);
33 int_4 _dlamc1 (int_4 *, int_4 *, logical_4 *, logical_4 *);
34 int_4 _dlamc2 (int_4 *, int_4 *, logical_4 *, real_8 *, int_4 *, real_8 *, int_4 *, real_8 *);
35 real_8 _dlamc3 (real_8 *, real_8 *);
36 int_4 _dlamc4 (int_4 *, real_8 *, int_4 *);
37 int_4 _dlamc5 (int_4 *, int_4 *, int_4 *, logical_4 *, int_4 *, real_8 *);
38
39 static real_8 zero_arg = 0.0;
40
41 logical_4 _lsame (char _p_ ca, char _p_ cb)
42 {
43 // -- LAPACK auxiliary routine (version 3.2) --
44 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
45 // November 2006
46 //
47 // .. Scalar Arguments ..
48 // ..
49 //
50 // Purpose
51 // =======
52 //
53 // LSAME returns .TRUE. if CA is the same letter as CB regardless of
54 // case.
55 //
56 // Arguments
57 // =========
58 //
59 // CA (input) CHARACTER*1
60 // CB (input) CHARACTER*1
61 // CA and CB specify the single characters to be compared.
62 //
63 // This is the VIF version of the LAPACK routine.
64 //
65 if (ca != NO_TEXT && cb != NO_TEXT) {
66 return tolower (*ca) == tolower (*cb);
67 } else {
68 return FALSE;
69 }
70 }
71
72 logical_4 _lsamen (int_4 *n, char _p_ ca, char _p_ cb)
73 {
74 //
75 // -- LAPACK auxiliary routine (version 2.0) --
76 // Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
77 // Courant Institute, Argonne National Lab, and Rice University
78 // September 30, 1994
79 //
80 // Purpose
81 // =======
82 //
83 // LSAMEN tests if the first N letters of CA are the same as the
84 // first N letters of CB, regardless of case.
85 // LSAMEN returns .TRUE. if CA and CB are equivalent except for case
86 // and .FALSE. otherwise. LSAMEN also returns .FALSE. if LEN( CA )
87 // or LEN( CB ) is less than N.
88 //
89 // Arguments
90 // =========
91 //
92 // N (input) INTEGER
93 // The number of characters in CA and CB to be compared.
94 //
95 // CA (input) CHARACTER*(*)
96 // CB (input) CHARACTER*(*)
97 // CA and CB specify two character strings of length at least N.
98 // Only the first N characters of each string will be accessed.
99 //
100 // This is the VIF version of the LAPACK routine.
101 //
102 if (ca != NO_TEXT && cb != NO_TEXT) {
103 return strncasecmp (ca, cb, (size_t) *n) == 0;
104 } else {
105 return FALSE;
106 }
107 }
108
109 real_8 _dlamch (char _p_ cmach)
110 {
111 //
112 // -- LAPACK auxiliary routine (version 3.2) --
113 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
114 // November 2006
115 //
116 // .. Scalar Arguments ..
117 // ..
118 //
119 // Purpose
120 // =======
121 //
122 // DLAMCH determines double precision machine parameters.
123 //
124 // Arguments
125 // =========
126 //
127 // CMACH (input) CHARACTER*1
128 // Specifies the value to be returned by DLAMCH:
129 // = 'E' or 'e', DLAMCH := eps
130 // = 'S' or 's , DLAMCH := sfmin
131 // = 'B' or 'b', DLAMCH := base
132 // = 'P' or 'p', DLAMCH := eps*base
133 // = 'N' or 'n', DLAMCH := t
134 // = 'R' or 'r', DLAMCH := rnd
135 // = 'M' or 'm', DLAMCH := emin
136 // = 'U' or 'u', DLAMCH := rmin
137 // = 'L' or 'l', DLAMCH := emax
138 // = 'O' or 'o', DLAMCH := rmax
139 //
140 // where
141 //
142 // eps = relative machine precision
143 // sfmin = safe minimum, such that 1/sfmin does not overflow
144 // base = base of the machine
145 // prec = eps*base
146 // t = number of (base) digits in the mantissa
147 // rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
148 // emin = minimum exponent before (gradual) underflow
149 // rmin = underflow threshold - base**(emin-1)
150 // emax = largest exponent before overflow
151 // rmax = overflow threshold - (base**emax)*(1-eps)
152 //
153 int_4 beta, imin, imax, it, i__1;
154 logical_4 lrnd;
155 real_8 rmach = 0.0, small;
156 static logical_4 first = TRUE;
157 static real_8 emin, prec, emax, rmin, rmax, rnd, eps, base, sfmin, t;
158 if (first) {
159 _dlamc2 (&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
160 base = (real_8) beta;
161 t = (real_8) it;
162 if (lrnd) {
163 rnd = 1.0;
164 i__1 = 1 - it;
165 eps = _up_real_8 (base, i__1) / 2;
166 } else {
167 rnd = 0.0;
168 i__1 = 1 - it;
169 eps = _up_real_8 (base, i__1);
170 }
171 prec = eps * base;
172 emin = (real_8) imin;
173 emax = (real_8) imax;
174 sfmin = rmin;
175 small = 1.0 / rmax;
176 if (small >= sfmin) {
177 // Use SMALL plus a bit, to avoid the possibility of rounding
178 // causing overflow when computing 1/sfmin.
179 sfmin = small * (eps + 1.0);
180 }
181 }
182 //
183 if (_lsame (cmach, "E")) {
184 rmach = eps;
185 } else if (_lsame (cmach, "S")) {
186 rmach = sfmin;
187 } else if (_lsame (cmach, "B")) {
188 rmach = base;
189 } else if (_lsame (cmach, "P")) {
190 rmach = prec;
191 } else if (_lsame (cmach, "N")) {
192 rmach = t;
193 } else if (_lsame (cmach, "R")) {
194 rmach = rnd;
195 } else if (_lsame (cmach, "M")) {
196 rmach = emin;
197 } else if (_lsame (cmach, "U")) {
198 rmach = rmin;
199 } else if (_lsame (cmach, "L")) {
200 rmach = emax;
201 } else if (_lsame (cmach, "O")) {
202 rmach = rmax;
203 }
204 //
205 first = FALSE;
206 return rmach;
207 }
208
209 int_4 _dlamc1 (int_4 _p_ beta, int_4 _p_ t, logical_4 _p_ rnd, logical_4 _p_ ieee1)
210 {
211 //
212 //
213 // -- LAPACK auxiliary routine (version 3.2) --
214 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
215 // November 2006
216 //
217 // .. Scalar Arguments ..
218 // ..
219 //
220 // Purpose
221 // =======
222 //
223 // DLAMC1 determines the machine parameters given by BETA, T, RND, and
224 // IEEE1.
225 //
226 // Arguments
227 // =========
228 //
229 // BETA (output) INTEGER
230 // The base of the machine.
231 //
232 // T (output) INTEGER
233 // The number of ( BETA ) digits in the mantissa.
234 //
235 // RND (output) LOGICAL
236 // Specifies whether proper rounding ( RND = .TRUE. ) or
237 // chopping ( RND = .FALSE. ) occurs in addition. This may not
238 // be a reliable guide to the way in which the machine performs
239 // its arithmetic.
240 //
241 // IEEE1 (output) LOGICAL
242 // Specifies whether rounding appears to be done in the IEEE
243 // 'round to nearest' style.
244 //
245 // Further Details
246 // ===============
247 //
248 // The routine is based on the routine ENVRON by Malcolm and
249 // incorporates suggestions by Gentleman and Marovich. See
250 //
251 // Malcolm M. A. (1972) Algorithms to reveal properties of
252 // floating-point arithmetic. Comms. of the ACM, 15, 949-951.
253 //
254 // Gentleman W. M. and Marovich S. B. (1974) More on algorithms
255 // that reveal properties of floating point arithmetic units.
256 // Comms. of the ACM, 17, 276-277.
257 //
258 real_8 a, b, c__, f, t1, t2, d__1, d__2, one, qtr, savec;
259 static int_4 lbeta, lt;
260 static logical_4 first = TRUE, lieee1, lrnd;
261 if (first) {
262 one = 1.0;
263 // LBETA, LIEEE1, LT and LRND are the local values of BETA,
264 // IEEE1, T and RND.
265 //
266 // Throughout this routine we use the function DLAMC3 to ensure
267 // that relevant values are stored and not held in registers, or
268 // are not affected by optimizers.
269 //
270 // Compute a = 2.0**m with the smallest positive int_4 m such
271 // that
272 //
273 // fl( a + 1.0 ) = a.
274 //
275 a = 1.0;
276 c__ = 1.0;
277 // + WHILE( C.EQ.ONE )LOOP
278 while (c__ == one) {
279 a *= 2;
280 c__ = _dlamc3 (&a, &one);
281 d__1 = -a;
282 c__ = _dlamc3 (&c__, &d__1);
283 }
284 // + END WHILE
285 //
286 // Now compute b = 2.0**m with the smallest positive int_4 m
287 // such that
288 //
289 // fl( a + b ) .gt. a.
290 b = 1.0;
291 c__ = _dlamc3 (&a, &b);
292 // + WHILE( C.EQ.A )LOOP
293 while (c__ == a) {
294 b *= 2;
295 c__ = _dlamc3 (&a, &b);
296 }
297 // + END WHILE
298 //
299 // Now compute the base. a and c are neighbouring floating point
300 // numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
301 // their difference is beta. Adding 0.25 to c is to ensure that it
302 // is truncated to beta and not ( beta - 1 ).
303 qtr = one / 4;
304 savec = c__;
305 d__1 = -a;
306 c__ = _dlamc3 (&c__, &d__1);
307 lbeta = (int_4) (c__ + qtr);
308 // Now determine whether rounding or chopping occurs, by adding a
309 // bit less than beta/2 and a bit more than beta/2 to a.
310 b = (real_8) lbeta;
311 d__1 = b / 2;
312 d__2 = -b / 100;
313 f = _dlamc3 (&d__1, &d__2);
314 c__ = _dlamc3 (&f, &a);
315 if (c__ == a) {
316 lrnd = TRUE;
317 } else {
318 lrnd = FALSE;
319 }
320 d__1 = b / 2;
321 d__2 = b / 100;
322 f = _dlamc3 (&d__1, &d__2);
323 c__ = _dlamc3 (&f, &a);
324 if (lrnd && c__ == a) {
325 lrnd = FALSE;
326 }
327 // Try and decide whether rounding is done in the IEEE 'round to
328 // nearest' style. B/2 is half a unit in the last place of the two
329 // numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
330 // zero, and SAVEC is odd. Thus adding B/2 to A should not change
331 // A, but adding B/2 to SAVEC should change SAVEC.
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 // Now find the mantissa, t. It should be the int_4 part of
338 // log to the base beta of a, however it is safer to determine t
339 // by powering. So we find t as the smallest positive int_4 for
340 // which
341 //
342 // fl( beta**t + 1.0 ) = 1.0.
343 lt = 0;
344 a = 1.0;
345 c__ = 1.0;
346 // + WHILE( C.EQ.ONE )LOOP
347 while (c__ == one) {
348 ++lt;
349 a *= lbeta;
350 c__ = _dlamc3 (&a, &one);
351 d__1 = -a;
352 c__ = _dlamc3 (&c__, &d__1);
353 }
354 // + END WHILE
355 }
356 *beta = lbeta;
357 *t = lt;
358 *rnd = lrnd;
359 *ieee1 = lieee1;
360 first = FALSE;
361 return 0;
362 }
363
364 int_4 _dlamc2 (int_4 _p_ beta, int_4 _p_ t, logical_4 _p_ rnd, real_8 _p_ eps, int_4 _p_ emin,
365 real_8 _p_ rmin, int_4 _p_ emax, real_8 _p_ rmax)
366 {
367 //
368 // -- LAPACK auxiliary routine (version 3.2) --
369 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
370 // November 2006
371 //
372 // .. Scalar Arguments ..
373 // ..
374 //
375 // Purpose
376 // =======
377 //
378 // DLAMC2 determines the machine parameters specified in its argument
379 // list.
380 //
381 // Arguments
382 // =========
383 //
384 // BETA (output) INTEGER
385 // The base of the machine.
386 //
387 // T (output) INTEGER
388 // The number of ( BETA ) digits in the mantissa.
389 //
390 // RND (output) LOGICAL
391 // Specifies whether proper rounding ( RND = .TRUE. ) or
392 // chopping ( RND = .FALSE. ) occurs in addition. This may not
393 // be a reliable guide to the way in which the machine performs
394 // its arithmetic.
395 //
396 // EPS (output) DOUBLE PRECISION
397 // The smallest positive number such that
398 //
399 // fl( 1.0 - EPS ) .LT. 1.0,
400 //
401 // where fl denotes the computed value.
402 //
403 // EMIN (output) INTEGER
404 // The minimum exponent before (gradual) underflow occurs.
405 //
406 // RMIN (output) DOUBLE PRECISION
407 // The smallest normalized number for the machine, given by
408 // BASE**( EMIN - 1 ), where BASE is the floating point value
409 // of BETA.
410 //
411 // EMAX (output) INTEGER
412 // The maximum exponent before overflow occurs.
413 //
414 // RMAX (output) DOUBLE PRECISION
415 // The largest positive number for the machine, given by
416 // BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
417 // value of BETA.
418 //
419 // Further Details
420 // ===============
421 //
422 // The computation of EPS is based on a routine PARANOIA by
423 // W. Kahan of the University of California at Berkeley.
424 //
425 int_4 gnmin, gpmin, i__, i__1, ngnmin, ngpmin;
426 logical_4 ieee, lieee1, lrnd = TRUE;
427 real_8 a, b, c__, d__1, d__2, d__3, d__4, d__5;
428 real_8 half, one, two, rbase, sixth, small, third, zero;
429 static int_4 lbeta, lemin, lemax, lt;
430 static logical_4 first = TRUE, iwarn = FALSE;
431 static real_8 leps, lrmin, lrmax;
432 if (first) {
433 zero = 0.0;
434 one = 1.0;
435 two = 2.0;
436 // LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
437 // BETA, T, RND, EPS, EMIN and RMIN.
438 //
439 // Throughout this routine we use the function DLAMC3 to ensure
440 // that relevant values are stored and not held in registers, or
441 // are not affected by optimizers.
442 //
443 // DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
444 _dlamc1 (&lbeta, <, &lrnd, &lieee1);
445 // Start to find EPS.
446 b = (real_8) lbeta;
447 i__1 = -lt;
448 a = _up_real_8 (b, i__1);
449 leps = a;
450 // Try some tricks to see whether or not this is the correct EPS.
451 b = two / 3;
452 half = one / 2;
453 d__1 = -half;
454 sixth = _dlamc3 (&b, &d__1);
455 third = _dlamc3 (&sixth, &sixth);
456 d__1 = -half;
457 b = _dlamc3 (&third, &d__1);
458 b = _dlamc3 (&b, &sixth);
459 b = abs (b);
460 if (b < leps) {
461 b = leps;
462 }
463 leps = 1.0;
464 // + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
465 while (leps > b && b > zero) {
466 leps = b;
467 d__1 = half * leps;
468 // Computing 5th power
469 d__3 = two, d__4 = d__3, d__3 *= d__3;
470 // Computing 2nd power
471 d__5 = leps;
472 d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
473 c__ = _dlamc3 (&d__1, &d__2);
474 d__1 = -c__;
475 c__ = _dlamc3 (&half, &d__1);
476 b = _dlamc3 (&half, &c__);
477 d__1 = -b;
478 c__ = _dlamc3 (&half, &d__1);
479 b = _dlamc3 (&half, &c__);
480 }
481 // + END WHILE
482 if (a < leps) {
483 leps = a;
484 }
485 // Computation of EPS complete.
486 //
487 // Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
488 // Keep dividing A by BETA until (gradual) underflow occurs. This
489 // is detected when we cannot recover the previous A.
490 rbase = one / lbeta;
491 small = one;
492 for (i__ = 1; i__ <= 3; ++i__) {
493 d__1 = small * rbase;
494 small = _dlamc3 (&d__1, &zero);
495 }
496 a = _dlamc3 (&one, &small);
497 _dlamc4 (&ngpmin, &one, &lbeta);
498 d__1 = -one;
499 _dlamc4 (&ngnmin, &d__1, &lbeta);
500 _dlamc4 (&gpmin, &a, &lbeta);
501 d__1 = -a;
502 _dlamc4 (&gnmin, &d__1, &lbeta);
503 ieee = FALSE;
504 if (ngpmin == ngnmin && gpmin == gnmin) {
505 if (ngpmin == gpmin) {
506 lemin = ngpmin;
507 // Non twos-complement machines, no gradual underflow; e.g., VAX
508 } else if (gpmin - ngpmin == 3) {
509 lemin = ngpmin - 1 + lt;
510 ieee = TRUE;
511 // Non twos-complement machines, with gradual underflow; e.g., IEEE standard followers
512 } else {
513 lemin = _min (ngpmin, gpmin);
514 // A guess; no known machine
515 iwarn = TRUE;
516 }
517 } else if (ngpmin == gpmin && ngnmin == gnmin) {
518 if ((i__1 = ngpmin - ngnmin, abs (i__1)) == 1) {
519 lemin = _max (ngpmin, ngnmin);
520 // Twos-complement machines, no gradual underflow; e.g., CYBER 205
521 } else {
522 lemin = _min (ngpmin, ngnmin);
523 // A guess; no known machine
524 iwarn = TRUE;
525 }
526 } else if ((i__1 = ngpmin - ngnmin, abs (i__1)) == 1 && gpmin == gnmin) {
527 if (gpmin - _min (ngpmin, ngnmin) == 3) {
528 lemin = _max (ngpmin, ngnmin) - 1 + lt;
529 // Twos-complement machines with gradual underflow; no known machine
530 } else {
531 lemin = _min (ngpmin, ngnmin);
532 // A guess; no known machine
533 iwarn = TRUE;
534 }
535 } else {
536 // Computing MIN
537 i__1 = _min (ngpmin, ngnmin), i__1 = _min (i__1, gpmin);
538 lemin = _min (i__1, gnmin);
539 // A guess; no known machine
540 iwarn = TRUE;
541 }
542 first = FALSE;
543 // **
544 // Comment out this if block if EMIN is ok
545 if (iwarn) {
546 fprintf (stderr, "WARNING. The value EMIN may be incorrect: EMIN = %d",
547 lemin);
548 fprintf (stderr,
549 "If, after inspection, the value EMIN looks acceptable please");
550 fprintf (stderr,
551 "comment out the IF block as marked within the code of routine,");
552 fprintf (stderr, "DLAMC2, otherwise supply EMIN explicitly.");
553 first = TRUE;
554 }
555 // Assume IEEE arithmetic if we found denormalised numbers above,
556 // or if arithmetic seems to round in the IEEE style, determined
557 // in routine DLAMC1. A true IEEE machine should have both things
558 // true; however, faulty machines may have one or the other.
559 ieee = ieee || lieee1;
560 // Compute RMIN by successive division by BETA. We could compute
561 // RMIN as BASE**( EMIN - 1 ), but some machines underflow during
562 // this computation.
563 lrmin = 1.0;
564 i__1 = 1 - lemin;
565 for (i__ = 1; i__ <= i__1; ++i__) {
566 d__1 = lrmin * rbase;
567 lrmin = _dlamc3 (&d__1, &zero);
568 }
569 // Finally, call DLAMC5 to compute EMAX and RMAX.
570 _dlamc5 (&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
571 }
572 *beta = lbeta;
573 *t = lt;
574 *rnd = lrnd;
575 *eps = leps;
576 *emin = lemin;
577 *rmin = lrmin;
578 *emax = lemax;
579 *rmax = lrmax;
580 return 0;
581 }
582
583 real_8 _dlamc3 (real_8 _p_ a, real_8 _p_ b)
584 {
585 //
586 // -- LAPACK auxiliary routine (version 3.2) --
587 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
588 // November 2006
589 //
590 // .. Scalar Arguments ..
591 // ..
592 //
593 // Purpose
594 // =======
595 //
596 // DLAMC3 is intended to force A and B to be stored prior to doing
597 // the addition of A and B , for use in situations where optimizers
598 // might hold one of these in a register.
599 //
600 // Arguments
601 // =========
602 //
603 // A (input) DOUBLE PRECISION
604 // B (input) DOUBLE PRECISION
605 // The values A and B.
606 //
607 return (*a) + (*b);
608 }
609
610 int_4 _dlamc4 (int_4 _p_ emin, real_8 _p_ start, int_4 _p_ base)
611 {
612 //
613 // -- LAPACK auxiliary routine (version 3.2) --
614 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
615 // November 2006
616 //
617 // .. Scalar Arguments ..
618 // ..
619 //
620 // Purpose
621 // =======
622 //
623 // DLAMC4 is a service routine for DLAMC2.
624 //
625 // Arguments
626 // =========
627 //
628 // EMIN (output) INTEGER
629 // The minimum exponent before (gradual) underflow, computed by
630 // setting A = START and dividing by BASE until the previous A
631 // can not be recovered.
632 //
633 // START (input) DOUBLE PRECISION
634 // The starting point for determining EMIN.
635 //
636 // BASE (input) INTEGER
637 // The base of the machine.
638 //
639 int_4 i__, i__1;
640 real_8 a, b1, b2, c1, c2, d__1, d1, d2, one, zero, rbase;
641 a = *start;
642 one = 1.0;
643 rbase = one / *base;
644 zero = 0.0;
645 *emin = 1;
646 d__1 = a * rbase;
647 b1 = _dlamc3 (&d__1, &zero);
648 c1 = a;
649 c2 = a;
650 d1 = a;
651 d2 = a;
652 // + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
653 // $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
654 while (c1 == a && c2 == a && d1 == a && d2 == a) {
655 --(*emin);
656 a = b1;
657 d__1 = a / *base;
658 b1 = _dlamc3 (&d__1, &zero);
659 d__1 = b1 * *base;
660 c1 = _dlamc3 (&d__1, &zero);
661 d1 = zero;
662 i__1 = *base;
663 for (i__ = 1; i__ <= i__1; ++i__) {
664 d1 += b1;
665 }
666 d__1 = a * rbase;
667 b2 = _dlamc3 (&d__1, &zero);
668 d__1 = b2 / rbase;
669 c2 = _dlamc3 (&d__1, &zero);
670 d2 = zero;
671 i__1 = *base;
672 for (i__ = 1; i__ <= i__1; ++i__) {
673 d2 += b2;
674 }
675 }
676 // + END WHILE
677 return 0;
678 }
679
680 int_4 _dlamc5 (int_4 _p_ beta, int_4 _p_ p, int_4 _p_ emin, logical_4 _p_ ieee, int_4 _p_ emax, real_8 _p_ rmax)
681 {
682 //
683 // -- LAPACK auxiliary routine (version 3.2) --
684 // Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
685 // November 2006
686 //
687 // .. Scalar Arguments ..
688 // ..
689 //
690 // Purpose
691 // =======
692 //
693 // DLAMC5 attempts to compute RMAX, the largest machine floating-point
694 // number, without overflow. It assumes that EMAX + abs(EMIN) sum
695 // approximately to a power of 2. It will fail on machines where this
696 // assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
697 // EMAX = 28718). It will also fail if the value supplied for EMIN is
698 // too large (i.e. too close to zero), probably with overflow.
699 //
700 // Arguments
701 // =========
702 //
703 // BETA (input) INTEGER
704 // The base of floating-point arithmetic.
705 //
706 // P (input) INTEGER
707 // The number of base BETA digits in the mantissa of a
708 // floating-point value.
709 //
710 // EMIN (input) INTEGER
711 // The minimum exponent before (gradual) underflow.
712 //
713 // IEEE (input) LOGICAL
714 // A logical flag specifying whether or not the arithmetic
715 // system is thought to comply with the IEEE standard.
716 //
717 // EMAX (output) INTEGER
718 // The largest exponent before overflow
719 //
720 // RMAX (output) DOUBLE PRECISION
721 // The largest machine floating-point number.
722 //
723 // First compute LEXP and UEXP, two powers of 2 that bound
724 // abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
725 // approximately to the bound that is closest to abs(EMIN).
726 // (EMAX is the exponent of the required number RMAX).
727 //
728 int_4 exbits, expsum, i__, i__1, try, lexp, uexp, nbits;
729 real_8 d__1, oldy, recbas;
730 lexp = 1;
731 exbits = 1;
732 logical_4 cont = TRUE;
733 while (cont) {
734 try = lexp << 1;
735 if (try <= -(*emin)) {
736 lexp = try;
737 ++exbits;
738 } else {
739 cont = FALSE;
740 }
741 }
742 if (lexp == -(*emin)) {
743 uexp = lexp;
744 } else {
745 uexp = try;
746 ++exbits;
747 }
748 // Now -LEXP is less than or equal to EMIN, and -UEXP is greater
749 // than or equal to EMIN. EXBITS is the number of bits needed to
750 // store the exponent.
751 if (uexp + *emin > -lexp - *emin) {
752 expsum = lexp << 1;
753 } else {
754 expsum = uexp << 1;
755 }
756 // EXPSUM is the exponent range, approximately equal to
757 // EMAX - EMIN + 1 .
758 *emax = expsum + *emin - 1;
759 nbits = exbits + 1 + *p;
760 // NBITS is the total number of bits needed to store a
761 // floating-point number.
762 if (nbits % 2 == 1 && *beta == 2) {
763 // Either there are an odd number of bits used to store a
764 // floating-point number, which is unlikely, or some bits are
765 // not used in the representation of numbers, which is possible,
766 // (e.g. Cray machines) or the mantissa has an implicit bit,
767 // (e.g. IEEE machines, Dec Vax machines), which is perhaps the
768 // most likely. We have to assume the last alternative.
769 // If this is true, then we need to reduce EMAX by one because
770 // there must be some way of representing zero in an implicit-bit
771 // system. On machines like Cray, we are reducing EMAX by one
772 // unnecessarily.
773 --(*emax);
774 }
775 if (*ieee) {
776 // Assume we are on an IEEE machine which reserves one exponent
777 // for infinity and NaN.
778 --(*emax);
779 }
780 // Now create RMAX, the largest machine number, which should
781 // be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
782 //
783 // First compute 1.0 - BETA**(-P), being careful that the
784 // result is less than 1.0 .
785 recbas = 1.0 / *beta;
786 real_8 z__ = *beta - 1.0;
787 real_8 y = 0.0;
788 i__1 = *p;
789 for (i__ = 1; i__ <= i__1; ++i__) {
790 z__ *= recbas;
791 if (y < 1.0) {
792 oldy = y;
793 }
794 y = _dlamc3 (&y, &z__);
795 }
796 if (y >= 1.0) {
797 y = oldy;
798 }
799 // Now multiply by BETA**EMAX to get RMAX.
800 i__1 = *emax;
801 for (i__ = 1; i__ <= i__1; ++i__) {
802 d__1 = y * (*beta);
803 y = _dlamc3 (&d__1, &zero_arg);
804 }
805 *rmax = y;
806 return 0;
807 }
|
© 2002-2025 J.M. van der Veer (jmvdveer@xs4all.nl)
|
