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: /// DSBEV computes all the eigenvalues and, optionally, eigenvectors of
27: /// a real symmetric band matrix A.
28: ///
29: ///</summary>
30: public class DSBEV
31: {
32:
33:
34: #region Dependencies
35:
36: LSAME _lsame; DLAMCH _dlamch; DLANSB _dlansb; DLASCL _dlascl; DSBTRD _dsbtrd; DSCAL _dscal; DSTEQR _dsteqr;
37: DSTERF _dsterf;XERBLA _xerbla;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LOWER = false; bool WANTZ = false; int IINFO = 0; int IMAX = 0;
45: int INDE = 0;int INDWRK = 0; int ISCALE = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0; double RMAX = 0;
46: double RMIN = 0;double SAFMIN = 0; double SIGMA = 0; double SMLNUM = 0;
47:
48: #endregion
49:
50: public DSBEV(LSAME lsame, DLAMCH dlamch, DLANSB dlansb, DLASCL dlascl, DSBTRD dsbtrd, DSCAL dscal, DSTEQR dsteqr, DSTERF dsterf, XERBLA xerbla)
51: {
52:
53:
54: #region Set Dependencies
55:
56: this._lsame = lsame; this._dlamch = dlamch; this._dlansb = dlansb; this._dlascl = dlascl; this._dsbtrd = dsbtrd;
57: this._dscal = dscal;this._dsteqr = dsteqr; this._dsterf = dsterf; this._xerbla = xerbla;
58:
59: #endregion
60:
61: }
62:
63: public DSBEV()
64: {
65:
66:
67: #region Dependencies (Initialization)
68:
69: LSAME lsame = new LSAME();
70: DLAMC3 dlamc3 = new DLAMC3();
71: DLASSQ dlassq = new DLASSQ();
72: XERBLA xerbla = new XERBLA();
73: DLAR2V dlar2v = new DLAR2V();
74: DLARGV dlargv = new DLARGV();
75: DLARTV dlartv = new DLARTV();
76: DROT drot = new DROT();
77: DSCAL dscal = new DSCAL();
78: DLAPY2 dlapy2 = new DLAPY2();
79: DLAE2 dlae2 = new DLAE2();
80: DLAEV2 dlaev2 = new DLAEV2();
81: DSWAP dswap = new DSWAP();
82: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
83: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
84: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
85: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
86: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
87: DLANSB dlansb = new DLANSB(dlassq, lsame);
88: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
89: DLARTG dlartg = new DLARTG(dlamch);
90: DLASET dlaset = new DLASET(lsame);
91: DSBTRD dsbtrd = new DSBTRD(dlar2v, dlargv, dlartg, dlartv, dlaset, drot, xerbla, lsame);
92: DLANST dlanst = new DLANST(lsame, dlassq);
93: DLASR dlasr = new DLASR(lsame, xerbla);
94: DLASRT dlasrt = new DLASRT(lsame, xerbla);
95: DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
96: , dlasrt, dswap, xerbla);
97: DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
98:
99: #endregion
100:
101:
102: #region Set Dependencies
103:
104: this._lsame = lsame; this._dlamch = dlamch; this._dlansb = dlansb; this._dlascl = dlascl; this._dsbtrd = dsbtrd;
105: this._dscal = dscal;this._dsteqr = dsteqr; this._dsterf = dsterf; this._xerbla = xerbla;
106:
107: #endregion
108:
109: }
110: /// <summary>
111: /// Purpose
112: /// =======
113: ///
114: /// DSBEV computes all the eigenvalues and, optionally, eigenvectors of
115: /// a real symmetric band matrix A.
116: ///
117: ///</summary>
118: /// <param name="JOBZ">
119: /// (input) CHARACTER*1
120: /// = 'N': Compute eigenvalues only;
121: /// = 'V': Compute eigenvalues and eigenvectors.
122: ///</param>
123: /// <param name="UPLO">
124: /// (input) CHARACTER*1
125: /// = 'U': Upper triangle of A is stored;
126: /// = 'L': Lower triangle of A is stored.
127: ///</param>
128: /// <param name="N">
129: /// (input) INTEGER
130: /// The order of the matrix A. N .GE. 0.
131: ///</param>
132: /// <param name="KD">
133: /// (input) INTEGER
134: /// The number of superdiagonals of the matrix A if UPLO = 'U',
135: /// or the number of subdiagonals if UPLO = 'L'. KD .GE. 0.
136: ///</param>
137: /// <param name="AB">
138: /// (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
139: /// On entry, the upper or lower triangle of the symmetric band
140: /// matrix A, stored in the first KD+1 rows of the array. The
141: /// j-th column of A is stored in the j-th column of the array AB
142: /// as follows:
143: /// if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd).LE.i.LE.j;
144: /// if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j.LE.i.LE.min(n,j+kd).
145: ///
146: /// On exit, AB is overwritten by values generated during the
147: /// reduction to tridiagonal form. If UPLO = 'U', the first
148: /// superdiagonal and the diagonal of the tridiagonal matrix T
149: /// are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
150: /// the diagonal and first subdiagonal of T are returned in the
151: /// first two rows of AB.
152: ///</param>
153: /// <param name="LDAB">
154: /// (input) INTEGER
155: /// The leading dimension of the array AB. LDAB .GE. KD + 1.
156: ///</param>
157: /// <param name="W">
158: /// (output) DOUBLE PRECISION array, dimension (N)
159: /// If INFO = 0, the eigenvalues in ascending order.
160: ///</param>
161: /// <param name="Z">
162: /// (output) DOUBLE PRECISION array, dimension (LDZ, N)
163: /// If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
164: /// eigenvectors of the matrix A, with the i-th column of Z
165: /// holding the eigenvector associated with W(i).
166: /// If JOBZ = 'N', then Z is not referenced.
167: ///</param>
168: /// <param name="LDZ">
169: /// (input) INTEGER
170: /// The leading dimension of the array Z. LDZ .GE. 1, and if
171: /// JOBZ = 'V', LDZ .GE. max(1,N).
172: ///</param>
173: /// <param name="WORK">
174: /// (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
175: ///</param>
176: /// <param name="INFO">
177: /// (output) INTEGER
178: /// = 0: successful exit
179: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
180: /// .GT. 0: if INFO = i, the algorithm failed to converge; i
181: /// off-diagonal elements of an intermediate tridiagonal
182: /// form did not converge to zero.
183: ///</param>
184: public void Run(string JOBZ, string UPLO, int N, int KD, ref double[] AB, int offset_ab, int LDAB
185: , ref double[] W, int offset_w, ref double[] Z, int offset_z, int LDZ, ref double[] WORK, int offset_work, ref int INFO)
186: {
187:
188: #region Array Index Correction
189:
190: int o_ab = -1 - LDAB + offset_ab; int o_w = -1 + offset_w; int o_z = -1 - LDZ + offset_z;
191: int o_work = -1 + offset_work;
192:
193: #endregion
194:
195:
196: #region Strings
197:
198: JOBZ = JOBZ.Substring(0, 1); UPLO = UPLO.Substring(0, 1);
199:
200: #endregion
201:
202:
203: #region Prolog
204:
205: // *
206: // * -- LAPACK driver routine (version 3.1) --
207: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
208: // * November 2006
209: // *
210: // * .. Scalar Arguments ..
211: // * ..
212: // * .. Array Arguments ..
213: // * ..
214: // *
215: // * Purpose
216: // * =======
217: // *
218: // * DSBEV computes all the eigenvalues and, optionally, eigenvectors of
219: // * a real symmetric band matrix A.
220: // *
221: // * Arguments
222: // * =========
223: // *
224: // * JOBZ (input) CHARACTER*1
225: // * = 'N': Compute eigenvalues only;
226: // * = 'V': Compute eigenvalues and eigenvectors.
227: // *
228: // * UPLO (input) CHARACTER*1
229: // * = 'U': Upper triangle of A is stored;
230: // * = 'L': Lower triangle of A is stored.
231: // *
232: // * N (input) INTEGER
233: // * The order of the matrix A. N >= 0.
234: // *
235: // * KD (input) INTEGER
236: // * The number of superdiagonals of the matrix A if UPLO = 'U',
237: // * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
238: // *
239: // * AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N)
240: // * On entry, the upper or lower triangle of the symmetric band
241: // * matrix A, stored in the first KD+1 rows of the array. The
242: // * j-th column of A is stored in the j-th column of the array AB
243: // * as follows:
244: // * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
245: // * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
246: // *
247: // * On exit, AB is overwritten by values generated during the
248: // * reduction to tridiagonal form. If UPLO = 'U', the first
249: // * superdiagonal and the diagonal of the tridiagonal matrix T
250: // * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
251: // * the diagonal and first subdiagonal of T are returned in the
252: // * first two rows of AB.
253: // *
254: // * LDAB (input) INTEGER
255: // * The leading dimension of the array AB. LDAB >= KD + 1.
256: // *
257: // * W (output) DOUBLE PRECISION array, dimension (N)
258: // * If INFO = 0, the eigenvalues in ascending order.
259: // *
260: // * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
261: // * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
262: // * eigenvectors of the matrix A, with the i-th column of Z
263: // * holding the eigenvector associated with W(i).
264: // * If JOBZ = 'N', then Z is not referenced.
265: // *
266: // * LDZ (input) INTEGER
267: // * The leading dimension of the array Z. LDZ >= 1, and if
268: // * JOBZ = 'V', LDZ >= max(1,N).
269: // *
270: // * WORK (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
271: // *
272: // * INFO (output) INTEGER
273: // * = 0: successful exit
274: // * < 0: if INFO = -i, the i-th argument had an illegal value
275: // * > 0: if INFO = i, the algorithm failed to converge; i
276: // * off-diagonal elements of an intermediate tridiagonal
277: // * form did not converge to zero.
278: // *
279: // * =====================================================================
280: // *
281: // * .. Parameters ..
282: // * ..
283: // * .. Local Scalars ..
284: // * ..
285: // * .. External Functions ..
286: // * ..
287: // * .. External Subroutines ..
288: // * ..
289: // * .. Intrinsic Functions ..
290: // INTRINSIC SQRT;
291: // * ..
292: // * .. Executable Statements ..
293: // *
294: // * Test the input parameters.
295: // *
296:
297: #endregion
298:
299:
300: #region Body
301:
302: WANTZ = this._lsame.Run(JOBZ, "V");
303: LOWER = this._lsame.Run(UPLO, "L");
304: // *
305: INFO = 0;
306: if (!(WANTZ || this._lsame.Run(JOBZ, "N")))
307: {
308: INFO = - 1;
309: }
310: else
311: {
312: if (!(LOWER || this._lsame.Run(UPLO, "U")))
313: {
314: INFO = - 2;
315: }
316: else
317: {
318: if (N < 0)
319: {
320: INFO = - 3;
321: }
322: else
323: {
324: if (KD < 0)
325: {
326: INFO = - 4;
327: }
328: else
329: {
330: if (LDAB < KD + 1)
331: {
332: INFO = - 6;
333: }
334: else
335: {
336: if (LDZ < 1 || (WANTZ && LDZ < N))
337: {
338: INFO = - 9;
339: }
340: }
341: }
342: }
343: }
344: }
345: // *
346: if (INFO != 0)
347: {
348: this._xerbla.Run("DSBEV ", - INFO);
349: return;
350: }
351: // *
352: // * Quick return if possible
353: // *
354: if (N == 0) return;
355: // *
356: if (N == 1)
357: {
358: if (LOWER)
359: {
360: W[1 + o_w] = AB[1+1 * LDAB + o_ab];
361: }
362: else
363: {
364: W[1 + o_w] = AB[KD + 1+1 * LDAB + o_ab];
365: }
366: if (WANTZ) Z[1+1 * LDZ + o_z] = ONE;
367: return;
368: }
369: // *
370: // * Get machine constants.
371: // *
372: SAFMIN = this._dlamch.Run("Safe minimum");
373: EPS = this._dlamch.Run("Precision");
374: SMLNUM = SAFMIN / EPS;
375: BIGNUM = ONE / SMLNUM;
376: RMIN = Math.Sqrt(SMLNUM);
377: RMAX = Math.Sqrt(BIGNUM);
378: // *
379: // * Scale matrix to allowable range, if necessary.
380: // *
381: ANRM = this._dlansb.Run("M", UPLO, N, KD, AB, offset_ab, LDAB, ref WORK, offset_work);
382: ISCALE = 0;
383: if (ANRM > ZERO && ANRM < RMIN)
384: {
385: ISCALE = 1;
386: SIGMA = RMIN / ANRM;
387: }
388: else
389: {
390: if (ANRM > RMAX)
391: {
392: ISCALE = 1;
393: SIGMA = RMAX / ANRM;
394: }
395: }
396: if (ISCALE == 1)
397: {
398: if (LOWER)
399: {
400: this._dlascl.Run("B", KD, KD, ONE, SIGMA, N
401: , N, ref AB, offset_ab, LDAB, ref INFO);
402: }
403: else
404: {
405: this._dlascl.Run("Q", KD, KD, ONE, SIGMA, N
406: , N, ref AB, offset_ab, LDAB, ref INFO);
407: }
408: }
409: // *
410: // * Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
411: // *
412: INDE = 1;
413: INDWRK = INDE + N;
414: this._dsbtrd.Run(JOBZ, UPLO, N, KD, ref AB, offset_ab, LDAB
415: , ref W, offset_w, ref WORK, INDE + o_work, ref Z, offset_z, LDZ, ref WORK, INDWRK + o_work, ref IINFO);
416: // *
417: // * For eigenvalues only, call DSTERF. For eigenvectors, call SSTEQR.
418: // *
419: if (!WANTZ)
420: {
421: this._dsterf.Run(N, ref W, offset_w, ref WORK, INDE + o_work, ref INFO);
422: }
423: else
424: {
425: this._dsteqr.Run(JOBZ, N, ref W, offset_w, ref WORK, INDE + o_work, ref Z, offset_z, LDZ
426: , ref WORK, INDWRK + o_work, ref INFO);
427: }
428: // *
429: // * If matrix was scaled, then rescale eigenvalues appropriately.
430: // *
431: if (ISCALE == 1)
432: {
433: if (INFO == 0)
434: {
435: IMAX = N;
436: }
437: else
438: {
439: IMAX = INFO - 1;
440: }
441: this._dscal.Run(IMAX, ONE / SIGMA, ref W, offset_w, 1);
442: }
443: // *
444: return;
445: // *
446: // * End of DSBEV
447: // *
448:
449: #endregion
450:
451: }
452: }
453: }