14 typedef long long BLASLONG;
15 typedef unsigned long long BLASULONG;
17 typedef long BLASLONG;
18 typedef unsigned long BLASULONG;
22 typedef BLASLONG blasint;
24 #define blasabs(x) llabs(x)
26 #define blasabs(x) labs(x)
30 #define blasabs(x) abs(x)
33 typedef blasint integer;
35 typedef unsigned int uinteger;
36 typedef char *address;
37 typedef short int shortint;
39 typedef double doublereal;
40 typedef struct { real r, i; } complex;
41 typedef struct { doublereal r, i; } doublecomplex;
43 static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
44 static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
45 static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
46 static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
48 static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
49 static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
50 static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
51 static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
53 #define pCf(z) (*_pCf(z))
54 #define pCd(z) (*_pCd(z))
56 typedef short int shortlogical;
57 typedef char logical1;
58 typedef char integer1;
63 /* Extern is for use with -E */
74 /*external read, write*/
83 /*internal read, write*/
113 /*rewind, backspace, endfile*/
125 ftnint *inex; /*parameters in standard's order*/
151 union Multitype { /* for multiple entry points */
162 typedef union Multitype Multitype;
164 struct Vardesc { /* for Namelist */
170 typedef struct Vardesc Vardesc;
177 typedef struct Namelist Namelist;
179 #define abs(x) ((x) >= 0 ? (x) : -(x))
180 #define dabs(x) (fabs(x))
181 #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
182 #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
183 #define dmin(a,b) (f2cmin(a,b))
184 #define dmax(a,b) (f2cmax(a,b))
185 #define bit_test(a,b) ((a) >> (b) & 1)
186 #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
187 #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
189 #define abort_() { sig_die("Fortran abort routine called", 1); }
190 #define c_abs(z) (cabsf(Cf(z)))
191 #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
193 #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
194 #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/df(b)._Val[1]);}
196 #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
197 #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
199 #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
200 #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
201 #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
202 //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
203 #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
204 #define d_abs(x) (fabs(*(x)))
205 #define d_acos(x) (acos(*(x)))
206 #define d_asin(x) (asin(*(x)))
207 #define d_atan(x) (atan(*(x)))
208 #define d_atn2(x, y) (atan2(*(x),*(y)))
209 #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
210 #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
211 #define d_cos(x) (cos(*(x)))
212 #define d_cosh(x) (cosh(*(x)))
213 #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
214 #define d_exp(x) (exp(*(x)))
215 #define d_imag(z) (cimag(Cd(z)))
216 #define r_imag(z) (cimagf(Cf(z)))
217 #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
218 #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
219 #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
220 #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
221 #define d_log(x) (log(*(x)))
222 #define d_mod(x, y) (fmod(*(x), *(y)))
223 #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
224 #define d_nint(x) u_nint(*(x))
225 #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
226 #define d_sign(a,b) u_sign(*(a),*(b))
227 #define r_sign(a,b) u_sign(*(a),*(b))
228 #define d_sin(x) (sin(*(x)))
229 #define d_sinh(x) (sinh(*(x)))
230 #define d_sqrt(x) (sqrt(*(x)))
231 #define d_tan(x) (tan(*(x)))
232 #define d_tanh(x) (tanh(*(x)))
233 #define i_abs(x) abs(*(x))
234 #define i_dnnt(x) ((integer)u_nint(*(x)))
235 #define i_len(s, n) (n)
236 #define i_nint(x) ((integer)u_nint(*(x)))
237 #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
238 #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
239 #define pow_si(B,E) spow_ui(*(B),*(E))
240 #define pow_ri(B,E) spow_ui(*(B),*(E))
241 #define pow_di(B,E) dpow_ui(*(B),*(E))
242 #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
243 #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
244 #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
245 #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
246 #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
247 #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
248 #define sig_die(s, kill) { exit(1); }
249 #define s_stop(s, n) {exit(0);}
250 static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
251 #define z_abs(z) (cabs(Cd(z)))
252 #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
253 #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
254 #define myexit_() break;
255 #define mycycle() continue;
256 #define myceiling(w) {ceil(w)}
257 #define myhuge(w) {HUGE_VAL}
258 //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
259 #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
261 /* procedure parameter types for -A and -C++ */
263 #define F2C_proc_par_types 1
265 typedef logical (*L_fp)(...);
267 typedef logical (*L_fp)();
270 static float spow_ui(float x, integer n) {
271 float pow=1.0; unsigned long int u;
273 if(n < 0) n = -n, x = 1/x;
282 static double dpow_ui(double x, integer n) {
283 double pow=1.0; unsigned long int u;
285 if(n < 0) n = -n, x = 1/x;
295 static _Fcomplex cpow_ui(complex x, integer n) {
296 complex pow={1.0,0.0}; unsigned long int u;
298 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
300 if(u & 01) pow.r *= x.r, pow.i *= x.i;
301 if(u >>= 1) x.r *= x.r, x.i *= x.i;
305 _Fcomplex p={pow.r, pow.i};
309 static _Complex float cpow_ui(_Complex float x, integer n) {
310 _Complex float pow=1.0; unsigned long int u;
312 if(n < 0) n = -n, x = 1/x;
323 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
324 _Dcomplex pow={1.0,0.0}; unsigned long int u;
326 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
328 if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
329 if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
333 _Dcomplex p = {pow._Val[0], pow._Val[1]};
337 static _Complex double zpow_ui(_Complex double x, integer n) {
338 _Complex double pow=1.0; unsigned long int u;
340 if(n < 0) n = -n, x = 1/x;
350 static integer pow_ii(integer x, integer n) {
351 integer pow; unsigned long int u;
353 if (n == 0 || x == 1) pow = 1;
354 else if (x != -1) pow = x == 0 ? 1/x : 0;
357 if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
367 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
369 double m; integer i, mi;
370 for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
371 if (w[i-1]>m) mi=i ,m=w[i-1];
374 static integer smaxloc_(float *w, integer s, integer e, integer *n)
376 float m; integer i, mi;
377 for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
378 if (w[i-1]>m) mi=i ,m=w[i-1];
381 static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
382 integer n = *n_, incx = *incx_, incy = *incy_, i;
384 _Fcomplex zdotc = {0.0, 0.0};
385 if (incx == 1 && incy == 1) {
386 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
387 zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
388 zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
391 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
392 zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
393 zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
399 _Complex float zdotc = 0.0;
400 if (incx == 1 && incy == 1) {
401 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
402 zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
405 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
406 zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
412 static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
413 integer n = *n_, incx = *incx_, incy = *incy_, i;
415 _Dcomplex zdotc = {0.0, 0.0};
416 if (incx == 1 && incy == 1) {
417 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
418 zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
419 zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
422 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
423 zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
424 zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
430 _Complex double zdotc = 0.0;
431 if (incx == 1 && incy == 1) {
432 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
433 zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
436 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
437 zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
443 static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
444 integer n = *n_, incx = *incx_, incy = *incy_, i;
446 _Fcomplex zdotc = {0.0, 0.0};
447 if (incx == 1 && incy == 1) {
448 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
449 zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
450 zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
453 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
454 zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
455 zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
461 _Complex float zdotc = 0.0;
462 if (incx == 1 && incy == 1) {
463 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
464 zdotc += Cf(&x[i]) * Cf(&y[i]);
467 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
468 zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
474 static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
475 integer n = *n_, incx = *incx_, incy = *incy_, i;
477 _Dcomplex zdotc = {0.0, 0.0};
478 if (incx == 1 && incy == 1) {
479 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
480 zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
481 zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
484 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
485 zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
486 zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
492 _Complex double zdotc = 0.0;
493 if (incx == 1 && incy == 1) {
494 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
495 zdotc += Cd(&x[i]) * Cd(&y[i]);
498 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
499 zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
505 /* -- translated by f2c (version 20000121).
506 You must link the resulting object file with the libraries:
507 -lf2c -lm (in that order)
513 /* Table of constant values */
515 static integer c__1 = 1;
516 static doublereal c_b32 = 0.;
518 /* > \brief \b DLAMCHF77 deprecated */
520 /* =========== DOCUMENTATION =========== */
522 /* Online html documentation available at */
523 /* http://www.netlib.org/lapack/explore-html/ */
528 /* DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) */
530 /* CHARACTER CMACH */
533 /* > \par Purpose: */
538 /* > DLAMCHF77 determines double precision machine parameters. */
544 /* > \param[in] CMACH */
546 /* > Specifies the value to be returned by DLAMCH: */
547 /* > = 'E' or 'e', DLAMCH := eps */
548 /* > = 'S' or 's , DLAMCH := sfmin */
549 /* > = 'B' or 'b', DLAMCH := base */
550 /* > = 'P' or 'p', DLAMCH := eps*base */
551 /* > = 'N' or 'n', DLAMCH := t */
552 /* > = 'R' or 'r', DLAMCH := rnd */
553 /* > = 'M' or 'm', DLAMCH := emin */
554 /* > = 'U' or 'u', DLAMCH := rmin */
555 /* > = 'L' or 'l', DLAMCH := emax */
556 /* > = 'O' or 'o', DLAMCH := rmax */
558 /* > eps = relative machine precision */
559 /* > sfmin = safe minimum, such that 1/sfmin does not overflow */
560 /* > base = base of the machine */
561 /* > prec = eps*base */
562 /* > t = number of (base) digits in the mantissa */
563 /* > rnd = 1.0 when rounding occurs in addition, 0.0 otherwise */
564 /* > emin = minimum exponent before (gradual) underflow */
565 /* > rmin = underflow threshold - base**(emin-1) */
566 /* > emax = largest exponent before overflow */
567 /* > rmax = overflow threshold - (base**emax)*(1-eps) */
573 /* > \author Univ. of Tennessee */
574 /* > \author Univ. of California Berkeley */
575 /* > \author Univ. of Colorado Denver */
576 /* > \author NAG Ltd. */
578 /* > \date April 2012 */
580 /* > \ingroup auxOTHERauxiliary */
582 /* ===================================================================== */
583 doublereal dlamch_(char *cmach)
585 /* Initialized data */
587 static logical first = TRUE_;
589 /* System generated locals */
593 /* Local variables */
594 static doublereal base;
596 static doublereal emin, prec, emax;
599 static doublereal rmin, rmax, t;
601 extern logical lsame_(char *, char *);
603 static doublereal sfmin;
604 extern /* Subroutine */ int dlamc2_(integer *, integer *, logical *,
605 doublereal *, integer *, doublereal *, integer *, doublereal *);
607 static doublereal rnd, eps;
610 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
611 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
612 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
617 dlamc2_(&beta, &it, &lrnd, &eps, &imin, &rmin, &imax, &rmax);
618 base = (doublereal) beta;
623 eps = pow_di(&base, &i__1) / 2;
627 eps = pow_di(&base, &i__1);
630 emin = (doublereal) imin;
631 emax = (doublereal) imax;
634 if (small >= sfmin) {
636 /* Use SMALL plus a bit, to avoid the possibility of rounding */
637 /* causing overflow when computing 1/sfmin. */
639 sfmin = small * (eps + 1.);
643 if (lsame_(cmach, "E")) {
645 } else if (lsame_(cmach, "S")) {
647 } else if (lsame_(cmach, "B")) {
649 } else if (lsame_(cmach, "P")) {
651 } else if (lsame_(cmach, "N")) {
653 } else if (lsame_(cmach, "R")) {
655 } else if (lsame_(cmach, "M")) {
657 } else if (lsame_(cmach, "U")) {
659 } else if (lsame_(cmach, "L")) {
661 } else if (lsame_(cmach, "O")) {
674 /* *********************************************************************** */
676 /* > \brief \b DLAMC1 */
680 /* > DLAMC1 determines the machine parameters given by BETA, T, RND, and */
684 /* > \param[out] BETA */
686 /* > The base of the machine. */
689 /* > \param[out] T */
691 /* > The number of ( BETA ) digits in the mantissa. */
694 /* > \param[out] RND */
696 /* > Specifies whether proper rounding ( RND = .TRUE. ) or */
697 /* > chopping ( RND = .FALSE. ) occurs in addition. This may not */
698 /* > be a reliable guide to the way in which the machine performs */
699 /* > its arithmetic. */
702 /* > \param[out] IEEE1 */
704 /* > Specifies whether rounding appears to be done in the IEEE */
705 /* > 'round to nearest' style. */
707 /* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ.
708 of Colorado Denver and NAG Ltd.. */
709 /* > \date April 2012 */
710 /* > \ingroup auxOTHERauxiliary */
712 /* > \details \b Further \b Details */
715 /* > The routine is based on the routine ENVRON by Malcolm and */
716 /* > incorporates suggestions by Gentleman and Marovich. See */
718 /* > Malcolm M. A. (1972) Algorithms to reveal properties of */
719 /* > floating-point arithmetic. Comms. of the ACM, 15, 949-951. */
721 /* > Gentleman W. M. and Marovich S. B. (1974) More on algorithms */
722 /* > that reveal properties of floating point arithmetic units. */
723 /* > Comms. of the ACM, 17, 276-277. */
726 /* Subroutine */ int dlamc1_(integer *beta, integer *t, logical *rnd, logical
729 /* Initialized data */
731 static logical first = TRUE_;
733 /* System generated locals */
734 doublereal d__1, d__2;
736 /* Local variables */
738 doublereal a, b, c__, f;
739 static integer lbeta;
741 extern doublereal dlamc3_(doublereal *, doublereal *);
742 static logical lieee1;
748 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
749 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
752 /* ===================================================================== */
758 /* LBETA, LIEEE1, LT and LRND are the local values of BETA, */
759 /* IEEE1, T and RND. */
761 /* Throughout this routine we use the function DLAMC3 to ensure */
762 /* that relevant values are stored and not held in registers, or */
763 /* are not affected by optimizers. */
765 /* Compute a = 2.0**m with the smallest positive integer m such */
768 /* fl( a + 1.0 ) = a. */
773 /* + WHILE( C.EQ.ONE )LOOP */
777 c__ = dlamc3_(&a, &one);
779 c__ = dlamc3_(&c__, &d__1);
784 /* Now compute b = 2.0**m with the smallest positive integer m */
787 /* fl( a + b ) .gt. a. */
790 c__ = dlamc3_(&a, &b);
792 /* + WHILE( C.EQ.A )LOOP */
796 c__ = dlamc3_(&a, &b);
801 /* Now compute the base. a and c are neighbouring floating point */
802 /* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so */
803 /* their difference is beta. Adding 0.25 to c is to ensure that it */
804 /* is truncated to beta and not ( beta - 1 ). */
809 c__ = dlamc3_(&c__, &d__1);
810 lbeta = (integer) (c__ + qtr);
812 /* Now determine whether rounding or chopping occurs, by adding a */
813 /* bit less than beta/2 and a bit more than beta/2 to a. */
815 b = (doublereal) lbeta;
818 f = dlamc3_(&d__1, &d__2);
819 c__ = dlamc3_(&f, &a);
827 f = dlamc3_(&d__1, &d__2);
828 c__ = dlamc3_(&f, &a);
829 if (lrnd && c__ == a) {
833 /* Try and decide whether rounding is done in the IEEE 'round to */
834 /* nearest' style. B/2 is half a unit in the last place of the two */
835 /* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit */
836 /* zero, and SAVEC is odd. Thus adding B/2 to A should not change */
837 /* A, but adding B/2 to SAVEC should change SAVEC. */
840 t1 = dlamc3_(&d__1, &a);
842 t2 = dlamc3_(&d__1, &savec);
843 lieee1 = t1 == a && t2 > savec && lrnd;
845 /* Now find the mantissa, t. It should be the integer part of */
846 /* log to the base beta of a, however it is safer to determine t */
847 /* by powering. So we find t as the smallest positive integer for */
850 /* fl( beta**t + 1.0 ) = 1.0. */
856 /* + WHILE( C.EQ.ONE )LOOP */
861 c__ = dlamc3_(&a, &one);
863 c__ = dlamc3_(&c__, &d__1);
882 /* *********************************************************************** */
884 /* > \brief \b DLAMC2 */
888 /* > DLAMC2 determines the machine parameters specified in its argument */
891 /* > \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ.
892 of Colorado Denver and NAG Ltd.. */
893 /* > \date April 2012 */
894 /* > \ingroup auxOTHERauxiliary */
896 /* > \param[out] BETA */
898 /* > The base of the machine. */
901 /* > \param[out] T */
903 /* > The number of ( BETA ) digits in the mantissa. */
906 /* > \param[out] RND */
908 /* > Specifies whether proper rounding ( RND = .TRUE. ) or */
909 /* > chopping ( RND = .FALSE. ) occurs in addition. This may not */
910 /* > be a reliable guide to the way in which the machine performs */
911 /* > its arithmetic. */
914 /* > \param[out] EPS */
916 /* > The smallest positive number such that */
917 /* > fl( 1.0 - EPS ) .LT. 1.0, */
918 /* > where fl denotes the computed value. */
921 /* > \param[out] EMIN */
923 /* > The minimum exponent before (gradual) underflow occurs. */
926 /* > \param[out] RMIN */
928 /* > The smallest normalized number for the machine, given by */
929 /* > BASE**( EMIN - 1 ), where BASE is the floating point value */
933 /* > \param[out] EMAX */
935 /* > The maximum exponent before overflow occurs. */
938 /* > \param[out] RMAX */
940 /* > The largest positive number for the machine, given by */
941 /* > BASE**EMAX * ( 1 - EPS ), where BASE is the floating point */
942 /* > value of BETA. */
945 /* > \details \b Further \b Details */
948 /* > The computation of EPS is based on a routine PARANOIA by */
949 /* > W. Kahan of the University of California at Berkeley. */
951 /* Subroutine */ int dlamc2_(integer *beta, integer *t, logical *rnd,
952 doublereal *eps, integer *emin, doublereal *rmin, integer *emax,
955 /* Initialized data */
957 static logical first = TRUE_;
958 static logical iwarn = FALSE_;
961 static char fmt_9999[] = "(//\002 WARNING. The value EMIN may be incorre"
962 "ct:-\002,\002 EMIN = \002,i8,/\002 If, after inspection, the va"
963 "lue EMIN looks\002,\002 acceptable please comment out \002,/\002"
964 " the IF block as marked within the code of routine\002,\002 DLAM"
965 "C2,\002,/\002 otherwise supply EMIN explicitly.\002,/)";
967 /* System generated locals */
969 doublereal d__1, d__2, d__3, d__4, d__5;
971 /* Local variables */
975 static doublereal leps;
976 doublereal zero, a, b, c__;
978 static integer lbeta;
980 static integer lemin, lemax;
985 static doublereal lrmin, lrmax;
987 extern /* Subroutine */ int dlamc1_(integer *, integer *, logical *,
989 extern doublereal dlamc3_(doublereal *, doublereal *);
991 extern /* Subroutine */ int dlamc4_(integer *, doublereal *, integer *),
992 dlamc5_(integer *, integer *, integer *, logical *, integer *,
995 integer ngnmin, ngpmin;
998 /* Fortran I/O blocks */
999 static cilist io___58 = { 0, 6, 0, fmt_9999, 0 };
1003 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
1004 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
1007 /* ===================================================================== */
1015 /* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of */
1016 /* BETA, T, RND, EPS, EMIN and RMIN. */
1018 /* Throughout this routine we use the function DLAMC3 to ensure */
1019 /* that relevant values are stored and not held in registers, or */
1020 /* are not affected by optimizers. */
1022 /* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. */
1024 dlamc1_(&lbeta, <, &lrnd, &lieee1);
1026 /* Start to find EPS. */
1028 b = (doublereal) lbeta;
1030 a = pow_di(&b, &i__1);
1033 /* Try some tricks to see whether or not this is the correct EPS. */
1038 sixth = dlamc3_(&b, &d__1);
1039 third = dlamc3_(&sixth, &sixth);
1041 b = dlamc3_(&third, &d__1);
1042 b = dlamc3_(&b, &sixth);
1050 /* + WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP */
1052 if (leps > b && b > zero) {
1055 /* Computing 5th power */
1056 d__3 = two, d__4 = d__3, d__3 *= d__3;
1057 /* Computing 2nd power */
1059 d__2 = d__4 * (d__3 * d__3) * (d__5 * d__5);
1060 c__ = dlamc3_(&d__1, &d__2);
1062 c__ = dlamc3_(&half, &d__1);
1063 b = dlamc3_(&half, &c__);
1065 c__ = dlamc3_(&half, &d__1);
1066 b = dlamc3_(&half, &c__);
1075 /* Computation of EPS complete. */
1077 /* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). */
1078 /* Keep dividing A by BETA until (gradual) underflow occurs. This */
1079 /* is detected when we cannot recover the previous A. */
1081 rbase = one / lbeta;
1083 for (i__ = 1; i__ <= 3; ++i__) {
1084 d__1 = small * rbase;
1085 small = dlamc3_(&d__1, &zero);
1088 a = dlamc3_(&one, &small);
1089 dlamc4_(&ngpmin, &one, &lbeta);
1091 dlamc4_(&ngnmin, &d__1, &lbeta);
1092 dlamc4_(&gpmin, &a, &lbeta);
1094 dlamc4_(&gnmin, &d__1, &lbeta);
1097 if (ngpmin == ngnmin && gpmin == gnmin) {
1098 if (ngpmin == gpmin) {
1100 /* ( Non twos-complement machines, no gradual underflow; */
1102 } else if (gpmin - ngpmin == 3) {
1103 lemin = ngpmin - 1 + lt;
1105 /* ( Non twos-complement machines, with gradual underflow; */
1106 /* e.g., IEEE standard followers ) */
1108 lemin = f2cmin(ngpmin,gpmin);
1109 /* ( A guess; no known machine ) */
1113 } else if (ngpmin == gpmin && ngnmin == gnmin) {
1114 if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1) {
1115 lemin = f2cmax(ngpmin,ngnmin);
1116 /* ( Twos-complement machines, no gradual underflow; */
1117 /* e.g., CYBER 205 ) */
1119 lemin = f2cmin(ngpmin,ngnmin);
1120 /* ( A guess; no known machine ) */
1124 } else if ((i__1 = ngpmin - ngnmin, abs(i__1)) == 1 && gpmin == gnmin)
1126 if (gpmin - f2cmin(ngpmin,ngnmin) == 3) {
1127 lemin = f2cmax(ngpmin,ngnmin) - 1 + lt;
1128 /* ( Twos-complement machines with gradual underflow; */
1129 /* no known machine ) */
1131 lemin = f2cmin(ngpmin,ngnmin);
1132 /* ( A guess; no known machine ) */
1138 i__1 = f2cmin(ngpmin,ngnmin), i__1 = f2cmin(i__1,gpmin);
1139 lemin = f2cmin(i__1,gnmin);
1140 /* ( A guess; no known machine ) */
1145 /* Comment out this if block if EMIN is ok */
1150 do_fio(&c__1, (char *)&lemin, (ftnlen)sizeof(integer));
1156 /* Assume IEEE arithmetic if we found denormalised numbers above, */
1157 /* or if arithmetic seems to round in the IEEE style, determined */
1158 /* in routine DLAMC1. A true IEEE machine should have both things */
1159 /* true; however, faulty machines may have one or the other. */
1161 ieee = ieee || lieee1;
1163 /* Compute RMIN by successive division by BETA. We could compute */
1164 /* RMIN as BASE**( EMIN - 1 ), but some machines underflow during */
1165 /* this computation. */
1169 for (i__ = 1; i__ <= i__1; ++i__) {
1170 d__1 = lrmin * rbase;
1171 lrmin = dlamc3_(&d__1, &zero);
1175 /* Finally, call DLAMC5 to compute EMAX and RMAX. */
1177 dlamc5_(&lbeta, <, &lemin, &ieee, &lemax, &lrmax);
1197 /* *********************************************************************** */
1199 /* > \brief \b DLAMC3 */
1203 /* > DLAMC3 is intended to force A and B to be stored prior to doing */
1204 /* > the addition of A and B , for use in situations where optimizers */
1205 /* > might hold one of these in a register. */
1206 /* > \endverbatim */
1208 /* > \param[in] A */
1210 /* > \param[in] B */
1212 /* > The values A and B. */
1213 /* > \endverbatim */
1214 doublereal dlamc3_(doublereal *a, doublereal *b)
1216 /* System generated locals */
1220 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
1221 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
1224 /* ===================================================================== */
1236 /* *********************************************************************** */
1238 /* > \brief \b DLAMC4 */
1242 /* > DLAMC4 is a service routine for DLAMC2. */
1243 /* > \endverbatim */
1245 /* > \param[out] EMIN */
1247 /* > The minimum exponent before (gradual) underflow, computed by */
1248 /* > setting A = START and dividing by BASE until the previous A */
1249 /* > can not be recovered. */
1250 /* > \endverbatim */
1252 /* > \param[in] START */
1254 /* > The starting point for determining EMIN. */
1255 /* > \endverbatim */
1257 /* > \param[in] BASE */
1259 /* > The base of the machine. */
1260 /* > \endverbatim */
1262 /* Subroutine */ int dlamc4_(integer *emin, doublereal *start, integer *base)
1264 /* System generated locals */
1268 /* Local variables */
1271 doublereal rbase, b1, b2, c1, c2, d1, d2;
1272 extern doublereal dlamc3_(doublereal *, doublereal *);
1276 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
1277 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
1280 /* ===================================================================== */
1285 rbase = one / *base;
1289 b1 = dlamc3_(&d__1, &zero);
1294 /* + WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. */
1295 /* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP */
1297 if (c1 == a && c2 == a && d1 == a && d2 == a) {
1301 b1 = dlamc3_(&d__1, &zero);
1303 c1 = dlamc3_(&d__1, &zero);
1306 for (i__ = 1; i__ <= i__1; ++i__) {
1311 b2 = dlamc3_(&d__1, &zero);
1313 c2 = dlamc3_(&d__1, &zero);
1316 for (i__ = 1; i__ <= i__1; ++i__) {
1331 /* *********************************************************************** */
1333 /* > \brief \b DLAMC5 */
1337 /* > DLAMC5 attempts to compute RMAX, the largest machine floating-point */
1338 /* > number, without overflow. It assumes that EMAX + abs(EMIN) sum */
1339 /* > approximately to a power of 2. It will fail on machines where this */
1340 /* > assumption does not hold, for example, the Cyber 205 (EMIN = -28625, */
1341 /* > EMAX = 28718). It will also fail if the value supplied for EMIN is */
1342 /* > too large (i.e. too close to zero), probably with overflow. */
1343 /* > \endverbatim */
1345 /* > \param[in] BETA */
1347 /* > The base of floating-point arithmetic. */
1348 /* > \endverbatim */
1350 /* > \param[in] P */
1352 /* > The number of base BETA digits in the mantissa of a */
1353 /* > floating-point value. */
1354 /* > \endverbatim */
1356 /* > \param[in] EMIN */
1358 /* > The minimum exponent before (gradual) underflow. */
1359 /* > \endverbatim */
1361 /* > \param[in] IEEE */
1363 /* > A logical flag specifying whether or not the arithmetic */
1364 /* > system is thought to comply with the IEEE standard. */
1365 /* > \endverbatim */
1367 /* > \param[out] EMAX */
1369 /* > The largest exponent before overflow */
1370 /* > \endverbatim */
1372 /* > \param[out] RMAX */
1374 /* > The largest machine floating-point number. */
1375 /* > \endverbatim */
1377 /* Subroutine */ int dlamc5_(integer *beta, integer *p, integer *emin,
1378 logical *ieee, integer *emax, doublereal *rmax)
1380 /* System generated locals */
1384 /* Local variables */
1390 extern doublereal dlamc3_(doublereal *, doublereal *);
1392 integer exbits, expsum, try__;
1395 /* -- LAPACK auxiliary routine (version 3.7.0) -- */
1396 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
1399 /* ===================================================================== */
1402 /* First compute LEXP and UEXP, two powers of 2 that bound */
1403 /* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum */
1404 /* approximately to the bound that is closest to abs(EMIN). */
1405 /* (EMAX is the exponent of the required number RMAX). */
1411 if (try__ <= -(*emin)) {
1416 if (lexp == -(*emin)) {
1423 /* Now -LEXP is less than or equal to EMIN, and -UEXP is greater */
1424 /* than or equal to EMIN. EXBITS is the number of bits needed to */
1425 /* store the exponent. */
1427 if (uexp + *emin > -lexp - *emin) {
1433 /* EXPSUM is the exponent range, approximately equal to */
1434 /* EMAX - EMIN + 1 . */
1436 *emax = expsum + *emin - 1;
1437 nbits = exbits + 1 + *p;
1439 /* NBITS is the total number of bits needed to store a */
1440 /* floating-point number. */
1442 if (nbits % 2 == 1 && *beta == 2) {
1444 /* Either there are an odd number of bits used to store a */
1445 /* floating-point number, which is unlikely, or some bits are */
1446 /* not used in the representation of numbers, which is possible, */
1447 /* (e.g. Cray machines) or the mantissa has an implicit bit, */
1448 /* (e.g. IEEE machines, Dec Vax machines), which is perhaps the */
1449 /* most likely. We have to assume the last alternative. */
1450 /* If this is true, then we need to reduce EMAX by one because */
1451 /* there must be some way of representing zero in an implicit-bit */
1452 /* system. On machines like Cray, we are reducing EMAX by one */
1453 /* unnecessarily. */
1460 /* Assume we are on an IEEE machine which reserves one exponent */
1461 /* for infinity and NaN. */
1466 /* Now create RMAX, the largest machine number, which should */
1467 /* be equal to (1.0 - BETA**(-P)) * BETA**EMAX . */
1469 /* First compute 1.0 - BETA**(-P), being careful that the */
1470 /* result is less than 1.0 . */
1472 recbas = 1. / *beta;
1476 for (i__ = 1; i__ <= i__1; ++i__) {
1481 y = dlamc3_(&y, &z__);
1488 /* Now multiply by BETA**EMAX to get RMAX. */
1491 for (i__ = 1; i__ <= i__1; ++i__) {
1493 y = dlamc3_(&d__1, &c_b32);