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: /// DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
27: /// [ A B ]
28: /// [ B C ].
29: /// On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
30: /// eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
31: /// eigenvector for RT1, giving the decomposition
32: ///
33: /// [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
34: /// [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
35: ///
36: ///</summary>
37: public class DLAEV2
38: {
39:
40:
41: #region Fields
42:
43: const double ONE = 1.0E0; const double TWO = 2.0E0; const double ZERO = 0.0E0; const double HALF = 0.5E0; int SGN1 = 0;
44: int SGN2 = 0;double AB = 0; double ACMN = 0; double ACMX = 0; double ACS = 0; double ADF = 0; double CS = 0;
45: double CT = 0;double DF = 0; double RT = 0; double SM = 0; double TB = 0; double TN = 0;
46:
47: #endregion
48:
49: public DLAEV2()
50: {
51:
52: }
53:
54: /// <summary>
55: /// Purpose
56: /// =======
57: ///
58: /// DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
59: /// [ A B ]
60: /// [ B C ].
61: /// On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
62: /// eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
63: /// eigenvector for RT1, giving the decomposition
64: ///
65: /// [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
66: /// [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
67: ///
68: ///</summary>
69: /// <param name="A">
70: /// (input) DOUBLE PRECISION
71: /// The (1,1) element of the 2-by-2 matrix.
72: ///</param>
73: /// <param name="B">
74: /// (input) DOUBLE PRECISION
75: /// The (1,2) element and the conjugate of the (2,1) element of
76: /// the 2-by-2 matrix.
77: ///</param>
78: /// <param name="C">
79: /// (input) DOUBLE PRECISION
80: /// The (2,2) element of the 2-by-2 matrix.
81: ///</param>
82: /// <param name="RT1">
83: /// (output) DOUBLE PRECISION
84: /// The eigenvalue of larger absolute value.
85: ///</param>
86: /// <param name="RT2">
87: /// (output) DOUBLE PRECISION
88: /// The eigenvalue of smaller absolute value.
89: ///</param>
90: /// <param name="CS1">
91: /// (output) DOUBLE PRECISION
92: ///</param>
93: /// <param name="SN1">
94: /// (output) DOUBLE PRECISION
95: /// The vector (CS1, SN1) is a unit right eigenvector for RT1.
96: ///</param>
97: public void Run(double A, double B, double C, ref double RT1, ref double RT2, ref double CS1
98: , ref double SN1)
99: {
100:
101: #region Prolog
102:
103: // *
104: // * -- LAPACK auxiliary routine (version 3.1) --
105: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
106: // * November 2006
107: // *
108: // * .. Scalar Arguments ..
109: // * ..
110: // *
111: // * Purpose
112: // * =======
113: // *
114: // * DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
115: // * [ A B ]
116: // * [ B C ].
117: // * On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
118: // * eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
119: // * eigenvector for RT1, giving the decomposition
120: // *
121: // * [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
122: // * [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
123: // *
124: // * Arguments
125: // * =========
126: // *
127: // * A (input) DOUBLE PRECISION
128: // * The (1,1) element of the 2-by-2 matrix.
129: // *
130: // * B (input) DOUBLE PRECISION
131: // * The (1,2) element and the conjugate of the (2,1) element of
132: // * the 2-by-2 matrix.
133: // *
134: // * C (input) DOUBLE PRECISION
135: // * The (2,2) element of the 2-by-2 matrix.
136: // *
137: // * RT1 (output) DOUBLE PRECISION
138: // * The eigenvalue of larger absolute value.
139: // *
140: // * RT2 (output) DOUBLE PRECISION
141: // * The eigenvalue of smaller absolute value.
142: // *
143: // * CS1 (output) DOUBLE PRECISION
144: // * SN1 (output) DOUBLE PRECISION
145: // * The vector (CS1, SN1) is a unit right eigenvector for RT1.
146: // *
147: // * Further Details
148: // * ===============
149: // *
150: // * RT1 is accurate to a few ulps barring over/underflow.
151: // *
152: // * RT2 may be inaccurate if there is massive cancellation in the
153: // * determinant A*C-B*B; higher precision or correctly rounded or
154: // * correctly truncated arithmetic would be needed to compute RT2
155: // * accurately in all cases.
156: // *
157: // * CS1 and SN1 are accurate to a few ulps barring over/underflow.
158: // *
159: // * Overflow is possible only if RT1 is within a factor of 5 of overflow.
160: // * Underflow is harmless if the input data is 0 or exceeds
161: // * underflow_threshold / macheps.
162: // *
163: // * =====================================================================
164: // *
165: // * .. Parameters ..
166: // * ..
167: // * .. Local Scalars ..
168: // * ..
169: // * .. Intrinsic Functions ..
170: // INTRINSIC ABS, SQRT;
171: // * ..
172: // * .. Executable Statements ..
173: // *
174: // * Compute the eigenvalues
175: // *
176:
177: #endregion
178:
179:
180: #region Body
181:
182: SM = A + C;
183: DF = A - C;
184: ADF = Math.Abs(DF);
185: TB = B + B;
186: AB = Math.Abs(TB);
187: if (Math.Abs(A) > Math.Abs(C))
188: {
189: ACMX = A;
190: ACMN = C;
191: }
192: else
193: {
194: ACMX = C;
195: ACMN = A;
196: }
197: if (ADF > AB)
198: {
199: RT = ADF * Math.Sqrt(ONE + Math.Pow(AB / ADF,2));
200: }
201: else
202: {
203: if (ADF < AB)
204: {
205: RT = AB * Math.Sqrt(ONE + Math.Pow(ADF / AB,2));
206: }
207: else
208: {
209: // *
210: // * Includes case AB=ADF=0
211: // *
212: RT = AB * Math.Sqrt(TWO);
213: }
214: }
215: if (SM < ZERO)
216: {
217: RT1 = HALF * (SM - RT);
218: SGN1 = - 1;
219: // *
220: // * Order of execution important.
221: // * To get fully accurate smaller eigenvalue,
222: // * next line needs to be executed in higher precision.
223: // *
224: RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
225: }
226: else
227: {
228: if (SM > ZERO)
229: {
230: RT1 = HALF * (SM + RT);
231: SGN1 = 1;
232: // *
233: // * Order of execution important.
234: // * To get fully accurate smaller eigenvalue,
235: // * next line needs to be executed in higher precision.
236: // *
237: RT2 = (ACMX / RT1) * ACMN - (B / RT1) * B;
238: }
239: else
240: {
241: // *
242: // * Includes case RT1 = RT2 = 0
243: // *
244: RT1 = HALF * RT;
245: RT2 = - HALF * RT;
246: SGN1 = 1;
247: }
248: }
249: // *
250: // * Compute the eigenvector
251: // *
252: if (DF >= ZERO)
253: {
254: CS = DF + RT;
255: SGN2 = 1;
256: }
257: else
258: {
259: CS = DF - RT;
260: SGN2 = - 1;
261: }
262: ACS = Math.Abs(CS);
263: if (ACS > AB)
264: {
265: CT = - TB / CS;
266: SN1 = ONE / Math.Sqrt(ONE + CT * CT);
267: CS1 = CT * SN1;
268: }
269: else
270: {
271: if (AB == ZERO)
272: {
273: CS1 = ONE;
274: SN1 = ZERO;
275: }
276: else
277: {
278: TN = - CS / TB;
279: CS1 = ONE / Math.Sqrt(ONE + TN * TN);
280: SN1 = TN * CS1;
281: }
282: }
283: if (SGN1 == SGN2)
284: {
285: TN = CS1;
286: CS1 = - SN1;
287: SN1 = TN;
288: }
289: return;
290: // *
291: // * End of DLAEV2
292: // *
293:
294: #endregion
295:
296: }
297: }
298: }