1: #region Translated by Jose Antonio De Santiago-Castillo.
2:
3: //Translated by Jose Antonio De Santiago-Castillo.
4: //E-mail:JAntonioDeSantiago@gmail.com
5: //Web: www.DotNumerics.com
6: //
7: //Fortran to C# Translation.
8: //Translated by:
9: //F2CSharp Version 0.71 (November 10, 2009)
10: //Code Optimizations: None
11: //
12: #endregion
13:
14: using System;
15: using DotNumerics.FortranLibrary;
16:
17: namespace DotNumerics.CSLapack
18: {
19:
20: #region The Class: DLAMCH
21:
22: /// <summary>
23: /// -- LAPACK auxiliary routine (version 3.1) --
24: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
25: /// November 2006
26: /// Purpose
27: /// =======
28: ///
29: /// DLAMCH determines double precision machine parameters.
30: ///
31: ///</summary>
32: public class DLAMCH
33: {
34:
35:
36: #region Dependencies
37:
38: LSAME _lsame; DLAMC2 _dlamc2;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; bool FIRST = false; bool LRND = false; int BETA = 0; int IMAX = 0;
46: int IMIN = 0;int IT = 0; double BASE = 0; double EMAX = 0; double EMIN = 0; double EPS = 0; double PREC = 0;
47: double RMACH = 0;double RMAX = 0; double RMIN = 0; double RND = 0; double SFMIN = 0; double SMALL = 0; double T = 0;
48:
49: #endregion
50:
51: public DLAMCH(LSAME lsame, DLAMC2 dlamc2)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._lsame = lsame; this._dlamc2 = dlamc2;
58:
59: #endregion
60:
61:
62: #region Data Inicializacion
63:
64: //FIRST/.TRUE.
65: FIRST = true;
66:
67: #endregion
68:
69: }
70:
71: public DLAMCH()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: LSAME lsame = new LSAME();
78: DLAMC3 dlamc3 = new DLAMC3();
79: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
80: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
81: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
82: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
83:
84: #endregion
85:
86:
87: #region Set Dependencies
88:
89: this._lsame = lsame; this._dlamc2 = dlamc2;
90:
91: #endregion
92:
93:
94: #region Data Inicializacion
95:
96: //FIRST/.TRUE.
97: FIRST = true;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DLAMCH determines double precision machine parameters.
107: ///
108: ///</summary>
109: /// <param name="CMACH">
110: /// (input) CHARACTER*1
111: /// Specifies the value to be returned by DLAMCH:
112: /// = 'E' or 'e', DLAMCH := eps
113: /// = 'S' or 's , DLAMCH := sfmin
114: /// = 'B' or 'b', DLAMCH := base
115: /// = 'P' or 'p', DLAMCH := eps*base
116: /// = 'N' or 'n', DLAMCH := t
117: /// = 'R' or 'r', DLAMCH := rnd
118: /// = 'M' or 'm', DLAMCH := emin
119: /// = 'U' or 'u', DLAMCH := rmin
120: /// = 'L' or 'l', DLAMCH := emax
121: /// = 'O' or 'o', DLAMCH := rmax
122: ///
123: /// where
124: ///
125: /// eps = relative machine precision
126: /// sfmin = safe minimum, such that 1/sfmin does not overflow
127: /// base = base of the machine
128: /// prec = eps*base
129: /// t = number of (base) digits in the mantissa
130: /// rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
131: /// emin = minimum exponent before (gradual) underflow
132: /// rmin = underflow threshold - base**(emin-1)
133: /// emax = largest exponent before overflow
134: /// rmax = overflow threshold - (base**emax)*(1-eps)
135: ///</param>
136: public double Run(string CMACH)
137: {
138: double dlamch = 0;
139:
140: #region Strings
141:
142: CMACH = CMACH.Substring(0, 1);
143:
144: #endregion
145:
146:
147: #region Prolog
148:
149: // *
150: // * -- LAPACK auxiliary routine (version 3.1) --
151: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
152: // * November 2006
153: // *
154: // * .. Scalar Arguments ..
155: // * ..
156: // *
157: // * Purpose
158: // * =======
159: // *
160: // * DLAMCH determines double precision machine parameters.
161: // *
162: // * Arguments
163: // * =========
164: // *
165: // * CMACH (input) CHARACTER*1
166: // * Specifies the value to be returned by DLAMCH:
167: // * = 'E' or 'e', DLAMCH := eps
168: // * = 'S' or 's , DLAMCH := sfmin
169: // * = 'B' or 'b', DLAMCH := base
170: // * = 'P' or 'p', DLAMCH := eps*base
171: // * = 'N' or 'n', DLAMCH := t
172: // * = 'R' or 'r', DLAMCH := rnd
173: // * = 'M' or 'm', DLAMCH := emin
174: // * = 'U' or 'u', DLAMCH := rmin
175: // * = 'L' or 'l', DLAMCH := emax
176: // * = 'O' or 'o', DLAMCH := rmax
177: // *
178: // * where
179: // *
180: // * eps = relative machine precision
181: // * sfmin = safe minimum, such that 1/sfmin does not overflow
182: // * base = base of the machine
183: // * prec = eps*base
184: // * t = number of (base) digits in the mantissa
185: // * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
186: // * emin = minimum exponent before (gradual) underflow
187: // * rmin = underflow threshold - base**(emin-1)
188: // * emax = largest exponent before overflow
189: // * rmax = overflow threshold - (base**emax)*(1-eps)
190: // *
191: // * =====================================================================
192: // *
193: // * .. Parameters ..
194: // * ..
195: // * .. Local Scalars ..
196: // * ..
197: // * .. External Functions ..
198: // * ..
199: // * .. External Subroutines ..
200: // * ..
201: // * .. Save statement ..
202: // * ..
203: // * .. Data statements ..
204: // * ..
205: // * .. Executable Statements ..
206: // *
207:
208: #endregion
209:
210:
211: #region Body
212:
213: if (FIRST)
214: {
215: this._dlamc2.Run(ref BETA, ref IT, ref LRND, ref EPS, ref IMIN, ref RMIN
216: , ref IMAX, ref RMAX);
217: BASE = BETA;
218: T = IT;
219: if (LRND)
220: {
221: RND = ONE;
222: EPS = (Math.Pow(BASE,1 - IT)) / 2;
223: }
224: else
225: {
226: RND = ZERO;
227: EPS = Math.Pow(BASE,1 - IT);
228: }
229: PREC = EPS * BASE;
230: EMIN = IMIN;
231: EMAX = IMAX;
232: SFMIN = RMIN;
233: SMALL = ONE / RMAX;
234: if (SMALL >= SFMIN)
235: {
236: // *
237: // * Use SMALL plus a bit, to avoid the possibility of rounding
238: // * causing overflow when computing 1/sfmin.
239: // *
240: SFMIN = SMALL * (ONE + EPS);
241: }
242: }
243: // *
244: if (this._lsame.Run(CMACH, "E"))
245: {
246: RMACH = EPS;
247: }
248: else
249: {
250: if (this._lsame.Run(CMACH, "S"))
251: {
252: RMACH = SFMIN;
253: }
254: else
255: {
256: if (this._lsame.Run(CMACH, "B"))
257: {
258: RMACH = BASE;
259: }
260: else
261: {
262: if (this._lsame.Run(CMACH, "P"))
263: {
264: RMACH = PREC;
265: }
266: else
267: {
268: if (this._lsame.Run(CMACH, "N"))
269: {
270: RMACH = T;
271: }
272: else
273: {
274: if (this._lsame.Run(CMACH, "R"))
275: {
276: RMACH = RND;
277: }
278: else
279: {
280: if (this._lsame.Run(CMACH, "M"))
281: {
282: RMACH = EMIN;
283: }
284: else
285: {
286: if (this._lsame.Run(CMACH, "U"))
287: {
288: RMACH = RMIN;
289: }
290: else
291: {
292: if (this._lsame.Run(CMACH, "L"))
293: {
294: RMACH = EMAX;
295: }
296: else
297: {
298: if (this._lsame.Run(CMACH, "O"))
299: {
300: RMACH = RMAX;
301: }
302: }
303: }
304: }
305: }
306: }
307: }
308: }
309: }
310: }
311: // *
312: dlamch = RMACH;
313: FIRST = false;
314: return dlamch;
315: // *
316: // * End of DLAMCH
317: // *
318:
319: #endregion
320:
321: }
322: }
323:
324: #endregion
325:
326:
327: #region The Class: DLAMC1
328:
329: // *
330: // ************************************************************************
331: // *
332: /// <summary>
333: /// -- LAPACK auxiliary routine (version 3.1) --
334: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
335: /// November 2006
336: /// Purpose
337: /// =======
338: ///
339: /// DLAMC1 determines the machine parameters given by BETA, T, RND, and
340: /// IEEE1.
341: ///
342: ///</summary>
343: public class DLAMC1
344: {
345:
346:
347: #region Dependencies
348:
349: DLAMC3 _dlamc3;
350:
351: #endregion
352:
353:
354: #region Fields
355:
356: bool FIRST = false; bool LIEEE1 = false; bool LRND = false; int LBETA = 0; int LT = 0; double A = 0; double B = 0;
357: double C = 0;double F = 0; double ONE = 0; double QTR = 0; double SAVEC = 0; double T1 = 0; double T2 = 0;
358:
359: #endregion
360:
361: public DLAMC1(DLAMC3 dlamc3)
362: {
363:
364:
365: #region Set Dependencies
366:
367: this._dlamc3 = dlamc3;
368:
369: #endregion
370:
371:
372: #region Data Inicializacion
373:
374: //FIRST/.TRUE.
375: FIRST = true;
376:
377: #endregion
378:
379: }
380:
381: public DLAMC1()
382: {
383:
384:
385: #region Dependencies (Initialization)
386:
387: DLAMC3 dlamc3 = new DLAMC3();
388:
389: #endregion
390:
391:
392: #region Set Dependencies
393:
394: this._dlamc3 = dlamc3;
395:
396: #endregion
397:
398:
399: #region Data Inicializacion
400:
401: //FIRST/.TRUE.
402: FIRST = true;
403:
404: #endregion
405:
406: }
407: /// <summary>
408: /// Purpose
409: /// =======
410: ///
411: /// DLAMC1 determines the machine parameters given by BETA, T, RND, and
412: /// IEEE1.
413: ///
414: ///</summary>
415: /// <param name="BETA">
416: /// (output) INTEGER
417: /// The base of the machine.
418: ///</param>
419: /// <param name="T">
420: /// (output) INTEGER
421: /// The number of ( BETA ) digits in the mantissa.
422: ///</param>
423: /// <param name="RND">
424: /// (output) LOGICAL
425: /// Specifies whether proper rounding ( RND = .TRUE. ) or
426: /// chopping ( RND = .FALSE. ) occurs in addition. This may not
427: /// be a reliable guide to the way in which the machine performs
428: /// its arithmetic.
429: ///</param>
430: /// <param name="IEEE1">
431: /// (output) LOGICAL
432: /// Specifies whether rounding appears to be done in the IEEE
433: /// 'round to nearest' style.
434: ///</param>
435: public void Run(ref int BETA, ref int T, ref bool RND, ref bool IEEE1)
436: {
437:
438: #region Prolog
439:
440: // *
441: // * -- LAPACK auxiliary routine (version 3.1) --
442: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
443: // * November 2006
444: // *
445: // * .. Scalar Arguments ..
446: // * ..
447: // *
448: // * Purpose
449: // * =======
450: // *
451: // * DLAMC1 determines the machine parameters given by BETA, T, RND, and
452: // * IEEE1.
453: // *
454: // * Arguments
455: // * =========
456: // *
457: // * BETA (output) INTEGER
458: // * The base of the machine.
459: // *
460: // * T (output) INTEGER
461: // * The number of ( BETA ) digits in the mantissa.
462: // *
463: // * RND (output) LOGICAL
464: // * Specifies whether proper rounding ( RND = .TRUE. ) or
465: // * chopping ( RND = .FALSE. ) occurs in addition. This may not
466: // * be a reliable guide to the way in which the machine performs
467: // * its arithmetic.
468: // *
469: // * IEEE1 (output) LOGICAL
470: // * Specifies whether rounding appears to be done in the IEEE
471: // * 'round to nearest' style.
472: // *
473: // * Further Details
474: // * ===============
475: // *
476: // * The routine is based on the routine ENVRON by Malcolm and
477: // * incorporates suggestions by Gentleman and Marovich. See
478: // *
479: // * Malcolm M. A. (1972) Algorithms to reveal properties of
480: // * floating-point arithmetic. Comms. of the ACM, 15, 949-951.
481: // *
482: // * Gentleman W. M. and Marovich S. B. (1974) More on algorithms
483: // * that reveal properties of floating point arithmetic units.
484: // * Comms. of the ACM, 17, 276-277.
485: // *
486: // * =====================================================================
487: // *
488: // * .. Local Scalars ..
489: // * ..
490: // * .. External Functions ..
491: // * ..
492: // * .. Save statement ..
493: // * ..
494: // * .. Data statements ..
495: // * ..
496: // * .. Executable Statements ..
497: // *
498:
499: #endregion
500:
501:
502: #region Body
503:
504: if (FIRST)
505: {
506: ONE = 1;
507: // *
508: // * LBETA, LIEEE1, LT and LRND are the local values of BETA,
509: // * IEEE1, T and RND.
510: // *
511: // * Throughout this routine we use the function DLAMC3 to ensure
512: // * that relevant values are stored and not held in registers, or
513: // * are not affected by optimizers.
514: // *
515: // * Compute a = 2.0**m with the smallest positive integer m such
516: // * that
517: // *
518: // * fl( a + 1.0 ) = a.
519: // *
520: A = 1;
521: C = 1;
522: // *
523: // *+ WHILE( C.EQ.ONE )LOOP
524: LABEL10:;
525: if (C == ONE)
526: {
527: A = 2 * A;
528: C = this._dlamc3.Run(A, ONE);
529: C = this._dlamc3.Run(C, - A);
530: goto LABEL10;
531: }
532: // *+ END WHILE
533: // *
534: // * Now compute b = 2.0**m with the smallest positive integer m
535: // * such that
536: // *
537: // * fl( a + b ) .gt. a.
538: // *
539: B = 1;
540: C = this._dlamc3.Run(A, B);
541: // *
542: // *+ WHILE( C.EQ.A )LOOP
543: LABEL20:;
544: if (C == A)
545: {
546: B = 2 * B;
547: C = this._dlamc3.Run(A, B);
548: goto LABEL20;
549: }
550: // *+ END WHILE
551: // *
552: // * Now compute the base. a and c are neighbouring floating point
553: // * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
554: // * their difference is beta. Adding 0.25 to c is to ensure that it
555: // * is truncated to beta and not ( beta - 1 ).
556: // *
557: QTR = ONE / 4;
558: C = this._dlamc3.Run(C, - A);
559: LBETA = Convert.ToInt32(C + QTR);
560: // *
561: // * Now determine whether rounding or chopping occurs, by adding a
562: // * bit less than beta/2 and a bit more than beta/2 to a.
563: // *
564: B = LBETA;
565: F = this._dlamc3.Run(B / 2, - B / 100);
566: C = this._dlamc3.Run(F, A);
567: if (C == A)
568: {
569: LRND = true;
570: }
571: else
572: {
573: LRND = false;
574: }
575: F = this._dlamc3.Run(B / 2, B / 100);
576: C = this._dlamc3.Run(F, A);
577: if ((LRND) && (C == A)) LRND = false;
578: // *
579: // * Try and decide whether rounding is done in the IEEE 'round to
580: // * nearest' style. B/2 is half a unit in the last place of the two
581: // * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
582: // * zero, and SAVEC is odd. Thus adding B/2 to A should not change
583: // * A, but adding B/2 to SAVEC should change SAVEC.
584: // *
585: T1 = this._dlamc3.Run(B / 2, A);
586: T2 = this._dlamc3.Run(B / 2, SAVEC);
587: LIEEE1 = (T1 == A) && (T2 > SAVEC) && LRND;
588: // *
589: // * Now find the mantissa, t. It should be the integer part of
590: // * log to the base beta of a, however it is safer to determine t
591: // * by powering. So we find t as the smallest positive integer for
592: // * which
593: // *
594: // * fl( beta**t + 1.0 ) = 1.0.
595: // *
596: LT = 0;
597: A = 1;
598: C = 1;
599: // *
600: // *+ WHILE( C.EQ.ONE )LOOP
601: LABEL30:;
602: if (C == ONE)
603: {
604: LT = LT + 1;
605: A = A * LBETA;
606: C = this._dlamc3.Run(A, ONE);
607: C = this._dlamc3.Run(C, - A);
608: goto LABEL30;
609: }
610: // *+ END WHILE
611: // *
612: }
613: // *
614: BETA = LBETA;
615: T = LT;
616: RND = LRND;
617: IEEE1 = LIEEE1;
618: FIRST = false;
619: return;
620: // *
621: // * End of DLAMC1
622: // *
623:
624: #endregion
625:
626: }
627: }
628:
629: #endregion
630:
631:
632: #region The Class: DLAMC2
633:
634: // *
635: // ************************************************************************
636: // *
637: /// <summary>
638: /// -- LAPACK auxiliary routine (version 3.1) --
639: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
640: /// November 2006
641: /// Purpose
642: /// =======
643: ///
644: /// DLAMC2 determines the machine parameters specified in its argument
645: /// list.
646: ///
647: ///</summary>
648: public class DLAMC2
649: {
650:
651:
652: #region Dependencies
653:
654: DLAMC3 _dlamc3; DLAMC1 _dlamc1; DLAMC4 _dlamc4; DLAMC5 _dlamc5;
655:
656: #endregion
657:
658:
659: #region Fields
660:
661: bool FIRST = false; bool IEEE = false; bool IWARN = false; bool LIEEE1 = false; bool LRND = false; int GNMIN = 0;
662: int GPMIN = 0;int I = 0; int LBETA = 0; int LEMAX = 0; int LEMIN = 0; int LT = 0; int NGNMIN = 0; int NGPMIN = 0;
663: double A = 0;double B = 0; double C = 0; double HALF = 0; double LEPS = 0; double LRMAX = 0; double LRMIN = 0;
664: double ONE = 0;double RBASE = 0; double SIXTH = 0; double SMALL = 0; double THIRD = 0; double TWO = 0; double ZERO = 0;
665:
666: #endregion
667:
668: public DLAMC2(DLAMC3 dlamc3, DLAMC1 dlamc1, DLAMC4 dlamc4, DLAMC5 dlamc5)
669: {
670:
671:
672: #region Set Dependencies
673:
674: this._dlamc3 = dlamc3; this._dlamc1 = dlamc1; this._dlamc4 = dlamc4; this._dlamc5 = dlamc5;
675:
676: #endregion
677:
678:
679: #region Data Inicializacion
680:
681: //FIRST/.TRUE.
682: FIRST = true;
683: //IWARN/.FALSE.
684: IWARN = false;
685:
686: #endregion
687:
688: }
689:
690: public DLAMC2()
691: {
692:
693:
694: #region Dependencies (Initialization)
695:
696: DLAMC3 dlamc3 = new DLAMC3();
697: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
698: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
699: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
700:
701: #endregion
702:
703:
704: #region Set Dependencies
705:
706: this._dlamc3 = dlamc3; this._dlamc1 = dlamc1; this._dlamc4 = dlamc4; this._dlamc5 = dlamc5;
707:
708: #endregion
709:
710:
711: #region Data Inicializacion
712:
713: //FIRST/.TRUE.
714: FIRST = true;
715: //IWARN/.FALSE.
716: IWARN = false;
717:
718: #endregion
719:
720: }
721: /// <summary>
722: /// Purpose
723: /// =======
724: ///
725: /// DLAMC2 determines the machine parameters specified in its argument
726: /// list.
727: ///
728: ///</summary>
729: /// <param name="BETA">
730: /// (output) INTEGER
731: /// The base of the machine.
732: ///</param>
733: /// <param name="T">
734: /// (output) INTEGER
735: /// The number of ( BETA ) digits in the mantissa.
736: ///</param>
737: /// <param name="RND">
738: /// (output) LOGICAL
739: /// Specifies whether proper rounding ( RND = .TRUE. ) or
740: /// chopping ( RND = .FALSE. ) occurs in addition. This may not
741: /// be a reliable guide to the way in which the machine performs
742: /// its arithmetic.
743: ///</param>
744: /// <param name="EPS">
745: /// (output) DOUBLE PRECISION
746: /// The smallest positive number such that
747: ///
748: /// fl( 1.0 - EPS ) .LT. 1.0,
749: ///
750: /// where fl denotes the computed value.
751: ///</param>
752: /// <param name="EMIN">
753: /// (output) INTEGER
754: /// The minimum exponent before (gradual) underflow occurs.
755: ///</param>
756: /// <param name="RMIN">
757: /// (output) DOUBLE PRECISION
758: /// The smallest normalized number for the machine, given by
759: /// BASE**( EMIN - 1 ), where BASE is the floating point value
760: /// of BETA.
761: ///</param>
762: /// <param name="EMAX">
763: /// (output) INTEGER
764: /// The maximum exponent before overflow occurs.
765: ///</param>
766: /// <param name="RMAX">
767: /// (output) DOUBLE PRECISION
768: /// The largest positive number for the machine, given by
769: /// BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
770: /// value of BETA.
771: ///</param>
772: public void Run(ref int BETA, ref int T, ref bool RND, ref double EPS, ref int EMIN, ref double RMIN
773: , ref int EMAX, ref double RMAX)
774: {
775:
776: #region Prolog
777:
778: // *
779: // * -- LAPACK auxiliary routine (version 3.1) --
780: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
781: // * November 2006
782: // *
783: // * .. Scalar Arguments ..
784: // * ..
785: // *
786: // * Purpose
787: // * =======
788: // *
789: // * DLAMC2 determines the machine parameters specified in its argument
790: // * list.
791: // *
792: // * Arguments
793: // * =========
794: // *
795: // * BETA (output) INTEGER
796: // * The base of the machine.
797: // *
798: // * T (output) INTEGER
799: // * The number of ( BETA ) digits in the mantissa.
800: // *
801: // * RND (output) LOGICAL
802: // * Specifies whether proper rounding ( RND = .TRUE. ) or
803: // * chopping ( RND = .FALSE. ) occurs in addition. This may not
804: // * be a reliable guide to the way in which the machine performs
805: // * its arithmetic.
806: // *
807: // * EPS (output) DOUBLE PRECISION
808: // * The smallest positive number such that
809: // *
810: // * fl( 1.0 - EPS ) .LT. 1.0,
811: // *
812: // * where fl denotes the computed value.
813: // *
814: // * EMIN (output) INTEGER
815: // * The minimum exponent before (gradual) underflow occurs.
816: // *
817: // * RMIN (output) DOUBLE PRECISION
818: // * The smallest normalized number for the machine, given by
819: // * BASE**( EMIN - 1 ), where BASE is the floating point value
820: // * of BETA.
821: // *
822: // * EMAX (output) INTEGER
823: // * The maximum exponent before overflow occurs.
824: // *
825: // * RMAX (output) DOUBLE PRECISION
826: // * The largest positive number for the machine, given by
827: // * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
828: // * value of BETA.
829: // *
830: // * Further Details
831: // * ===============
832: // *
833: // * The computation of EPS is based on a routine PARANOIA by
834: // * W. Kahan of the University of California at Berkeley.
835: // *
836: // * =====================================================================
837: // *
838: // * .. Local Scalars ..
839: // * ..
840: // * .. External Functions ..
841: // * ..
842: // * .. External Subroutines ..
843: // * ..
844: // * .. Intrinsic Functions ..
845: // INTRINSIC ABS, MAX, MIN;
846: // * ..
847: // * .. Save statement ..
848: // * ..
849: // * .. Data statements ..
850: // * ..
851: // * .. Executable Statements ..
852: // *
853:
854: #endregion
855:
856:
857: #region Body
858:
859: if (FIRST)
860: {
861: ZERO = 0;
862: ONE = 1;
863: TWO = 2;
864: // *
865: // * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
866: // * BETA, T, RND, EPS, EMIN and RMIN.
867: // *
868: // * Throughout this routine we use the function DLAMC3 to ensure
869: // * that relevant values are stored and not held in registers, or
870: // * are not affected by optimizers.
871: // *
872: // * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
873: // *
874: this._dlamc1.Run(ref LBETA, ref LT, ref LRND, ref LIEEE1);
875: // *
876: // * Start to find EPS.
877: // *
878: B = LBETA;
879: A = Math.Pow(B, - LT);
880: LEPS = A;
881: // *
882: // * Try some tricks to see whether or not this is the correct EPS.
883: // *
884: B = TWO / 3;
885: HALF = ONE / 2;
886: SIXTH = this._dlamc3.Run(B, - HALF);
887: THIRD = this._dlamc3.Run(SIXTH, SIXTH);
888: B = this._dlamc3.Run(THIRD, - HALF);
889: B = this._dlamc3.Run(B, SIXTH);
890: B = Math.Abs(B);
891: if (B < LEPS) B = LEPS;
892: // *
893: LEPS = 1;
894: // *
895: // *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
896: LABEL10:;
897: if ((LEPS > B) && (B > ZERO))
898: {
899: LEPS = B;
900: C = this._dlamc3.Run(HALF * LEPS, (Math.Pow(TWO,5)) * (Math.Pow(LEPS,2)));
901: C = this._dlamc3.Run(HALF, - C);
902: B = this._dlamc3.Run(HALF, C);
903: C = this._dlamc3.Run(HALF, - B);
904: B = this._dlamc3.Run(HALF, C);
905: goto LABEL10;
906: }
907: // *+ END WHILE
908: // *
909: if (A < LEPS) LEPS = A;
910: // *
911: // * Computation of EPS complete.
912: // *
913: // * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
914: // * Keep dividing A by BETA until (gradual) underflow occurs. This
915: // * is detected when we cannot recover the previous A.
916: // *
917: RBASE = ONE / LBETA;
918: SMALL = ONE;
919: for (I = 1; I <= 3; I++)
920: {
921: SMALL = this._dlamc3.Run(SMALL * RBASE, ZERO);
922: }
923: A = this._dlamc3.Run(ONE, SMALL);
924: this._dlamc4.Run(ref NGPMIN, ONE, LBETA);
925: this._dlamc4.Run(ref NGNMIN, - ONE, LBETA);
926: this._dlamc4.Run(ref GPMIN, A, LBETA);
927: this._dlamc4.Run(ref GNMIN, - A, LBETA);
928: IEEE = false;
929: // *
930: if ((NGPMIN == NGNMIN) && (GPMIN == GNMIN))
931: {
932: if (NGPMIN == GPMIN)
933: {
934: LEMIN = NGPMIN;
935: // * ( Non twos-complement machines, no gradual underflow;
936: // * e.g., VAX )
937: }
938: else
939: {
940: if ((GPMIN - NGPMIN) == 3)
941: {
942: LEMIN = NGPMIN - 1 + LT;
943: IEEE = true;
944: // * ( Non twos-complement machines, with gradual underflow;
945: // * e.g., IEEE standard followers )
946: }
947: else
948: {
949: LEMIN = Math.Min(NGPMIN, GPMIN);
950: // * ( A guess; no known machine )
951: IWARN = true;
952: }
953: }
954: // *
955: }
956: else
957: {
958: if ((NGPMIN == GPMIN) && (NGNMIN == GNMIN))
959: {
960: if (Math.Abs(NGPMIN - NGNMIN) == 1)
961: {
962: LEMIN = Math.Max(NGPMIN, NGNMIN);
963: // * ( Twos-complement machines, no gradual underflow;
964: // * e.g., CYBER 205 )
965: }
966: else
967: {
968: LEMIN = Math.Min(NGPMIN, NGNMIN);
969: // * ( A guess; no known machine )
970: IWARN = true;
971: }
972: // *
973: }
974: else
975: {
976: if ((Math.Abs(NGPMIN - NGNMIN) == 1) && (GPMIN == GNMIN))
977: {
978: if ((GPMIN - Math.Min(NGPMIN, NGNMIN)) == 3)
979: {
980: LEMIN = Math.Max(NGPMIN, NGNMIN) - 1 + LT;
981: // * ( Twos-complement machines with gradual underflow;
982: // * no known machine )
983: }
984: else
985: {
986: LEMIN = Math.Min(NGPMIN, NGNMIN);
987: // * ( A guess; no known machine )
988: IWARN = true;
989: }
990: // *
991: }
992: else
993: {
994: LEMIN = Math.Min(NGPMIN, Math.Min(NGNMIN, Math.Min(GPMIN, GNMIN)));
995: // * ( A guess; no known machine )
996: IWARN = true;
997: }
998: }
999: }
1000: FIRST = false;
1001: // ***
1002: // * Comment out this if block if EMIN is ok
1003: if (IWARN)
1004: {
1005: FIRST = true;
1006: //ERROR-ERROR WRITE( 6, FMT = 9999 )LEMIN;
1007: }
1008: // ***
1009: // *
1010: // * Assume IEEE arithmetic if we found denormalised numbers above,
1011: // * or if arithmetic seems to round in the IEEE style, determined
1012: // * in routine DLAMC1. A true IEEE machine should have both things
1013: // * true; however, faulty machines may have one or the other.
1014: // *
1015: IEEE = IEEE || LIEEE1;
1016: // *
1017: // * Compute RMIN by successive division by BETA. We could compute
1018: // * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
1019: // * this computation.
1020: // *
1021: LRMIN = 1;
1022: for (I = 1; I <= 1 - LEMIN; I++)
1023: {
1024: LRMIN = this._dlamc3.Run(LRMIN * RBASE, ZERO);
1025: }
1026: // *
1027: // * Finally, call DLAMC5 to compute EMAX and RMAX.
1028: // *
1029: this._dlamc5.Run(LBETA, LT, LEMIN, IEEE, ref LEMAX, ref LRMAX);
1030: }
1031: // *
1032: BETA = LBETA;
1033: T = LT;
1034: RND = LRND;
1035: EPS = LEPS;
1036: EMIN = LEMIN;
1037: RMIN = LRMIN;
1038: EMAX = LEMAX;
1039: RMAX = LRMAX;
1040: // *
1041: return;
1042: // *
1043: // *
1044: // * End of DLAMC2
1045: // *
1046:
1047: #endregion
1048:
1049: }
1050: }
1051:
1052: #endregion
1053:
1054:
1055: #region The Class: DLAMC3
1056:
1057: // *
1058: // ************************************************************************
1059: // *
1060: /// <summary>
1061: /// -- LAPACK auxiliary routine (version 3.1) --
1062: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1063: /// November 2006
1064: /// Purpose
1065: /// =======
1066: ///
1067: /// DLAMC3 is intended to force A and B to be stored prior to doing
1068: /// the addition of A and B , for use in situations where optimizers
1069: /// might hold one of these in a register.
1070: ///
1071: ///</summary>
1072: public class DLAMC3
1073: {
1074:
1075: public DLAMC3()
1076: {
1077:
1078: }
1079:
1080: /// <summary>
1081: /// Purpose
1082: /// =======
1083: ///
1084: /// DLAMC3 is intended to force A and B to be stored prior to doing
1085: /// the addition of A and B , for use in situations where optimizers
1086: /// might hold one of these in a register.
1087: ///
1088: ///</summary>
1089: /// <param name="A">
1090: /// (input) DOUBLE PRECISION
1091: ///</param>
1092: /// <param name="B">
1093: /// (input) DOUBLE PRECISION
1094: /// The values A and B.
1095: ///</param>
1096: public double Run(double A, double B)
1097: {
1098: double dlamc3 = 0;
1099:
1100: #region Prolog
1101:
1102: // *
1103: // * -- LAPACK auxiliary routine (version 3.1) --
1104: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1105: // * November 2006
1106: // *
1107: // * .. Scalar Arguments ..
1108: // * ..
1109: // *
1110: // * Purpose
1111: // * =======
1112: // *
1113: // * DLAMC3 is intended to force A and B to be stored prior to doing
1114: // * the addition of A and B , for use in situations where optimizers
1115: // * might hold one of these in a register.
1116: // *
1117: // * Arguments
1118: // * =========
1119: // *
1120: // * A (input) DOUBLE PRECISION
1121: // * B (input) DOUBLE PRECISION
1122: // * The values A and B.
1123: // *
1124: // * =====================================================================
1125: // *
1126: // * .. Executable Statements ..
1127: // *
1128:
1129: #endregion
1130:
1131: dlamc3 = A + B;
1132: // *
1133: return dlamc3;
1134: // *
1135: // * End of DLAMC3
1136: // *
1137: }
1138: }
1139:
1140: #endregion
1141:
1142:
1143: #region The Class: DLAMC4
1144:
1145: // *
1146: // ************************************************************************
1147: // *
1148: /// <summary>
1149: /// -- LAPACK auxiliary routine (version 3.1) --
1150: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1151: /// November 2006
1152: /// Purpose
1153: /// =======
1154: ///
1155: /// DLAMC4 is a service routine for DLAMC2.
1156: ///
1157: ///</summary>
1158: public class DLAMC4
1159: {
1160:
1161:
1162: #region Dependencies
1163:
1164: DLAMC3 _dlamc3;
1165:
1166: #endregion
1167:
1168:
1169: #region Fields
1170:
1171: int I = 0; double A = 0; double B1 = 0; double B2 = 0; double C1 = 0; double C2 = 0; double D1 = 0; double D2 = 0;
1172: double ONE = 0;double RBASE = 0; double ZERO = 0;
1173:
1174: #endregion
1175:
1176: public DLAMC4(DLAMC3 dlamc3)
1177: {
1178:
1179:
1180: #region Set Dependencies
1181:
1182: this._dlamc3 = dlamc3;
1183:
1184: #endregion
1185:
1186: }
1187:
1188: public DLAMC4()
1189: {
1190:
1191:
1192: #region Dependencies (Initialization)
1193:
1194: DLAMC3 dlamc3 = new DLAMC3();
1195:
1196: #endregion
1197:
1198:
1199: #region Set Dependencies
1200:
1201: this._dlamc3 = dlamc3;
1202:
1203: #endregion
1204:
1205: }
1206: /// <summary>
1207: /// Purpose
1208: /// =======
1209: ///
1210: /// DLAMC4 is a service routine for DLAMC2.
1211: ///
1212: ///</summary>
1213: /// <param name="EMIN">
1214: /// (output) INTEGER
1215: /// The minimum exponent before (gradual) underflow, computed by
1216: /// setting A = START and dividing by BASE until the previous A
1217: /// can not be recovered.
1218: ///</param>
1219: /// <param name="START">
1220: /// (input) DOUBLE PRECISION
1221: /// The starting point for determining EMIN.
1222: ///</param>
1223: /// <param name="BASE">
1224: /// (input) INTEGER
1225: /// The base of the machine.
1226: ///</param>
1227: public void Run(ref int EMIN, double START, int BASE)
1228: {
1229:
1230: #region Prolog
1231:
1232: // *
1233: // * -- LAPACK auxiliary routine (version 3.1) --
1234: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1235: // * November 2006
1236: // *
1237: // * .. Scalar Arguments ..
1238: // * ..
1239: // *
1240: // * Purpose
1241: // * =======
1242: // *
1243: // * DLAMC4 is a service routine for DLAMC2.
1244: // *
1245: // * Arguments
1246: // * =========
1247: // *
1248: // * EMIN (output) INTEGER
1249: // * The minimum exponent before (gradual) underflow, computed by
1250: // * setting A = START and dividing by BASE until the previous A
1251: // * can not be recovered.
1252: // *
1253: // * START (input) DOUBLE PRECISION
1254: // * The starting point for determining EMIN.
1255: // *
1256: // * BASE (input) INTEGER
1257: // * The base of the machine.
1258: // *
1259: // * =====================================================================
1260: // *
1261: // * .. Local Scalars ..
1262: // * ..
1263: // * .. External Functions ..
1264: // * ..
1265: // * .. Executable Statements ..
1266: // *
1267:
1268: #endregion
1269:
1270:
1271: #region Body
1272:
1273: A = START;
1274: ONE = 1;
1275: RBASE = ONE / BASE;
1276: ZERO = 0;
1277: EMIN = 1;
1278: B1 = this._dlamc3.Run(A * RBASE, ZERO);
1279: C1 = A;
1280: C2 = A;
1281: D1 = A;
1282: D2 = A;
1283: // *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
1284: // * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
1285: LABEL10:;
1286: if ((C1 == A) && (C2 == A) && (D1 == A) && (D2 == A))
1287: {
1288: EMIN = EMIN - 1;
1289: A = B1;
1290: B1 = this._dlamc3.Run(A / BASE, ZERO);
1291: C1 = this._dlamc3.Run(B1 * BASE, ZERO);
1292: D1 = ZERO;
1293: for (I = 1; I <= BASE; I++)
1294: {
1295: D1 = D1 + B1;
1296: }
1297: B2 = this._dlamc3.Run(A * RBASE, ZERO);
1298: C2 = this._dlamc3.Run(B2 / RBASE, ZERO);
1299: D2 = ZERO;
1300: for (I = 1; I <= BASE; I++)
1301: {
1302: D2 = D2 + B2;
1303: }
1304: goto LABEL10;
1305: }
1306: // *+ END WHILE
1307: // *
1308: return;
1309: // *
1310: // * End of DLAMC4
1311: // *
1312:
1313: #endregion
1314:
1315: }
1316: }
1317:
1318: #endregion
1319:
1320:
1321: #region The Class: DLAMC5
1322:
1323: // *
1324: // ************************************************************************
1325: // *
1326: /// <summary>
1327: /// -- LAPACK auxiliary routine (version 3.1) --
1328: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1329: /// November 2006
1330: /// Purpose
1331: /// =======
1332: ///
1333: /// DLAMC5 attempts to compute RMAX, the largest machine floating-point
1334: /// number, without overflow. It assumes that EMAX + abs(EMIN) sum
1335: /// approximately to a power of 2. It will fail on machines where this
1336: /// assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
1337: /// EMAX = 28718). It will also fail if the value supplied for EMIN is
1338: /// too large (i.e. too close to zero), probably with overflow.
1339: ///
1340: ///</summary>
1341: public class DLAMC5
1342: {
1343:
1344:
1345: #region Dependencies
1346:
1347: DLAMC3 _dlamc3;
1348:
1349: #endregion
1350:
1351:
1352: #region Fields
1353:
1354: const double ZERO = 0.0E0; const double ONE = 1.0E0; int EXBITS = 0; int EXPSUM = 0; int I = 0; int LEXP = 0;
1355: int NBITS = 0;int TRY = 0; int UEXP = 0; double OLDY = 0; double RECBAS = 0; double Y = 0; double Z = 0;
1356:
1357: #endregion
1358:
1359: public DLAMC5(DLAMC3 dlamc3)
1360: {
1361:
1362:
1363: #region Set Dependencies
1364:
1365: this._dlamc3 = dlamc3;
1366:
1367: #endregion
1368:
1369: }
1370:
1371: public DLAMC5()
1372: {
1373:
1374:
1375: #region Dependencies (Initialization)
1376:
1377: DLAMC3 dlamc3 = new DLAMC3();
1378:
1379: #endregion
1380:
1381:
1382: #region Set Dependencies
1383:
1384: this._dlamc3 = dlamc3;
1385:
1386: #endregion
1387:
1388: }
1389: /// <summary>
1390: /// Purpose
1391: /// =======
1392: ///
1393: /// DLAMC5 attempts to compute RMAX, the largest machine floating-point
1394: /// number, without overflow. It assumes that EMAX + abs(EMIN) sum
1395: /// approximately to a power of 2. It will fail on machines where this
1396: /// assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
1397: /// EMAX = 28718). It will also fail if the value supplied for EMIN is
1398: /// too large (i.e. too close to zero), probably with overflow.
1399: ///
1400: ///</summary>
1401: /// <param name="BETA">
1402: /// (input) INTEGER
1403: /// The base of floating-point arithmetic.
1404: ///</param>
1405: /// <param name="P">
1406: /// (input) INTEGER
1407: /// The number of base BETA digits in the mantissa of a
1408: /// floating-point value.
1409: ///</param>
1410: /// <param name="EMIN">
1411: /// (input) INTEGER
1412: /// The minimum exponent before (gradual) underflow.
1413: ///</param>
1414: /// <param name="IEEE">
1415: /// (input) LOGICAL
1416: /// A logical flag specifying whether or not the arithmetic
1417: /// system is thought to comply with the IEEE standard.
1418: ///</param>
1419: /// <param name="EMAX">
1420: /// (output) INTEGER
1421: /// The largest exponent before overflow
1422: ///</param>
1423: /// <param name="RMAX">
1424: /// (output) DOUBLE PRECISION
1425: /// The largest machine floating-point number.
1426: ///</param>
1427: public void Run(int BETA, int P, int EMIN, bool IEEE, ref int EMAX, ref double RMAX)
1428: {
1429:
1430: #region Prolog
1431:
1432: // *
1433: // * -- LAPACK auxiliary routine (version 3.1) --
1434: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
1435: // * November 2006
1436: // *
1437: // * .. Scalar Arguments ..
1438: // * ..
1439: // *
1440: // * Purpose
1441: // * =======
1442: // *
1443: // * DLAMC5 attempts to compute RMAX, the largest machine floating-point
1444: // * number, without overflow. It assumes that EMAX + abs(EMIN) sum
1445: // * approximately to a power of 2. It will fail on machines where this
1446: // * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
1447: // * EMAX = 28718). It will also fail if the value supplied for EMIN is
1448: // * too large (i.e. too close to zero), probably with overflow.
1449: // *
1450: // * Arguments
1451: // * =========
1452: // *
1453: // * BETA (input) INTEGER
1454: // * The base of floating-point arithmetic.
1455: // *
1456: // * P (input) INTEGER
1457: // * The number of base BETA digits in the mantissa of a
1458: // * floating-point value.
1459: // *
1460: // * EMIN (input) INTEGER
1461: // * The minimum exponent before (gradual) underflow.
1462: // *
1463: // * IEEE (input) LOGICAL
1464: // * A logical flag specifying whether or not the arithmetic
1465: // * system is thought to comply with the IEEE standard.
1466: // *
1467: // * EMAX (output) INTEGER
1468: // * The largest exponent before overflow
1469: // *
1470: // * RMAX (output) DOUBLE PRECISION
1471: // * The largest machine floating-point number.
1472: // *
1473: // * =====================================================================
1474: // *
1475: // * .. Parameters ..
1476: // * ..
1477: // * .. Local Scalars ..
1478: // * ..
1479: // * .. External Functions ..
1480: // * ..
1481: // * .. Intrinsic Functions ..
1482: // INTRINSIC MOD;
1483: // * ..
1484: // * .. Executable Statements ..
1485: // *
1486: // * First compute LEXP and UEXP, two powers of 2 that bound
1487: // * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
1488: // * approximately to the bound that is closest to abs(EMIN).
1489: // * (EMAX is the exponent of the required number RMAX).
1490: // *
1491:
1492: #endregion
1493:
1494:
1495: #region Body
1496:
1497: LEXP = 1;
1498: EXBITS = 1;
1499: LABEL10:;
1500: TRY = LEXP * 2;
1501: if (TRY <= ( - EMIN))
1502: {
1503: LEXP = TRY;
1504: EXBITS = EXBITS + 1;
1505: goto LABEL10;
1506: }
1507: if (LEXP == - EMIN)
1508: {
1509: UEXP = LEXP;
1510: }
1511: else
1512: {
1513: UEXP = TRY;
1514: EXBITS = EXBITS + 1;
1515: }
1516: // *
1517: // * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
1518: // * than or equal to EMIN. EXBITS is the number of bits needed to
1519: // * store the exponent.
1520: // *
1521: if ((UEXP + EMIN) > ( - LEXP - EMIN))
1522: {
1523: EXPSUM = 2 * LEXP;
1524: }
1525: else
1526: {
1527: EXPSUM = 2 * UEXP;
1528: }
1529: // *
1530: // * EXPSUM is the exponent range, approximately equal to
1531: // * EMAX - EMIN + 1 .
1532: // *
1533: EMAX = EXPSUM + EMIN - 1;
1534: NBITS = 1 + EXBITS + P;
1535: // *
1536: // * NBITS is the total number of bits needed to store a
1537: // * floating-point number.
1538: // *
1539: if ((FortranLib.Mod(NBITS,2) == 1) && (BETA == 2))
1540: {
1541: // *
1542: // * Either there are an odd number of bits used to store a
1543: // * floating-point number, which is unlikely, or some bits are
1544: // * not used in the representation of numbers, which is possible,
1545: // * (e.g. Cray machines) or the mantissa has an implicit bit,
1546: // * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
1547: // * most likely. We have to assume the last alternative.
1548: // * If this is true, then we need to reduce EMAX by one because
1549: // * there must be some way of representing zero in an implicit-bit
1550: // * system. On machines like Cray, we are reducing EMAX by one
1551: // * unnecessarily.
1552: // *
1553: EMAX = EMAX - 1;
1554: }
1555: // *
1556: if (IEEE)
1557: {
1558: // *
1559: // * Assume we are on an IEEE machine which reserves one exponent
1560: // * for infinity and NaN.
1561: // *
1562: EMAX = EMAX - 1;
1563: }
1564: // *
1565: // * Now create RMAX, the largest machine number, which should
1566: // * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
1567: // *
1568: // * First compute 1.0 - BETA**(-P), being careful that the
1569: // * result is less than 1.0 .
1570: // *
1571: RECBAS = ONE / BETA;
1572: Z = BETA - ONE;
1573: Y = ZERO;
1574: for (I = 1; I <= P; I++)
1575: {
1576: Z = Z * RECBAS;
1577: if (Y < ONE) OLDY = Y;
1578: Y = this._dlamc3.Run(Y, Z);
1579: }
1580: if (Y >= ONE) Y = OLDY;
1581: // *
1582: // * Now multiply by BETA**EMAX to get RMAX.
1583: // *
1584: for (I = 1; I <= EMAX; I++)
1585: {
1586: Y = this._dlamc3.Run(Y * BETA, ZERO);
1587: }
1588: // *
1589: RMAX = Y;
1590: return;
1591: // *
1592: // * End of DLAMC5
1593: // *
1594:
1595: #endregion
1596:
1597: }
1598: }
1599:
1600: #endregion
1601:
1602: }