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