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: /// DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
27: /// using the Pal-Walker-Kahan variant of the QL or QR algorithm.
28: ///
29: ///</summary>
30: public class DSTERF
31: {
32:
33:
34: #region Dependencies
35:
36: DLAMCH _dlamch; DLANST _dlanst; DLAPY2 _dlapy2; DLAE2 _dlae2; DLASCL _dlascl; DLASRT _dlasrt; XERBLA _xerbla;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double THREE = 3.0E0;
44: const int MAXIT = 30;int I = 0; int ISCALE = 0; int JTOT = 0; int L = 0; int L1 = 0; int LEND = 0; int LENDSV = 0;
45: int LSV = 0;int M = 0; int NMAXIT = 0; double ALPHA = 0; double ANORM = 0; double BB = 0; double C = 0; double EPS = 0;
46: double EPS2 = 0;double GAMMA = 0; double OLDC = 0; double OLDGAM = 0; double P = 0; double R = 0; double RT1 = 0;
47: double RT2 = 0;double RTE = 0; double S = 0; double SAFMAX = 0; double SAFMIN = 0; double SIGMA = 0; double SSFMAX = 0;
48: double SSFMIN = 0;
49:
50: #endregion
51:
52: public DSTERF(DLAMCH dlamch, DLANST dlanst, DLAPY2 dlapy2, DLAE2 dlae2, DLASCL dlascl, DLASRT dlasrt, XERBLA xerbla)
53: {
54:
55:
56: #region Set Dependencies
57:
58: this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2; this._dlascl = dlascl;
59: this._dlasrt = dlasrt;this._xerbla = xerbla;
60:
61: #endregion
62:
63: }
64:
65: public DSTERF()
66: {
67:
68:
69: #region Dependencies (Initialization)
70:
71: LSAME lsame = new LSAME();
72: DLAMC3 dlamc3 = new DLAMC3();
73: DLASSQ dlassq = new DLASSQ();
74: DLAPY2 dlapy2 = new DLAPY2();
75: DLAE2 dlae2 = new DLAE2();
76: XERBLA xerbla = new XERBLA();
77: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
78: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
79: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
80: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
81: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
82: DLANST dlanst = new DLANST(lsame, dlassq);
83: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
84: DLASRT dlasrt = new DLASRT(lsame, xerbla);
85:
86: #endregion
87:
88:
89: #region Set Dependencies
90:
91: this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2; this._dlascl = dlascl;
92: this._dlasrt = dlasrt;this._xerbla = xerbla;
93:
94: #endregion
95:
96: }
97: /// <summary>
98: /// Purpose
99: /// =======
100: ///
101: /// DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
102: /// using the Pal-Walker-Kahan variant of the QL or QR algorithm.
103: ///
104: ///</summary>
105: /// <param name="N">
106: /// (input) INTEGER
107: /// The order of the matrix. N .GE. 0.
108: ///</param>
109: /// <param name="D">
110: /// (input/output) DOUBLE PRECISION array, dimension (N)
111: /// On entry, the n diagonal elements of the tridiagonal matrix.
112: /// On exit, if INFO = 0, the eigenvalues in ascending order.
113: ///</param>
114: /// <param name="E">
115: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
116: /// On entry, the (n-1) subdiagonal elements of the tridiagonal
117: /// matrix.
118: /// On exit, E has been destroyed.
119: ///</param>
120: /// <param name="INFO">
121: /// (output) INTEGER
122: /// = 0: successful exit
123: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
124: /// .GT. 0: the algorithm failed to find all of the eigenvalues in
125: /// a total of 30*N iterations; if INFO = i, then i
126: /// elements of E have not converged to zero.
127: ///</param>
128: public void Run(int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref int INFO)
129: {
130:
131: #region Array Index Correction
132:
133: int o_d = -1 + offset_d; int o_e = -1 + offset_e;
134:
135: #endregion
136:
137:
138: #region Prolog
139:
140: // *
141: // * -- LAPACK routine (version 3.1) --
142: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
143: // * November 2006
144: // *
145: // * .. Scalar Arguments ..
146: // * ..
147: // * .. Array Arguments ..
148: // * ..
149: // *
150: // * Purpose
151: // * =======
152: // *
153: // * DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
154: // * using the Pal-Walker-Kahan variant of the QL or QR algorithm.
155: // *
156: // * Arguments
157: // * =========
158: // *
159: // * N (input) INTEGER
160: // * The order of the matrix. N >= 0.
161: // *
162: // * D (input/output) DOUBLE PRECISION array, dimension (N)
163: // * On entry, the n diagonal elements of the tridiagonal matrix.
164: // * On exit, if INFO = 0, the eigenvalues in ascending order.
165: // *
166: // * E (input/output) DOUBLE PRECISION array, dimension (N-1)
167: // * On entry, the (n-1) subdiagonal elements of the tridiagonal
168: // * matrix.
169: // * On exit, E has been destroyed.
170: // *
171: // * INFO (output) INTEGER
172: // * = 0: successful exit
173: // * < 0: if INFO = -i, the i-th argument had an illegal value
174: // * > 0: the algorithm failed to find all of the eigenvalues in
175: // * a total of 30*N iterations; if INFO = i, then i
176: // * elements of E have not converged to zero.
177: // *
178: // * =====================================================================
179: // *
180: // * .. Parameters ..
181: // * ..
182: // * .. Local Scalars ..
183: // * ..
184: // * .. External Functions ..
185: // * ..
186: // * .. External Subroutines ..
187: // * ..
188: // * .. Intrinsic Functions ..
189: // INTRINSIC ABS, SIGN, SQRT;
190: // * ..
191: // * .. Executable Statements ..
192: // *
193: // * Test the input parameters.
194: // *
195:
196: #endregion
197:
198:
199: #region Body
200:
201: INFO = 0;
202: // *
203: // * Quick return if possible
204: // *
205: if (N < 0)
206: {
207: INFO = - 1;
208: this._xerbla.Run("DSTERF", - INFO);
209: return;
210: }
211: if (N <= 1) return;
212: // *
213: // * Determine the unit roundoff for this environment.
214: // *
215: EPS = this._dlamch.Run("E");
216: EPS2 = Math.Pow(EPS,2);
217: SAFMIN = this._dlamch.Run("S");
218: SAFMAX = ONE / SAFMIN;
219: SSFMAX = Math.Sqrt(SAFMAX) / THREE;
220: SSFMIN = Math.Sqrt(SAFMIN) / EPS2;
221: // *
222: // * Compute the eigenvalues of the tridiagonal matrix.
223: // *
224: NMAXIT = N * MAXIT;
225: SIGMA = ZERO;
226: JTOT = 0;
227: // *
228: // * Determine where the matrix splits and choose QL or QR iteration
229: // * for each block, according to whether top or bottom diagonal
230: // * element is smaller.
231: // *
232: L1 = 1;
233: // *
234: LABEL10:;
235: if (L1 > N) goto LABEL170;
236: if (L1 > 1) E[L1 - 1 + o_e] = ZERO;
237: for (M = L1; M <= N - 1; M++)
238: {
239: if (Math.Abs(E[M + o_e]) <= (Math.Sqrt(Math.Abs(D[M + o_d])) * Math.Sqrt(Math.Abs(D[M + 1 + o_d]))) * EPS)
240: {
241: E[M + o_e] = ZERO;
242: goto LABEL30;
243: }
244: }
245: M = N;
246: // *
247: LABEL30:;
248: L = L1;
249: LSV = L;
250: LEND = M;
251: LENDSV = LEND;
252: L1 = M + 1;
253: if (LEND == L) goto LABEL10;
254: // *
255: // * Scale submatrix in rows and columns L to LEND
256: // *
257: ANORM = this._dlanst.Run("I", LEND - L + 1, D, L + o_d, E, L + o_e);
258: ISCALE = 0;
259: if (ANORM > SSFMAX)
260: {
261: ISCALE = 1;
262: this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L + 1
263: , 1, ref D, L + o_d, N, ref INFO);
264: this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L
265: , 1, ref E, L + o_e, N, ref INFO);
266: }
267: else
268: {
269: if (ANORM < SSFMIN)
270: {
271: ISCALE = 2;
272: this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L + 1
273: , 1, ref D, L + o_d, N, ref INFO);
274: this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L
275: , 1, ref E, L + o_e, N, ref INFO);
276: }
277: }
278: // *
279: for (I = L; I <= LEND - 1; I++)
280: {
281: E[I + o_e] = Math.Pow(E[I + o_e],2);
282: }
283: // *
284: // * Choose between QL and QR iteration
285: // *
286: if (Math.Abs(D[LEND + o_d]) < Math.Abs(D[L + o_d]))
287: {
288: LEND = LSV;
289: L = LENDSV;
290: }
291: // *
292: if (LEND >= L)
293: {
294: // *
295: // * QL Iteration
296: // *
297: // * Look for small subdiagonal element.
298: // *
299: LABEL50:;
300: if (L != LEND)
301: {
302: for (M = L; M <= LEND - 1; M++)
303: {
304: if (Math.Abs(E[M + o_e]) <= EPS2 * Math.Abs(D[M + o_d] * D[M + 1 + o_d])) goto LABEL70;
305: }
306: }
307: M = LEND;
308: // *
309: LABEL70:;
310: if (M < LEND) E[M + o_e] = ZERO;
311: P = D[L + o_d];
312: if (M == L) goto LABEL90;
313: // *
314: // * If remaining matrix is 2 by 2, use DLAE2 to compute its
315: // * eigenvalues.
316: // *
317: if (M == L + 1)
318: {
319: RTE = Math.Sqrt(E[L + o_e]);
320: this._dlae2.Run(D[L + o_d], RTE, D[L + 1 + o_d], ref RT1, ref RT2);
321: D[L + o_d] = RT1;
322: D[L + 1 + o_d] = RT2;
323: E[L + o_e] = ZERO;
324: L = L + 2;
325: if (L <= LEND) goto LABEL50;
326: goto LABEL150;
327: }
328: // *
329: if (JTOT == NMAXIT) goto LABEL150;
330: JTOT = JTOT + 1;
331: // *
332: // * Form shift.
333: // *
334: RTE = Math.Sqrt(E[L + o_e]);
335: SIGMA = (D[L + 1 + o_d] - P) / (TWO * RTE);
336: R = this._dlapy2.Run(SIGMA, ONE);
337: SIGMA = P - (RTE / (SIGMA + FortranLib.Sign(R,SIGMA)));
338: // *
339: C = ONE;
340: S = ZERO;
341: GAMMA = D[M + o_d] - SIGMA;
342: P = GAMMA * GAMMA;
343: // *
344: // * Inner loop
345: // *
346: for (I = M - 1; I >= L; I += - 1)
347: {
348: BB = E[I + o_e];
349: R = P + BB;
350: if (I != M - 1) E[I + 1 + o_e] = S * R;
351: OLDC = C;
352: C = P / R;
353: S = BB / R;
354: OLDGAM = GAMMA;
355: ALPHA = D[I + o_d];
356: GAMMA = C * (ALPHA - SIGMA) - S * OLDGAM;
357: D[I + 1 + o_d] = OLDGAM + (ALPHA - GAMMA);
358: if (C != ZERO)
359: {
360: P = (GAMMA * GAMMA) / C;
361: }
362: else
363: {
364: P = OLDC * BB;
365: }
366: }
367: // *
368: E[L + o_e] = S * P;
369: D[L + o_d] = SIGMA + GAMMA;
370: goto LABEL50;
371: // *
372: // * Eigenvalue found.
373: // *
374: LABEL90:;
375: D[L + o_d] = P;
376: // *
377: L = L + 1;
378: if (L <= LEND) goto LABEL50;
379: goto LABEL150;
380: // *
381: }
382: else
383: {
384: // *
385: // * QR Iteration
386: // *
387: // * Look for small superdiagonal element.
388: // *
389: LABEL100:;
390: for (M = L; M >= LEND + 1; M += - 1)
391: {
392: if (Math.Abs(E[M - 1 + o_e]) <= EPS2 * Math.Abs(D[M + o_d] * D[M - 1 + o_d])) goto LABEL120;
393: }
394: M = LEND;
395: // *
396: LABEL120:;
397: if (M > LEND) E[M - 1 + o_e] = ZERO;
398: P = D[L + o_d];
399: if (M == L) goto LABEL140;
400: // *
401: // * If remaining matrix is 2 by 2, use DLAE2 to compute its
402: // * eigenvalues.
403: // *
404: if (M == L - 1)
405: {
406: RTE = Math.Sqrt(E[L - 1 + o_e]);
407: this._dlae2.Run(D[L + o_d], RTE, D[L - 1 + o_d], ref RT1, ref RT2);
408: D[L + o_d] = RT1;
409: D[L - 1 + o_d] = RT2;
410: E[L - 1 + o_e] = ZERO;
411: L = L - 2;
412: if (L >= LEND) goto LABEL100;
413: goto LABEL150;
414: }
415: // *
416: if (JTOT == NMAXIT) goto LABEL150;
417: JTOT = JTOT + 1;
418: // *
419: // * Form shift.
420: // *
421: RTE = Math.Sqrt(E[L - 1 + o_e]);
422: SIGMA = (D[L - 1 + o_d] - P) / (TWO * RTE);
423: R = this._dlapy2.Run(SIGMA, ONE);
424: SIGMA = P - (RTE / (SIGMA + FortranLib.Sign(R,SIGMA)));
425: // *
426: C = ONE;
427: S = ZERO;
428: GAMMA = D[M + o_d] - SIGMA;
429: P = GAMMA * GAMMA;
430: // *
431: // * Inner loop
432: // *
433: for (I = M; I <= L - 1; I++)
434: {
435: BB = E[I + o_e];
436: R = P + BB;
437: if (I != M) E[I - 1 + o_e] = S * R;
438: OLDC = C;
439: C = P / R;
440: S = BB / R;
441: OLDGAM = GAMMA;
442: ALPHA = D[I + 1 + o_d];
443: GAMMA = C * (ALPHA - SIGMA) - S * OLDGAM;
444: D[I + o_d] = OLDGAM + (ALPHA - GAMMA);
445: if (C != ZERO)
446: {
447: P = (GAMMA * GAMMA) / C;
448: }
449: else
450: {
451: P = OLDC * BB;
452: }
453: }
454: // *
455: E[L - 1 + o_e] = S * P;
456: D[L + o_d] = SIGMA + GAMMA;
457: goto LABEL100;
458: // *
459: // * Eigenvalue found.
460: // *
461: LABEL140:;
462: D[L + o_d] = P;
463: // *
464: L = L - 1;
465: if (L >= LEND) goto LABEL100;
466: goto LABEL150;
467: // *
468: }
469: // *
470: // * Undo scaling if necessary
471: // *
472: LABEL150:;
473: if (ISCALE == 1)
474: {
475: this._dlascl.Run("G", 0, 0, SSFMAX, ANORM, LENDSV - LSV + 1
476: , 1, ref D, LSV + o_d, N, ref INFO);
477: }
478: if (ISCALE == 2)
479: {
480: this._dlascl.Run("G", 0, 0, SSFMIN, ANORM, LENDSV - LSV + 1
481: , 1, ref D, LSV + o_d, N, ref INFO);
482: }
483: // *
484: // * Check for no convergence to an eigenvalue after a total
485: // * of N*MAXIT iterations.
486: // *
487: if (JTOT < NMAXIT) goto LABEL10;
488: for (I = 1; I <= N - 1; I++)
489: {
490: if (E[I + o_e] != ZERO) INFO = INFO + 1;
491: }
492: goto LABEL180;
493: // *
494: // * Sort eigenvalues in increasing order.
495: // *
496: LABEL170:;
497: this._dlasrt.Run("I", N, ref D, offset_d, ref INFO);
498: // *
499: LABEL180:;
500: return;
501: // *
502: // * End of DSTERF
503: // *
504:
505: #endregion
506:
507: }
508: }
509: }