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: /// DSYEV computes all eigenvalues and, optionally, eigenvectors of a
27: /// real symmetric matrix A.
28: ///
29: ///</summary>
30: public class DSYEV
31: {
32:
33:
34: #region Dependencies
35:
36: LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANSY _dlansy; DLASCL _dlascl; DORGTR _dorgtr; DSCAL _dscal;
37: DSTEQR _dsteqr;DSTERF _dsterf; DSYTRD _dsytrd; 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 LQUERY = false; bool WANTZ = false;
45: int IINFO = 0;int IMAX = 0; int INDE = 0; int INDTAU = 0; int INDWRK = 0; int ISCALE = 0; int LLWORK = 0; int LWKOPT = 0;
46: int NB = 0;double ANRM = 0; double BIGNUM = 0; double EPS = 0; double RMAX = 0; double RMIN = 0; double SAFMIN = 0;
47: double SIGMA = 0;double SMLNUM = 0;
48:
49: #endregion
50:
51: public DSYEV(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANSY dlansy, DLASCL dlascl, DORGTR dorgtr, DSCAL dscal, DSTEQR dsteqr, DSTERF dsterf, DSYTRD dsytrd
52: , XERBLA xerbla)
53: {
54:
55:
56: #region Set Dependencies
57:
58: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlansy = dlansy; this._dlascl = dlascl;
59: this._dorgtr = dorgtr;this._dscal = dscal; this._dsteqr = dsteqr; this._dsterf = dsterf; this._dsytrd = dsytrd;
60: this._xerbla = xerbla;
61:
62: #endregion
63:
64: }
65:
66: public DSYEV()
67: {
68:
69:
70: #region Dependencies (Initialization)
71:
72: LSAME lsame = new LSAME();
73: IEEECK ieeeck = new IEEECK();
74: IPARMQ iparmq = new IPARMQ();
75: DLAMC3 dlamc3 = new DLAMC3();
76: DLASSQ dlassq = new DLASSQ();
77: XERBLA xerbla = new XERBLA();
78: DCOPY dcopy = new DCOPY();
79: DSCAL dscal = new DSCAL();
80: DLAPY2 dlapy2 = new DLAPY2();
81: DLAE2 dlae2 = new DLAE2();
82: DLAEV2 dlaev2 = new DLAEV2();
83: DSWAP dswap = new DSWAP();
84: DAXPY daxpy = new DAXPY();
85: DNRM2 dnrm2 = new DNRM2();
86: DDOT ddot = new DDOT();
87: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
88: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
89: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
90: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
91: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
92: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
93: DLANSY dlansy = new DLANSY(dlassq, lsame);
94: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
95: DGEMM dgemm = new DGEMM(lsame, xerbla);
96: DTRMM dtrmm = new DTRMM(lsame, xerbla);
97: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
98: DGEMV dgemv = new DGEMV(lsame, xerbla);
99: DTRMV dtrmv = new DTRMV(lsame, xerbla);
100: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
101: DGER dger = new DGER(xerbla);
102: DLARF dlarf = new DLARF(dgemv, dger, lsame);
103: DORG2L dorg2l = new DORG2L(dlarf, dscal, xerbla);
104: DORGQL dorgql = new DORGQL(dlarfb, dlarft, dorg2l, xerbla, ilaenv);
105: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
106: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
107: DORGTR dorgtr = new DORGTR(lsame, ilaenv, dorgql, dorgqr, xerbla);
108: DLANST dlanst = new DLANST(lsame, dlassq);
109: DLARTG dlartg = new DLARTG(dlamch);
110: DLASET dlaset = new DLASET(lsame);
111: DLASR dlasr = new DLASR(lsame, xerbla);
112: DLASRT dlasrt = new DLASRT(lsame, xerbla);
113: DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
114: , dlasrt, dswap, xerbla);
115: DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
116: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
117: DSYMV dsymv = new DSYMV(lsame, xerbla);
118: DLATRD dlatrd = new DLATRD(daxpy, dgemv, dlarfg, dscal, dsymv, lsame, ddot);
119: DSYR2K dsyr2k = new DSYR2K(lsame, xerbla);
120: DSYR2 dsyr2 = new DSYR2(lsame, xerbla);
121: DSYTD2 dsytd2 = new DSYTD2(daxpy, dlarfg, dsymv, dsyr2, xerbla, lsame, ddot);
122: DSYTRD dsytrd = new DSYTRD(dlatrd, dsyr2k, dsytd2, xerbla, lsame, ilaenv);
123:
124: #endregion
125:
126:
127: #region Set Dependencies
128:
129: this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlansy = dlansy; this._dlascl = dlascl;
130: this._dorgtr = dorgtr;this._dscal = dscal; this._dsteqr = dsteqr; this._dsterf = dsterf; this._dsytrd = dsytrd;
131: this._xerbla = xerbla;
132:
133: #endregion
134:
135: }
136: /// <summary>
137: /// Purpose
138: /// =======
139: ///
140: /// DSYEV computes all eigenvalues and, optionally, eigenvectors of a
141: /// real symmetric matrix A.
142: ///
143: ///</summary>
144: /// <param name="JOBZ">
145: /// (input) CHARACTER*1
146: /// = 'N': Compute eigenvalues only;
147: /// = 'V': Compute eigenvalues and eigenvectors.
148: ///</param>
149: /// <param name="UPLO">
150: /// (input) CHARACTER*1
151: /// = 'U': Upper triangle of A is stored;
152: /// = 'L': Lower triangle of A is stored.
153: ///</param>
154: /// <param name="N">
155: /// (input) INTEGER
156: /// The order of the matrix A. N .GE. 0.
157: ///</param>
158: /// <param name="A">
159: /// (input/output) DOUBLE PRECISION array, dimension (LDA, N)
160: /// On entry, the symmetric matrix A. If UPLO = 'U', the
161: /// leading N-by-N upper triangular part of A contains the
162: /// upper triangular part of the matrix A. If UPLO = 'L',
163: /// the leading N-by-N lower triangular part of A contains
164: /// the lower triangular part of the matrix A.
165: /// On exit, if JOBZ = 'V', then if INFO = 0, A contains the
166: /// orthonormal eigenvectors of the matrix A.
167: /// If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
168: /// or the upper triangle (if UPLO='U') of A, including the
169: /// diagonal, is destroyed.
170: ///</param>
171: /// <param name="LDA">
172: /// (input) INTEGER
173: /// The leading dimension of the array A. LDA .GE. max(1,N).
174: ///</param>
175: /// <param name="W">
176: /// (output) DOUBLE PRECISION array, dimension (N)
177: /// If INFO = 0, the eigenvalues in ascending order.
178: ///</param>
179: /// <param name="WORK">
180: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182: ///</param>
183: /// <param name="LWORK">
184: /// (input) INTEGER
185: /// The length of the array WORK. LWORK .GE. max(1,3*N-1).
186: /// For optimal efficiency, LWORK .GE. (NB+2)*N,
187: /// where NB is the blocksize for DSYTRD returned by ILAENV.
188: ///
189: /// If LWORK = -1, then a workspace query is assumed; the routine
190: /// only calculates the optimal size of the WORK array, returns
191: /// this value as the first entry of the WORK array, and no error
192: /// message related to LWORK is issued by XERBLA.
193: ///</param>
194: /// <param name="INFO">
195: /// (output) INTEGER
196: /// = 0: successful exit
197: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
198: /// .GT. 0: if INFO = i, the algorithm failed to converge; i
199: /// off-diagonal elements of an intermediate tridiagonal
200: /// form did not converge to zero.
201: ///</param>
202: public void Run(string JOBZ, string UPLO, int N, ref double[] A, int offset_a, int LDA, ref double[] W, int offset_w
203: , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
204: {
205:
206: #region Array Index Correction
207:
208: int o_a = -1 - LDA + offset_a; int o_w = -1 + offset_w; int o_work = -1 + offset_work;
209:
210: #endregion
211:
212:
213: #region Strings
214:
215: JOBZ = JOBZ.Substring(0, 1); UPLO = UPLO.Substring(0, 1);
216:
217: #endregion
218:
219:
220: #region Prolog
221:
222: // *
223: // * -- LAPACK driver routine (version 3.1) --
224: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
225: // * November 2006
226: // *
227: // * .. Scalar Arguments ..
228: // * ..
229: // * .. Array Arguments ..
230: // * ..
231: // *
232: // * Purpose
233: // * =======
234: // *
235: // * DSYEV computes all eigenvalues and, optionally, eigenvectors of a
236: // * real symmetric matrix A.
237: // *
238: // * Arguments
239: // * =========
240: // *
241: // * JOBZ (input) CHARACTER*1
242: // * = 'N': Compute eigenvalues only;
243: // * = 'V': Compute eigenvalues and eigenvectors.
244: // *
245: // * UPLO (input) CHARACTER*1
246: // * = 'U': Upper triangle of A is stored;
247: // * = 'L': Lower triangle of A is stored.
248: // *
249: // * N (input) INTEGER
250: // * The order of the matrix A. N >= 0.
251: // *
252: // * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
253: // * On entry, the symmetric matrix A. If UPLO = 'U', the
254: // * leading N-by-N upper triangular part of A contains the
255: // * upper triangular part of the matrix A. If UPLO = 'L',
256: // * the leading N-by-N lower triangular part of A contains
257: // * the lower triangular part of the matrix A.
258: // * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
259: // * orthonormal eigenvectors of the matrix A.
260: // * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
261: // * or the upper triangle (if UPLO='U') of A, including the
262: // * diagonal, is destroyed.
263: // *
264: // * LDA (input) INTEGER
265: // * The leading dimension of the array A. LDA >= max(1,N).
266: // *
267: // * W (output) DOUBLE PRECISION array, dimension (N)
268: // * If INFO = 0, the eigenvalues in ascending order.
269: // *
270: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
271: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
272: // *
273: // * LWORK (input) INTEGER
274: // * The length of the array WORK. LWORK >= max(1,3*N-1).
275: // * For optimal efficiency, LWORK >= (NB+2)*N,
276: // * where NB is the blocksize for DSYTRD returned by ILAENV.
277: // *
278: // * If LWORK = -1, then a workspace query is assumed; the routine
279: // * only calculates the optimal size of the WORK array, returns
280: // * this value as the first entry of the WORK array, and no error
281: // * message related to LWORK is issued by XERBLA.
282: // *
283: // * INFO (output) INTEGER
284: // * = 0: successful exit
285: // * < 0: if INFO = -i, the i-th argument had an illegal value
286: // * > 0: if INFO = i, the algorithm failed to converge; i
287: // * off-diagonal elements of an intermediate tridiagonal
288: // * form did not converge to zero.
289: // *
290: // * =====================================================================
291: // *
292: // * .. Parameters ..
293: // * ..
294: // * .. Local Scalars ..
295: // * ..
296: // * .. External Functions ..
297: // * ..
298: // * .. External Subroutines ..
299: // * ..
300: // * .. Intrinsic Functions ..
301: // INTRINSIC MAX, SQRT;
302: // * ..
303: // * .. Executable Statements ..
304: // *
305: // * Test the input parameters.
306: // *
307:
308: #endregion
309:
310:
311: #region Body
312:
313: WANTZ = this._lsame.Run(JOBZ, "V");
314: LOWER = this._lsame.Run(UPLO, "L");
315: LQUERY = (LWORK == - 1);
316: // *
317: INFO = 0;
318: if (!(WANTZ || this._lsame.Run(JOBZ, "N")))
319: {
320: INFO = - 1;
321: }
322: else
323: {
324: if (!(LOWER || this._lsame.Run(UPLO, "U")))
325: {
326: INFO = - 2;
327: }
328: else
329: {
330: if (N < 0)
331: {
332: INFO = - 3;
333: }
334: else
335: {
336: if (LDA < Math.Max(1, N))
337: {
338: INFO = - 5;
339: }
340: }
341: }
342: }
343: // *
344: if (INFO == 0)
345: {
346: NB = this._ilaenv.Run(1, "DSYTRD", UPLO, N, - 1, - 1, - 1);
347: LWKOPT = Math.Max(1, (NB + 2) * N);
348: WORK[1 + o_work] = LWKOPT;
349: // *
350: if (LWORK < Math.Max(1, 3 * N - 1) && !LQUERY) INFO = - 8;
351: }
352: // *
353: if (INFO != 0)
354: {
355: this._xerbla.Run("DSYEV ", - INFO);
356: return;
357: }
358: else
359: {
360: if (LQUERY)
361: {
362: return;
363: }
364: }
365: // *
366: // * Quick return if possible
367: // *
368: if (N == 0)
369: {
370: return;
371: }
372: // *
373: if (N == 1)
374: {
375: W[1 + o_w] = A[1+1 * LDA + o_a];
376: WORK[1 + o_work] = 2;
377: if (WANTZ) A[1+1 * LDA + o_a] = ONE;
378: return;
379: }
380: // *
381: // * Get machine constants.
382: // *
383: SAFMIN = this._dlamch.Run("Safe minimum");
384: EPS = this._dlamch.Run("Precision");
385: SMLNUM = SAFMIN / EPS;
386: BIGNUM = ONE / SMLNUM;
387: RMIN = Math.Sqrt(SMLNUM);
388: RMAX = Math.Sqrt(BIGNUM);
389: // *
390: // * Scale matrix to allowable range, if necessary.
391: // *
392: ANRM = this._dlansy.Run("M", UPLO, N, A, offset_a, LDA, ref WORK, offset_work);
393: ISCALE = 0;
394: if (ANRM > ZERO && ANRM < RMIN)
395: {
396: ISCALE = 1;
397: SIGMA = RMIN / ANRM;
398: }
399: else
400: {
401: if (ANRM > RMAX)
402: {
403: ISCALE = 1;
404: SIGMA = RMAX / ANRM;
405: }
406: }
407: if (ISCALE == 1)
408: {
409: this._dlascl.Run(UPLO, 0, 0, ONE, SIGMA, N
410: , N, ref A, offset_a, LDA, ref INFO);
411: }
412: // *
413: // * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
414: // *
415: INDE = 1;
416: INDTAU = INDE + N;
417: INDWRK = INDTAU + N;
418: LLWORK = LWORK - INDWRK + 1;
419: this._dsytrd.Run(UPLO, N, ref A, offset_a, LDA, ref W, offset_w, ref WORK, INDE + o_work
420: , ref WORK, INDTAU + o_work, ref WORK, INDWRK + o_work, LLWORK, ref IINFO);
421: // *
422: // * For eigenvalues only, call DSTERF. For eigenvectors, first call
423: // * DORGTR to generate the orthogonal matrix, then call DSTEQR.
424: // *
425: if (!WANTZ)
426: {
427: this._dsterf.Run(N, ref W, offset_w, ref WORK, INDE + o_work, ref INFO);
428: }
429: else
430: {
431: this._dorgtr.Run(UPLO, N, ref A, offset_a, LDA, WORK, INDTAU + o_work, ref WORK, INDWRK + o_work
432: , LLWORK, ref IINFO);
433: this._dsteqr.Run(JOBZ, N, ref W, offset_w, ref WORK, INDE + o_work, ref A, offset_a, LDA
434: , ref WORK, INDTAU + o_work, ref INFO);
435: }
436: // *
437: // * If matrix was scaled, then rescale eigenvalues appropriately.
438: // *
439: if (ISCALE == 1)
440: {
441: if (INFO == 0)
442: {
443: IMAX = N;
444: }
445: else
446: {
447: IMAX = INFO - 1;
448: }
449: this._dscal.Run(IMAX, ONE / SIGMA, ref W, offset_w, 1);
450: }
451: // *
452: // * Set WORK(1) to optimal workspace size.
453: // *
454: WORK[1 + o_work] = LWKOPT;
455: // *
456: return;
457: // *
458: // * End of DSYEV
459: // *
460:
461: #endregion
462:
463: }
464: }
465: }