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 driver routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
27: /// matrix in standard form:
28: ///
29: /// [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
30: /// [ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
31: ///
32: /// where either
33: /// 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
34: /// 2) AA = DD and BB*CC .LT. 0, so that AA + or - sqrt(BB*CC) are complex
35: /// conjugate eigenvalues.
36: ///
37: ///</summary>
38: public class DLANV2
39: {
40:
41:
42: #region Dependencies
43:
44: DLAMCH _dlamch; DLAPY2 _dlapy2;
45:
46: #endregion
47:
48:
49: #region Fields
50:
51: const double ZERO = 0.0E+0; const double HALF = 0.5E+0; const double ONE = 1.0E+0; const double MULTPL = 4.0E+0;
52: double AA = 0;double BB = 0; double BCMAX = 0; double BCMIS = 0; double CC = 0; double CS1 = 0; double DD = 0;
53: double EPS = 0;double P = 0; double SAB = 0; double SAC = 0; double SCALE = 0; double SIGMA = 0; double SN1 = 0;
54: double TAU = 0;double TEMP = 0; double Z = 0;
55:
56: #endregion
57:
58: public DLANV2(DLAMCH dlamch, DLAPY2 dlapy2)
59: {
60:
61:
62: #region Set Dependencies
63:
64: this._dlamch = dlamch; this._dlapy2 = dlapy2;
65:
66: #endregion
67:
68: }
69:
70: public DLANV2()
71: {
72:
73:
74: #region Dependencies (Initialization)
75:
76: LSAME lsame = new LSAME();
77: DLAMC3 dlamc3 = new DLAMC3();
78: DLAPY2 dlapy2 = new DLAPY2();
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: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
84:
85: #endregion
86:
87:
88: #region Set Dependencies
89:
90: this._dlamch = dlamch; this._dlapy2 = dlapy2;
91:
92: #endregion
93:
94: }
95: /// <summary>
96: /// Purpose
97: /// =======
98: ///
99: /// DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
100: /// matrix in standard form:
101: ///
102: /// [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
103: /// [ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
104: ///
105: /// where either
106: /// 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
107: /// 2) AA = DD and BB*CC .LT. 0, so that AA + or - sqrt(BB*CC) are complex
108: /// conjugate eigenvalues.
109: ///
110: ///</summary>
111: /// <param name="A">
112: /// (input/output) DOUBLE PRECISION
113: ///</param>
114: /// <param name="B">
115: /// (input/output) DOUBLE PRECISION
116: ///</param>
117: /// <param name="C">
118: /// (input/output) DOUBLE PRECISION
119: ///</param>
120: /// <param name="D">
121: /// (input/output) DOUBLE PRECISION
122: /// On entry, the elements of the input matrix.
123: /// On exit, they are overwritten by the elements of the
124: /// standardised Schur form.
125: ///</param>
126: /// <param name="RT1R">
127: /// (output) DOUBLE PRECISION
128: ///</param>
129: /// <param name="RT1I">
130: /// (output) DOUBLE PRECISION
131: ///</param>
132: /// <param name="RT2R">
133: /// (output) DOUBLE PRECISION
134: ///</param>
135: /// <param name="RT2I">
136: /// (output) DOUBLE PRECISION
137: /// The real and imaginary parts of the eigenvalues. If the
138: /// eigenvalues are a complex conjugate pair, RT1I .GT. 0.
139: ///</param>
140: /// <param name="CS">
141: /// (output) DOUBLE PRECISION
142: ///</param>
143: /// <param name="SN">
144: /// (output) DOUBLE PRECISION
145: /// Parameters of the rotation matrix.
146: ///</param>
147: public void Run(ref double A, ref double B, ref double C, ref double D, ref double RT1R, ref double RT1I
148: , ref double RT2R, ref double RT2I, ref double CS, ref double SN)
149: {
150:
151: #region Prolog
152:
153: // *
154: // * -- LAPACK driver routine (version 3.1) --
155: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
156: // * November 2006
157: // *
158: // * .. Scalar Arguments ..
159: // * ..
160: // *
161: // * Purpose
162: // * =======
163: // *
164: // * DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
165: // * matrix in standard form:
166: // *
167: // * [ A B ] = [ CS -SN ] [ AA BB ] [ CS SN ]
168: // * [ C D ] [ SN CS ] [ CC DD ] [-SN CS ]
169: // *
170: // * where either
171: // * 1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
172: // * 2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
173: // * conjugate eigenvalues.
174: // *
175: // * Arguments
176: // * =========
177: // *
178: // * A (input/output) DOUBLE PRECISION
179: // * B (input/output) DOUBLE PRECISION
180: // * C (input/output) DOUBLE PRECISION
181: // * D (input/output) DOUBLE PRECISION
182: // * On entry, the elements of the input matrix.
183: // * On exit, they are overwritten by the elements of the
184: // * standardised Schur form.
185: // *
186: // * RT1R (output) DOUBLE PRECISION
187: // * RT1I (output) DOUBLE PRECISION
188: // * RT2R (output) DOUBLE PRECISION
189: // * RT2I (output) DOUBLE PRECISION
190: // * The real and imaginary parts of the eigenvalues. If the
191: // * eigenvalues are a complex conjugate pair, RT1I > 0.
192: // *
193: // * CS (output) DOUBLE PRECISION
194: // * SN (output) DOUBLE PRECISION
195: // * Parameters of the rotation matrix.
196: // *
197: // * Further Details
198: // * ===============
199: // *
200: // * Modified by V. Sima, Research Institute for Informatics, Bucharest,
201: // * Romania, to reduce the risk of cancellation errors,
202: // * when computing real eigenvalues, and to ensure, if possible, that
203: // * abs(RT1R) >= abs(RT2R).
204: // *
205: // * =====================================================================
206: // *
207: // * .. Parameters ..
208: // * ..
209: // * .. Local Scalars ..
210: // * ..
211: // * .. External Functions ..
212: // * ..
213: // * .. Intrinsic Functions ..
214: // INTRINSIC ABS, MAX, MIN, SIGN, SQRT;
215: // * ..
216: // * .. Executable Statements ..
217: // *
218:
219: #endregion
220:
221:
222: #region Body
223:
224: EPS = this._dlamch.Run("P");
225: if (C == ZERO)
226: {
227: CS = ONE;
228: SN = ZERO;
229: goto LABEL10;
230: // *
231: }
232: else
233: {
234: if (B == ZERO)
235: {
236: // *
237: // * Swap rows and columns
238: // *
239: CS = ZERO;
240: SN = ONE;
241: TEMP = D;
242: D = A;
243: A = TEMP;
244: B = - C;
245: C = ZERO;
246: goto LABEL10;
247: }
248: else
249: {
250: if ((A - D) == ZERO && FortranLib.Sign(ONE,B) != FortranLib.Sign(ONE,C))
251: {
252: CS = ONE;
253: SN = ZERO;
254: goto LABEL10;
255: }
256: else
257: {
258: // *
259: TEMP = A - D;
260: P = HALF * TEMP;
261: BCMAX = Math.Max(Math.Abs(B), Math.Abs(C));
262: BCMIS = Math.Min(Math.Abs(B), Math.Abs(C)) * FortranLib.Sign(ONE,B) * FortranLib.Sign(ONE,C);
263: SCALE = Math.Max(Math.Abs(P), BCMAX);
264: Z = (P / SCALE) * P + (BCMAX / SCALE) * BCMIS;
265: // *
266: // * If Z is of the order of the machine accuracy, postpone the
267: // * decision on the nature of eigenvalues
268: // *
269: if (Z >= MULTPL * EPS)
270: {
271: // *
272: // * Real eigenvalues. Compute A and D.
273: // *
274: Z = P + FortranLib.Sign(Math.Sqrt(SCALE) * Math.Sqrt(Z),P);
275: A = D + Z;
276: D = D - (BCMAX / Z) * BCMIS;
277: // *
278: // * Compute B and the rotation matrix
279: // *
280: TAU = this._dlapy2.Run(C, Z);
281: CS = Z / TAU;
282: SN = C / TAU;
283: B = B - C;
284: C = ZERO;
285: }
286: else
287: {
288: // *
289: // * Complex eigenvalues, or real (almost) equal eigenvalues.
290: // * Make diagonal elements equal.
291: // *
292: SIGMA = B + C;
293: TAU = this._dlapy2.Run(SIGMA, TEMP);
294: CS = Math.Sqrt(HALF * (ONE + Math.Abs(SIGMA) / TAU));
295: SN = - (P / (TAU * CS)) * FortranLib.Sign(ONE,SIGMA);
296: // *
297: // * Compute [ AA BB ] = [ A B ] [ CS -SN ]
298: // * [ CC DD ] [ C D ] [ SN CS ]
299: // *
300: AA = A * CS + B * SN;
301: BB = - A * SN + B * CS;
302: CC = C * CS + D * SN;
303: DD = - C * SN + D * CS;
304: // *
305: // * Compute [ A B ] = [ CS SN ] [ AA BB ]
306: // * [ C D ] [-SN CS ] [ CC DD ]
307: // *
308: A = AA * CS + CC * SN;
309: B = BB * CS + DD * SN;
310: C = - AA * SN + CC * CS;
311: D = - BB * SN + DD * CS;
312: // *
313: TEMP = HALF * (A + D);
314: A = TEMP;
315: D = TEMP;
316: // *
317: if (C != ZERO)
318: {
319: if (B != ZERO)
320: {
321: if (FortranLib.Sign(ONE,B) == FortranLib.Sign(ONE,C))
322: {
323: // *
324: // * Real eigenvalues: reduce to upper triangular form
325: // *
326: SAB = Math.Sqrt(Math.Abs(B));
327: SAC = Math.Sqrt(Math.Abs(C));
328: P = FortranLib.Sign(SAB * SAC,C);
329: TAU = ONE / Math.Sqrt(Math.Abs(B + C));
330: A = TEMP + P;
331: D = TEMP - P;
332: B = B - C;
333: C = ZERO;
334: CS1 = SAB * TAU;
335: SN1 = SAC * TAU;
336: TEMP = CS * CS1 - SN * SN1;
337: SN = CS * SN1 + SN * CS1;
338: CS = TEMP;
339: }
340: }
341: else
342: {
343: B = - C;
344: C = ZERO;
345: TEMP = CS;
346: CS = - SN;
347: SN = TEMP;
348: }
349: }
350: }
351: // *
352: }
353: }
354: }
355: // *
356: LABEL10:;
357: // *
358: // * Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
359: // *
360: RT1R = A;
361: RT2R = D;
362: if (C == ZERO)
363: {
364: RT1I = ZERO;
365: RT2I = ZERO;
366: }
367: else
368: {
369: RT1I = Math.Sqrt(Math.Abs(B)) * Math.Sqrt(Math.Abs(C));
370: RT2I = - RT1I;
371: }
372: return;
373: // *
374: // * End of DLANV2
375: // *
376:
377: #endregion
378:
379: }
380: }
381: }