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 routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// This subroutine computes the I-th eigenvalue of a symmetric rank-one
27: /// modification of a 2-by-2 diagonal matrix
28: ///
29: /// diag( D ) + RHO * Z * transpose(Z) .
30: ///
31: /// The diagonal elements in the array D are assumed to satisfy
32: ///
33: /// D(i) .LT. D(j) for i .LT. j .
34: ///
35: /// We also assume RHO .GT. 0 and that the Euclidean norm of the vector
36: /// Z is one.
37: ///
38: ///</summary>
39: public class DLAED5
40: {
41:
42:
43: #region Fields
44:
45: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double FOUR = 4.0E0; double B = 0;
46: double C = 0;double DEL = 0; double TAU = 0; double TEMP = 0; double W = 0;
47:
48: #endregion
49:
50: public DLAED5()
51: {
52:
53: }
54:
55: /// <summary>
56: /// Purpose
57: /// =======
58: ///
59: /// This subroutine computes the I-th eigenvalue of a symmetric rank-one
60: /// modification of a 2-by-2 diagonal matrix
61: ///
62: /// diag( D ) + RHO * Z * transpose(Z) .
63: ///
64: /// The diagonal elements in the array D are assumed to satisfy
65: ///
66: /// D(i) .LT. D(j) for i .LT. j .
67: ///
68: /// We also assume RHO .GT. 0 and that the Euclidean norm of the vector
69: /// Z is one.
70: ///
71: ///</summary>
72: /// <param name="I">
73: /// (input) INTEGER
74: /// The index of the eigenvalue to be computed. I = 1 or I = 2.
75: ///</param>
76: /// <param name="D">
77: /// (input) DOUBLE PRECISION array, dimension (2)
78: /// The original eigenvalues. We assume D(1) .LT. D(2).
79: ///</param>
80: /// <param name="Z">
81: /// (input) DOUBLE PRECISION array, dimension (2)
82: /// The components of the updating vector.
83: ///</param>
84: /// <param name="DELTA">
85: /// (output) DOUBLE PRECISION array, dimension (2)
86: /// The vector DELTA contains the information necessary
87: /// to construct the eigenvectors.
88: ///</param>
89: /// <param name="RHO">
90: /// (input) DOUBLE PRECISION
91: /// The scalar in the symmetric updating formula.
92: ///</param>
93: /// <param name="DLAM">
94: /// (output) DOUBLE PRECISION
95: /// The computed lambda_I, the I-th updated eigenvalue.
96: ///</param>
97: public void Run(int I, double[] D, int offset_d, double[] Z, int offset_z, ref double[] DELTA, int offset_delta, double RHO, ref double DLAM)
98: {
99:
100: #region Array Index Correction
101:
102: int o_d = -1 + offset_d; int o_z = -1 + offset_z; int o_delta = -1 + offset_delta;
103:
104: #endregion
105:
106:
107: #region Prolog
108:
109: // *
110: // * -- LAPACK routine (version 3.1) --
111: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
112: // * November 2006
113: // *
114: // * .. Scalar Arguments ..
115: // * ..
116: // * .. Array Arguments ..
117: // * ..
118: // *
119: // * Purpose
120: // * =======
121: // *
122: // * This subroutine computes the I-th eigenvalue of a symmetric rank-one
123: // * modification of a 2-by-2 diagonal matrix
124: // *
125: // * diag( D ) + RHO * Z * transpose(Z) .
126: // *
127: // * The diagonal elements in the array D are assumed to satisfy
128: // *
129: // * D(i) < D(j) for i < j .
130: // *
131: // * We also assume RHO > 0 and that the Euclidean norm of the vector
132: // * Z is one.
133: // *
134: // * Arguments
135: // * =========
136: // *
137: // * I (input) INTEGER
138: // * The index of the eigenvalue to be computed. I = 1 or I = 2.
139: // *
140: // * D (input) DOUBLE PRECISION array, dimension (2)
141: // * The original eigenvalues. We assume D(1) < D(2).
142: // *
143: // * Z (input) DOUBLE PRECISION array, dimension (2)
144: // * The components of the updating vector.
145: // *
146: // * DELTA (output) DOUBLE PRECISION array, dimension (2)
147: // * The vector DELTA contains the information necessary
148: // * to construct the eigenvectors.
149: // *
150: // * RHO (input) DOUBLE PRECISION
151: // * The scalar in the symmetric updating formula.
152: // *
153: // * DLAM (output) DOUBLE PRECISION
154: // * The computed lambda_I, the I-th updated eigenvalue.
155: // *
156: // * Further Details
157: // * ===============
158: // *
159: // * Based on contributions by
160: // * Ren-Cang Li, Computer Science Division, University of California
161: // * at Berkeley, USA
162: // *
163: // * =====================================================================
164: // *
165: // * .. Parameters ..
166: // * ..
167: // * .. Local Scalars ..
168: // * ..
169: // * .. Intrinsic Functions ..
170: // INTRINSIC ABS, SQRT;
171: // * ..
172: // * .. Executable Statements ..
173: // *
174:
175: #endregion
176:
177:
178: #region Body
179:
180: DEL = D[2 + o_d] - D[1 + o_d];
181: if (I == 1)
182: {
183: W = ONE + TWO * RHO * (Z[2 + o_z] * Z[2 + o_z] - Z[1 + o_z] * Z[1 + o_z]) / DEL;
184: if (W > ZERO)
185: {
186: B = DEL + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);
187: C = RHO * Z[1 + o_z] * Z[1 + o_z] * DEL;
188: // *
189: // * B > ZERO, always
190: // *
191: TAU = TWO * C / (B + Math.Sqrt(Math.Abs(B * B - FOUR * C)));
192: DLAM = D[1 + o_d] + TAU;
193: DELTA[1 + o_delta] = - Z[1 + o_z] / TAU;
194: DELTA[2 + o_delta] = Z[2 + o_z] / (DEL - TAU);
195: }
196: else
197: {
198: B = - DEL + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);
199: C = RHO * Z[2 + o_z] * Z[2 + o_z] * DEL;
200: if (B > ZERO)
201: {
202: TAU = - TWO * C / (B + Math.Sqrt(B * B + FOUR * C));
203: }
204: else
205: {
206: TAU = (B - Math.Sqrt(B * B + FOUR * C)) / TWO;
207: }
208: DLAM = D[2 + o_d] + TAU;
209: DELTA[1 + o_delta] = - Z[1 + o_z] / (DEL + TAU);
210: DELTA[2 + o_delta] = - Z[2 + o_z] / TAU;
211: }
212: TEMP = Math.Sqrt(DELTA[1 + o_delta] * DELTA[1 + o_delta] + DELTA[2 + o_delta] * DELTA[2 + o_delta]);
213: DELTA[1 + o_delta] = DELTA[1 + o_delta] / TEMP;
214: DELTA[2 + o_delta] = DELTA[2 + o_delta] / TEMP;
215: }
216: else
217: {
218: // *
219: // * Now I=2
220: // *
221: B = - DEL + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);
222: C = RHO * Z[2 + o_z] * Z[2 + o_z] * DEL;
223: if (B > ZERO)
224: {
225: TAU = (B + Math.Sqrt(B * B + FOUR * C)) / TWO;
226: }
227: else
228: {
229: TAU = TWO * C / ( - B + Math.Sqrt(B * B + FOUR * C));
230: }
231: DLAM = D[2 + o_d] + TAU;
232: DELTA[1 + o_delta] = - Z[1 + o_z] / (DEL + TAU);
233: DELTA[2 + o_delta] = - Z[2 + o_z] / TAU;
234: TEMP = Math.Sqrt(DELTA[1 + o_delta] * DELTA[1 + o_delta] + DELTA[2 + o_delta] * DELTA[2 + o_delta]);
235: DELTA[1 + o_delta] = DELTA[1 + o_delta] / TEMP;
236: DELTA[2 + o_delta] = DELTA[2 + o_delta] / TEMP;
237: }
238: return;
239: // *
240: // * End OF DLAED5
241: // *
242:
243: #endregion
244:
245: }
246: }
247: }