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: /// DLAED9 finds the roots of the secular equation, as defined by the
27: /// values in D, Z, and RHO, between KSTART and KSTOP. It makes the
28: /// appropriate calls to DLAED4 and then stores the new matrix of
29: /// eigenvectors for use in calculating the next level of Z vectors.
30: ///
31: ///</summary>
32: public class DLAED9
33: {
34:
35:
36: #region Dependencies
37:
38: DLAMC3 _dlamc3; DNRM2 _dnrm2; DCOPY _dcopy; DLAED4 _dlaed4; XERBLA _xerbla;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: int I = 0; int J = 0; double TEMP = 0;
46:
47: #endregion
48:
49: public DLAED9(DLAMC3 dlamc3, DNRM2 dnrm2, DCOPY dcopy, DLAED4 dlaed4, XERBLA xerbla)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dlaed4 = dlaed4; this._xerbla = xerbla;
56:
57: #endregion
58:
59: }
60:
61: public DLAED9()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: DLAMC3 dlamc3 = new DLAMC3();
68: DNRM2 dnrm2 = new DNRM2();
69: DCOPY dcopy = new DCOPY();
70: LSAME lsame = new LSAME();
71: DLAED5 dlaed5 = new DLAED5();
72: XERBLA xerbla = new XERBLA();
73: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
74: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
75: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
76: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
77: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
78: DLAED6 dlaed6 = new DLAED6(dlamch);
79: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
80:
81: #endregion
82:
83:
84: #region Set Dependencies
85:
86: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dlaed4 = dlaed4; this._xerbla = xerbla;
87:
88: #endregion
89:
90: }
91: /// <summary>
92: /// Purpose
93: /// =======
94: ///
95: /// DLAED9 finds the roots of the secular equation, as defined by the
96: /// values in D, Z, and RHO, between KSTART and KSTOP. It makes the
97: /// appropriate calls to DLAED4 and then stores the new matrix of
98: /// eigenvectors for use in calculating the next level of Z vectors.
99: ///
100: ///</summary>
101: /// <param name="K">
102: /// (input) INTEGER
103: /// The number of terms in the rational function to be solved by
104: /// DLAED4. K .GE. 0.
105: ///</param>
106: /// <param name="KSTART">
107: /// (input) INTEGER
108: ///</param>
109: /// <param name="KSTOP">
110: /// (input) INTEGER
111: /// The updated eigenvalues Lambda(I), KSTART .LE. I .LE. KSTOP
112: /// are to be computed. 1 .LE. KSTART .LE. KSTOP .LE. K.
113: ///</param>
114: /// <param name="N">
115: /// (input) INTEGER
116: /// The number of rows and columns in the Q matrix.
117: /// N .GE. K (delation may result in N .GT. K).
118: ///</param>
119: /// <param name="D">
120: /// (output) DOUBLE PRECISION array, dimension (N)
121: /// D(I) contains the updated eigenvalues
122: /// for KSTART .LE. I .LE. KSTOP.
123: ///</param>
124: /// <param name="Q">
125: /// (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
126: ///</param>
127: /// <param name="LDQ">
128: /// (input) INTEGER
129: /// The leading dimension of the array Q. LDQ .GE. max( 1, N ).
130: ///</param>
131: /// <param name="RHO">
132: /// (input) DOUBLE PRECISION
133: /// The value of the parameter in the rank one update equation.
134: /// RHO .GE. 0 required.
135: ///</param>
136: /// <param name="DLAMDA">
137: /// (input) DOUBLE PRECISION array, dimension (K)
138: /// The first K elements of this array contain the old roots
139: /// of the deflated updating problem. These are the poles
140: /// of the secular equation.
141: ///</param>
142: /// <param name="W">
143: /// (input) DOUBLE PRECISION array, dimension (K)
144: /// The first K elements of this array contain the components
145: /// of the deflation-adjusted updating vector.
146: ///</param>
147: /// <param name="S">
148: /// (output) DOUBLE PRECISION array, dimension (LDS, K)
149: /// Will contain the eigenvectors of the repaired matrix which
150: /// will be stored for subsequent Z vector calculation and
151: /// multiplied by the previously accumulated eigenvectors
152: /// to update the system.
153: ///</param>
154: /// <param name="LDS">
155: /// (input) INTEGER
156: /// The leading dimension of S. LDS .GE. max( 1, K ).
157: ///</param>
158: /// <param name="INFO">
159: /// (output) INTEGER
160: /// = 0: successful exit.
161: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
162: /// .GT. 0: if INFO = 1, an eigenvalue did not converge
163: ///</param>
164: public void Run(int K, int KSTART, int KSTOP, int N, ref double[] D, int offset_d, ref double[] Q, int offset_q
165: , int LDQ, double RHO, ref double[] DLAMDA, int offset_dlamda, ref double[] W, int offset_w, ref double[] S, int offset_s, int LDS
166: , ref int INFO)
167: {
168:
169: #region Array Index Correction
170:
171: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_dlamda = -1 + offset_dlamda;
172: int o_w = -1 + offset_w; int o_s = -1 - LDS + offset_s;
173:
174: #endregion
175:
176:
177: #region Prolog
178:
179: // *
180: // * -- LAPACK routine (version 3.1) --
181: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
182: // * November 2006
183: // *
184: // * .. Scalar Arguments ..
185: // * ..
186: // * .. Array Arguments ..
187: // * ..
188: // *
189: // * Purpose
190: // * =======
191: // *
192: // * DLAED9 finds the roots of the secular equation, as defined by the
193: // * values in D, Z, and RHO, between KSTART and KSTOP. It makes the
194: // * appropriate calls to DLAED4 and then stores the new matrix of
195: // * eigenvectors for use in calculating the next level of Z vectors.
196: // *
197: // * Arguments
198: // * =========
199: // *
200: // * K (input) INTEGER
201: // * The number of terms in the rational function to be solved by
202: // * DLAED4. K >= 0.
203: // *
204: // * KSTART (input) INTEGER
205: // * KSTOP (input) INTEGER
206: // * The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
207: // * are to be computed. 1 <= KSTART <= KSTOP <= K.
208: // *
209: // * N (input) INTEGER
210: // * The number of rows and columns in the Q matrix.
211: // * N >= K (delation may result in N > K).
212: // *
213: // * D (output) DOUBLE PRECISION array, dimension (N)
214: // * D(I) contains the updated eigenvalues
215: // * for KSTART <= I <= KSTOP.
216: // *
217: // * Q (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
218: // *
219: // * LDQ (input) INTEGER
220: // * The leading dimension of the array Q. LDQ >= max( 1, N ).
221: // *
222: // * RHO (input) DOUBLE PRECISION
223: // * The value of the parameter in the rank one update equation.
224: // * RHO >= 0 required.
225: // *
226: // * DLAMDA (input) DOUBLE PRECISION array, dimension (K)
227: // * The first K elements of this array contain the old roots
228: // * of the deflated updating problem. These are the poles
229: // * of the secular equation.
230: // *
231: // * W (input) DOUBLE PRECISION array, dimension (K)
232: // * The first K elements of this array contain the components
233: // * of the deflation-adjusted updating vector.
234: // *
235: // * S (output) DOUBLE PRECISION array, dimension (LDS, K)
236: // * Will contain the eigenvectors of the repaired matrix which
237: // * will be stored for subsequent Z vector calculation and
238: // * multiplied by the previously accumulated eigenvectors
239: // * to update the system.
240: // *
241: // * LDS (input) INTEGER
242: // * The leading dimension of S. LDS >= max( 1, K ).
243: // *
244: // * INFO (output) INTEGER
245: // * = 0: successful exit.
246: // * < 0: if INFO = -i, the i-th argument had an illegal value.
247: // * > 0: if INFO = 1, an eigenvalue did not converge
248: // *
249: // * Further Details
250: // * ===============
251: // *
252: // * Based on contributions by
253: // * Jeff Rutter, Computer Science Division, University of California
254: // * at Berkeley, USA
255: // *
256: // * =====================================================================
257: // *
258: // * .. Local Scalars ..
259: // * ..
260: // * .. External Functions ..
261: // * ..
262: // * .. External Subroutines ..
263: // * ..
264: // * .. Intrinsic Functions ..
265: // INTRINSIC MAX, SIGN, SQRT;
266: // * ..
267: // * .. Executable Statements ..
268: // *
269: // * Test the input parameters.
270: // *
271:
272: #endregion
273:
274:
275: #region Body
276:
277: INFO = 0;
278: // *
279: if (K < 0)
280: {
281: INFO = - 1;
282: }
283: else
284: {
285: if (KSTART < 1 || KSTART > Math.Max(1, K))
286: {
287: INFO = - 2;
288: }
289: else
290: {
291: if (Math.Max(1, KSTOP) < KSTART || KSTOP > Math.Max(1, K))
292: {
293: INFO = - 3;
294: }
295: else
296: {
297: if (N < K)
298: {
299: INFO = - 4;
300: }
301: else
302: {
303: if (LDQ < Math.Max(1, K))
304: {
305: INFO = - 7;
306: }
307: else
308: {
309: if (LDS < Math.Max(1, K))
310: {
311: INFO = - 12;
312: }
313: }
314: }
315: }
316: }
317: }
318: if (INFO != 0)
319: {
320: this._xerbla.Run("DLAED9", - INFO);
321: return;
322: }
323: // *
324: // * Quick return if possible
325: // *
326: if (K == 0) return;
327: // *
328: // * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
329: // * be computed with high relative accuracy (barring over/underflow).
330: // * This is a problem on machines without a guard digit in
331: // * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
332: // * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
333: // * which on any of these machines zeros out the bottommost
334: // * bit of DLAMDA(I) if it is 1; this makes the subsequent
335: // * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
336: // * occurs. On binary machines with a guard digit (almost all
337: // * machines) it does not change DLAMDA(I) at all. On hexadecimal
338: // * and decimal machines with a guard digit, it slightly
339: // * changes the bottommost bits of DLAMDA(I). It does not account
340: // * for hexadecimal or decimal machines without guard digits
341: // * (we know of none). We use a subroutine call to compute
342: // * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
343: // * this code.
344: // *
345: for (I = 1; I <= N; I++)
346: {
347: DLAMDA[I + o_dlamda] = this._dlamc3.Run(DLAMDA[I + o_dlamda], DLAMDA[I + o_dlamda]) - DLAMDA[I + o_dlamda];
348: }
349: // *
350: for (J = KSTART; J <= KSTOP; J++)
351: {
352: this._dlaed4.Run(K, J, DLAMDA, offset_dlamda, W, offset_w, ref Q, 1+J * LDQ + o_q, RHO
353: , ref D[J + o_d], ref INFO);
354: // *
355: // * If the zero finder fails, the computation is terminated.
356: // *
357: if (INFO != 0) goto LABEL120;
358: }
359: // *
360: if (K == 1 || K == 2)
361: {
362: for (I = 1; I <= K; I++)
363: {
364: for (J = 1; J <= K; J++)
365: {
366: S[J+I * LDS + o_s] = Q[J+I * LDQ + o_q];
367: }
368: }
369: goto LABEL120;
370: }
371: // *
372: // * Compute updated W.
373: // *
374: this._dcopy.Run(K, W, offset_w, 1, ref S, offset_s, 1);
375: // *
376: // * Initialize W(I) = Q(I,I)
377: // *
378: this._dcopy.Run(K, Q, offset_q, LDQ + 1, ref W, offset_w, 1);
379: for (J = 1; J <= K; J++)
380: {
381: for (I = 1; I <= J - 1; I++)
382: {
383: W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
384: }
385: for (I = J + 1; I <= K; I++)
386: {
387: W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
388: }
389: }
390: for (I = 1; I <= K; I++)
391: {
392: W[I + o_w] = FortranLib.Sign(Math.Sqrt( - W[I + o_w]),S[I+1 * LDS + o_s]);
393: }
394: // *
395: // * Compute eigenvectors of the modified rank-1 modification.
396: // *
397: for (J = 1; J <= K; J++)
398: {
399: for (I = 1; I <= K; I++)
400: {
401: Q[I+J * LDQ + o_q] = W[I + o_w] / Q[I+J * LDQ + o_q];
402: }
403: TEMP = this._dnrm2.Run(K, Q, 1+J * LDQ + o_q, 1);
404: for (I = 1; I <= K; I++)
405: {
406: S[I+J * LDS + o_s] = Q[I+J * LDQ + o_q] / TEMP;
407: }
408: }
409: // *
410: LABEL120:;
411: return;
412: // *
413: // * End of DLAED9
414: // *
415:
416: #endregion
417:
418: }
419: }
420: }