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