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: /// <summary>
20: /// -- LAPACK auxiliary routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
27: /// [ A B ]
28: /// [ B C ].
29: /// On return, RT1 is the eigenvalue of larger absolute value, and RT2
30: /// is the eigenvalue of smaller absolute value.
31: ///
32: ///</summary>
33: public class DLAE2
34: {
35:
36:
37: #region Fields
38:
39: const double ONE = 1.0E0; const double TWO = 2.0E0; const double ZERO = 0.0E0; const double HALF = 0.5E0; double AB = 0;
40: double ACMN = 0;double ACMX = 0; double ADF = 0; double DF = 0; double RT = 0; double SM = 0; double TB = 0;
41:
42: #endregion
43:
44: public DLAE2()
45: {
46:
47: }
48:
49: /// <summary>
50: /// Purpose
51: /// =======
52: ///
53: /// DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
54: /// [ A B ]
55: /// [ B C ].
56: /// On return, RT1 is the eigenvalue of larger absolute value, and RT2
57: /// is the eigenvalue of smaller absolute value.
58: ///
59: ///</summary>
60: /// <param name="A">
61: /// (input) DOUBLE PRECISION
62: /// The (1,1) element of the 2-by-2 matrix.
63: ///</param>
64: /// <param name="B">
65: /// (input) DOUBLE PRECISION
66: /// The (1,2) and (2,1) elements of the 2-by-2 matrix.
67: ///</param>
68: /// <param name="C">
69: /// (input) DOUBLE PRECISION
70: /// The (2,2) element of the 2-by-2 matrix.
71: ///</param>
72: /// <param name="RT1">
73: /// (output) DOUBLE PRECISION
74: /// The eigenvalue of larger absolute value.
75: ///</param>
76: /// <param name="RT2">
77: /// (output) DOUBLE PRECISION
78: /// The eigenvalue of smaller absolute value.
79: ///</param>
80: public void Run(double A, double B, double C, ref double RT1, ref double RT2)
81: {
82:
83: #region Prolog
84:
85: // *
86: // * -- LAPACK auxiliary routine (version 3.1) --
87: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
88: // * November 2006
89: // *
90: // * .. Scalar Arguments ..
91: // * ..
92: // *
93: // * Purpose
94: // * =======
95: // *
96: // * DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
97: // * [ A B ]
98: // * [ B C ].
99: // * On return, RT1 is the eigenvalue of larger absolute value, and RT2
100: // * is the eigenvalue of smaller absolute value.
101: // *
102: // * Arguments
103: // * =========
104: // *
105: // * A (input) DOUBLE PRECISION
106: // * The (1,1) element of the 2-by-2 matrix.
107: // *
108: // * B (input) DOUBLE PRECISION
109: // * The (1,2) and (2,1) elements of the 2-by-2 matrix.
110: // *
111: // * C (input) DOUBLE PRECISION
112: // * The (2,2) element of the 2-by-2 matrix.
113: // *
114: // * RT1 (output) DOUBLE PRECISION
115: // * The eigenvalue of larger absolute value.
116: // *
117: // * RT2 (output) DOUBLE PRECISION
118: // * The eigenvalue of smaller absolute value.
119: // *
120: // * Further Details
121: // * ===============
122: // *
123: // * RT1 is accurate to a few ulps barring over/underflow.
124: // *
125: // * RT2 may be inaccurate if there is massive cancellation in the
126: // * determinant A*C-B*B; higher precision or correctly rounded or
127: // * correctly truncated arithmetic would be needed to compute RT2
128: // * accurately in all cases.
129: // *
130: // * Overflow is possible only if RT1 is within a factor of 5 of overflow.
131: // * Underflow is harmless if the input data is 0 or exceeds
132: // * underflow_threshold / macheps.
133: // *
134: // * =====================================================================
135: // *
136: // * .. Parameters ..
137: // * ..
138: // * .. Local Scalars ..
139: // * ..
140: // * .. Intrinsic Functions ..
141: // INTRINSIC ABS, SQRT;
142: // * ..
143: // * .. Executable Statements ..
144: // *
145: // * Compute the eigenvalues
146: // *
147:
148: #endregion
149:
150:
151: #region Body
152:
153: SM = A + C;
154: DF = A - C;
155: ADF = Math.Abs(DF);
156: TB = B + B;
157: AB = Math.Abs(TB);
158: if (Math.Abs(A) > Math.Abs(C))
159: {
160: ACMX = A;
161: ACMN = C;
162: }
163: else
164: {
165: ACMX = C;
166: ACMN = A;
167: }
168: if (ADF > AB)
169: {
170: RT = ADF * Math.Sqrt(ONE + Math.Pow(AB / ADF,2));
171: }
172: else
173: {
174: if (ADF < AB)
175: {
176: RT = AB * Math.Sqrt(ONE + Math.Pow(ADF / AB,2));
177: }
178: else
179: {
180: // *
181: // * Includes case AB=ADF=0
182: // *
183: RT = AB * Math.Sqrt(TWO);
184: }
185: }
186: if (SM < ZERO)
187: {
188: RT1 = HALF * (SM - RT);
189: // *
190: // * Order of execution important.
191: // * To get fully accurate smaller eigenvalue,
192: // * next line needs to be executed in higher precision.
193: // *
194: RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
195: }
196: else
197: {
198: if (SM > ZERO)
199: {
200: RT1 = HALF * (SM + RT);
201: // *
202: // * Order of execution important.
203: // * To get fully accurate smaller eigenvalue,
204: // * next line needs to be executed in higher precision.
205: // *
206: RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
207: }
208: else
209: {
210: // *
211: // * Includes case RT1 = RT2 = 0
212: // *
213: RT1 = HALF * RT;
214: RT2 = - HALF * RT;
215: }
216: }
217: return;
218: // *
219: // * End of DLAE2
220: // *
221:
222: #endregion
223:
224: }
225: }
226: }