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: /// DGEEV computes for an N-by-N real nonsymmetric matrix A, the
27: /// eigenvalues and, optionally, the left and/or right eigenvectors.
28: ///
29: /// The right eigenvector v(j) of A satisfies
30: /// A * v(j) = lambda(j) * v(j)
31: /// where lambda(j) is its eigenvalue.
32: /// The left eigenvector u(j) of A satisfies
33: /// u(j)**H * A = lambda(j) * u(j)**H
34: /// where u(j)**H denotes the conjugate transpose of u(j).
35: ///
36: /// The computed eigenvectors are normalized to have Euclidean norm
37: /// equal to 1 and largest component real.
38: ///
39: ///</summary>
40: public class DGEEV
41: {
42:
43:
44: #region Dependencies
45:
46: DGEBAK _dgebak; DGEBAL _dgebal; DGEHRD _dgehrd; DHSEQR _dhseqr; DLABAD _dlabad; DLACPY _dlacpy; DLARTG _dlartg;
47: DLASCL _dlascl;DORGHR _dorghr; DROT _drot; DSCAL _dscal; DTREVC _dtrevc; XERBLA _xerbla; LSAME _lsame; IDAMAX _idamax;
48: ILAENV _ilaenv;DLAMCH _dlamch; DLANGE _dlange; DLAPY2 _dlapy2; DNRM2 _dnrm2;
49:
50: #endregion
51:
52:
53: #region Fields
54:
55: const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool SCALEA = false; bool WANTVL = false;
56: bool WANTVR = false;string SIDE = new string(' ', 1); int HSWORK = 0; int I = 0; int IBAL = 0; int IERR = 0; int IHI = 0;
57: int ILO = 0;int ITAU = 0; int IWRK = 0; int K = 0; int MAXWRK = 0; int MINWRK = 0; int NOUT = 0; double ANRM = 0;
58: double BIGNUM = 0;double CS = 0; double CSCALE = 0; double EPS = 0; double R = 0; double SCL = 0; double SMLNUM = 0;
59: double SN = 0;bool[] SELECT = new bool[1]; int offset_select = 0;
60: double[] DUM = new double[1]; int offset_dum = 0;
61:
62: #endregion
63:
64: public DGEEV(DGEBAK dgebak, DGEBAL dgebal, DGEHRD dgehrd, DHSEQR dhseqr, DLABAD dlabad, DLACPY dlacpy, DLARTG dlartg, DLASCL dlascl, DORGHR dorghr, DROT drot
65: , DSCAL dscal, DTREVC dtrevc, XERBLA xerbla, LSAME lsame, IDAMAX idamax, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange, DLAPY2 dlapy2, DNRM2 dnrm2)
66: {
67:
68:
69: #region Set Dependencies
70:
71: this._dgebak = dgebak; this._dgebal = dgebal; this._dgehrd = dgehrd; this._dhseqr = dhseqr; this._dlabad = dlabad;
72: this._dlacpy = dlacpy;this._dlartg = dlartg; this._dlascl = dlascl; this._dorghr = dorghr; this._drot = drot;
73: this._dscal = dscal;this._dtrevc = dtrevc; this._xerbla = xerbla; this._lsame = lsame; this._idamax = idamax;
74: this._ilaenv = ilaenv;this._dlamch = dlamch; this._dlange = dlange; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2;
75:
76: #endregion
77:
78: }
79:
80: public DGEEV()
81: {
82:
83:
84: #region Dependencies (Initialization)
85:
86: LSAME lsame = new LSAME();
87: DSCAL dscal = new DSCAL();
88: DSWAP dswap = new DSWAP();
89: XERBLA xerbla = new XERBLA();
90: IDAMAX idamax = new IDAMAX();
91: DLAMC3 dlamc3 = new DLAMC3();
92: DAXPY daxpy = new DAXPY();
93: DLAPY2 dlapy2 = new DLAPY2();
94: DNRM2 dnrm2 = new DNRM2();
95: DCOPY dcopy = new DCOPY();
96: IEEECK ieeeck = new IEEECK();
97: IPARMQ iparmq = new IPARMQ();
98: DLABAD dlabad = new DLABAD();
99: DROT drot = new DROT();
100: DLASSQ dlassq = new DLASSQ();
101: DLAQR1 dlaqr1 = new DLAQR1();
102: DDOT ddot = new DDOT();
103: DLADIV dladiv = new DLADIV();
104: DGEBAK dgebak = new DGEBAK(lsame, dscal, dswap, xerbla);
105: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
106: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
107: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
108: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
109: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
110: DGEBAL dgebal = new DGEBAL(lsame, idamax, dlamch, dscal, dswap, xerbla);
111: DGEMV dgemv = new DGEMV(lsame, xerbla);
112: DGER dger = new DGER(xerbla);
113: DLARF dlarf = new DLARF(dgemv, dger, lsame);
114: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
115: DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
116: DGEMM dgemm = new DGEMM(lsame, xerbla);
117: DLACPY dlacpy = new DLACPY(lsame);
118: DTRMM dtrmm = new DTRMM(lsame, xerbla);
119: DTRMV dtrmv = new DTRMV(lsame, xerbla);
120: DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
121: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
122: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
123: DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
124: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
125: DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
126: DLASET dlaset = new DLASET(lsame);
127: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
128: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
129: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
130: DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
131: DLANGE dlange = new DLANGE(dlassq, lsame);
132: DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
133: DLARTG dlartg = new DLARTG(dlamch);
134: DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
135: DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
136: DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
137: DLAQR2 dlaqr2 = new DLAQR2(dlamch, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlarf, dlarfg
138: , dlaset, dorghr, dtrexc);
139: DLAQR5 dlaqr5 = new DLAQR5(dlamch, dgemm, dlabad, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm);
140: DLAQR4 dlaqr4 = new DLAQR4(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5);
141: DLAQR3 dlaqr3 = new DLAQR3(dlamch, ilaenv, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlaqr4
142: , dlarf, dlarfg, dlaset, dorghr, dtrexc);
143: DLAQR0 dlaqr0 = new DLAQR0(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr3, dlaqr4, dlaqr5);
144: DHSEQR dhseqr = new DHSEQR(ilaenv, lsame, dlacpy, dlahqr, dlaqr0, dlaset, xerbla);
145: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
146: DLALN2 dlaln2 = new DLALN2(dlamch, dladiv);
147: DTREVC dtrevc = new DTREVC(lsame, idamax, ddot, dlamch, daxpy, dcopy, dgemv, dlaln2, dscal, xerbla
148: , dlabad);
149:
150: #endregion
151:
152:
153: #region Set Dependencies
154:
155: this._dgebak = dgebak; this._dgebal = dgebal; this._dgehrd = dgehrd; this._dhseqr = dhseqr; this._dlabad = dlabad;
156: this._dlacpy = dlacpy;this._dlartg = dlartg; this._dlascl = dlascl; this._dorghr = dorghr; this._drot = drot;
157: this._dscal = dscal;this._dtrevc = dtrevc; this._xerbla = xerbla; this._lsame = lsame; this._idamax = idamax;
158: this._ilaenv = ilaenv;this._dlamch = dlamch; this._dlange = dlange; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2;
159:
160: #endregion
161:
162: }
163: /// <summary>
164: /// Purpose
165: /// =======
166: ///
167: /// DGEEV computes for an N-by-N real nonsymmetric matrix A, the
168: /// eigenvalues and, optionally, the left and/or right eigenvectors.
169: ///
170: /// The right eigenvector v(j) of A satisfies
171: /// A * v(j) = lambda(j) * v(j)
172: /// where lambda(j) is its eigenvalue.
173: /// The left eigenvector u(j) of A satisfies
174: /// u(j)**H * A = lambda(j) * u(j)**H
175: /// where u(j)**H denotes the conjugate transpose of u(j).
176: ///
177: /// The computed eigenvectors are normalized to have Euclidean norm
178: /// equal to 1 and largest component real.
179: ///
180: ///</summary>
181: /// <param name="JOBVL">
182: /// (input) CHARACTER*1
183: /// = 'N': left eigenvectors of A are not computed;
184: /// = 'V': left eigenvectors of A are computed.
185: ///</param>
186: /// <param name="JOBVR">
187: /// (input) CHARACTER*1
188: /// = 'N': right eigenvectors of A are not computed;
189: /// = 'V': right eigenvectors of A are computed.
190: ///</param>
191: /// <param name="N">
192: /// (input) INTEGER
193: /// The order of the matrix A. N .GE. 0.
194: ///</param>
195: /// <param name="A">
196: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
197: /// On entry, the N-by-N matrix A.
198: /// On exit, A has been overwritten.
199: ///</param>
200: /// <param name="LDA">
201: /// (input) INTEGER
202: /// The leading dimension of the array A. LDA .GE. max(1,N).
203: ///</param>
204: /// <param name="WR">
205: /// (output) DOUBLE PRECISION array, dimension (N)
206: ///</param>
207: /// <param name="WI">
208: /// (output) DOUBLE PRECISION array, dimension (N)
209: /// WR and WI contain the real and imaginary parts,
210: /// respectively, of the computed eigenvalues. Complex
211: /// conjugate pairs of eigenvalues appear consecutively
212: /// with the eigenvalue having the positive imaginary part
213: /// first.
214: ///</param>
215: /// <param name="VL">
216: /// (output) DOUBLE PRECISION array, dimension (LDVL,N)
217: /// If JOBVL = 'V', the left eigenvectors u(j) are stored one
218: /// after another in the columns of VL, in the same order
219: /// as their eigenvalues.
220: /// If JOBVL = 'N', VL is not referenced.
221: /// If the j-th eigenvalue is real, then u(j) = VL(:,j),
222: /// the j-th column of VL.
223: /// If the j-th and (j+1)-st eigenvalues form a complex
224: /// conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
225: /// u(j+1) = VL(:,j) - i*VL(:,j+1).
226: ///</param>
227: /// <param name="LDVL">
228: /// (input) INTEGER
229: /// The leading dimension of the array VL. LDVL .GE. 1; if
230: /// JOBVL = 'V', LDVL .GE. N.
231: ///</param>
232: /// <param name="VR">
233: /// (output) DOUBLE PRECISION array, dimension (LDVR,N)
234: /// If JOBVR = 'V', the right eigenvectors v(j) are stored one
235: /// after another in the columns of VR, in the same order
236: /// as their eigenvalues.
237: /// If JOBVR = 'N', VR is not referenced.
238: /// If the j-th eigenvalue is real, then v(j) = VR(:,j),
239: /// the j-th column of VR.
240: /// If the j-th and (j+1)-st eigenvalues form a complex
241: /// conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
242: /// v(j+1) = VR(:,j) - i*VR(:,j+1).
243: ///</param>
244: /// <param name="LDVR">
245: /// (input) INTEGER
246: /// The leading dimension of the array VR. LDVR .GE. 1; if
247: /// JOBVR = 'V', LDVR .GE. N.
248: ///</param>
249: /// <param name="WORK">
250: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
251: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
252: ///</param>
253: /// <param name="LWORK">
254: /// (input) INTEGER
255: /// The dimension of the array WORK. LWORK .GE. max(1,3*N), and
256: /// if JOBVL = 'V' or JOBVR = 'V', LWORK .GE. 4*N. For good
257: /// performance, LWORK must generally be larger.
258: ///
259: /// If LWORK = -1, then a workspace query is assumed; the routine
260: /// only calculates the optimal size of the WORK array, returns
261: /// this value as the first entry of the WORK array, and no error
262: /// message related to LWORK is issued by XERBLA.
263: ///</param>
264: /// <param name="INFO">
265: /// (output) INTEGER
266: /// = 0: successful exit
267: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
268: /// .GT. 0: if INFO = i, the QR algorithm failed to compute all the
269: /// eigenvalues, and no eigenvectors have been computed;
270: /// elements i+1:N of WR and WI contain eigenvalues which
271: /// have converged.
272: ///</param>
273: public void Run(string JOBVL, string JOBVR, int N, ref double[] A, int offset_a, int LDA, ref double[] WR, int offset_wr
274: , ref double[] WI, int offset_wi, ref double[] VL, int offset_vl, int LDVL, ref double[] VR, int offset_vr, int LDVR, ref double[] WORK, int offset_work
275: , int LWORK, ref int INFO)
276: {
277:
278: #region Array Index Correction
279:
280: int o_a = -1 - LDA + offset_a; int o_wr = -1 + offset_wr; int o_wi = -1 + offset_wi;
281: int o_vl = -1 - LDVL + offset_vl; int o_vr = -1 - LDVR + offset_vr; int o_work = -1 + offset_work;
282:
283: #endregion
284:
285:
286: #region Strings
287:
288: JOBVL = JOBVL.Substring(0, 1); JOBVR = JOBVR.Substring(0, 1);
289:
290: #endregion
291:
292:
293: #region Prolog
294:
295: // *
296: // * -- LAPACK driver routine (version 3.1) --
297: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
298: // * November 2006
299: // *
300: // * .. Scalar Arguments ..
301: // * ..
302: // * .. Array Arguments ..
303: // * ..
304: // *
305: // * Purpose
306: // * =======
307: // *
308: // * DGEEV computes for an N-by-N real nonsymmetric matrix A, the
309: // * eigenvalues and, optionally, the left and/or right eigenvectors.
310: // *
311: // * The right eigenvector v(j) of A satisfies
312: // * A * v(j) = lambda(j) * v(j)
313: // * where lambda(j) is its eigenvalue.
314: // * The left eigenvector u(j) of A satisfies
315: // * u(j)**H * A = lambda(j) * u(j)**H
316: // * where u(j)**H denotes the conjugate transpose of u(j).
317: // *
318: // * The computed eigenvectors are normalized to have Euclidean norm
319: // * equal to 1 and largest component real.
320: // *
321: // * Arguments
322: // * =========
323: // *
324: // * JOBVL (input) CHARACTER*1
325: // * = 'N': left eigenvectors of A are not computed;
326: // * = 'V': left eigenvectors of A are computed.
327: // *
328: // * JOBVR (input) CHARACTER*1
329: // * = 'N': right eigenvectors of A are not computed;
330: // * = 'V': right eigenvectors of A are computed.
331: // *
332: // * N (input) INTEGER
333: // * The order of the matrix A. N >= 0.
334: // *
335: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
336: // * On entry, the N-by-N matrix A.
337: // * On exit, A has been overwritten.
338: // *
339: // * LDA (input) INTEGER
340: // * The leading dimension of the array A. LDA >= max(1,N).
341: // *
342: // * WR (output) DOUBLE PRECISION array, dimension (N)
343: // * WI (output) DOUBLE PRECISION array, dimension (N)
344: // * WR and WI contain the real and imaginary parts,
345: // * respectively, of the computed eigenvalues. Complex
346: // * conjugate pairs of eigenvalues appear consecutively
347: // * with the eigenvalue having the positive imaginary part
348: // * first.
349: // *
350: // * VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
351: // * If JOBVL = 'V', the left eigenvectors u(j) are stored one
352: // * after another in the columns of VL, in the same order
353: // * as their eigenvalues.
354: // * If JOBVL = 'N', VL is not referenced.
355: // * If the j-th eigenvalue is real, then u(j) = VL(:,j),
356: // * the j-th column of VL.
357: // * If the j-th and (j+1)-st eigenvalues form a complex
358: // * conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
359: // * u(j+1) = VL(:,j) - i*VL(:,j+1).
360: // *
361: // * LDVL (input) INTEGER
362: // * The leading dimension of the array VL. LDVL >= 1; if
363: // * JOBVL = 'V', LDVL >= N.
364: // *
365: // * VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
366: // * If JOBVR = 'V', the right eigenvectors v(j) are stored one
367: // * after another in the columns of VR, in the same order
368: // * as their eigenvalues.
369: // * If JOBVR = 'N', VR is not referenced.
370: // * If the j-th eigenvalue is real, then v(j) = VR(:,j),
371: // * the j-th column of VR.
372: // * If the j-th and (j+1)-st eigenvalues form a complex
373: // * conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
374: // * v(j+1) = VR(:,j) - i*VR(:,j+1).
375: // *
376: // * LDVR (input) INTEGER
377: // * The leading dimension of the array VR. LDVR >= 1; if
378: // * JOBVR = 'V', LDVR >= N.
379: // *
380: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
381: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
382: // *
383: // * LWORK (input) INTEGER
384: // * The dimension of the array WORK. LWORK >= max(1,3*N), and
385: // * if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
386: // * performance, LWORK must generally be larger.
387: // *
388: // * If LWORK = -1, then a workspace query is assumed; the routine
389: // * only calculates the optimal size of the WORK array, returns
390: // * this value as the first entry of the WORK array, and no error
391: // * message related to LWORK is issued by XERBLA.
392: // *
393: // * INFO (output) INTEGER
394: // * = 0: successful exit
395: // * < 0: if INFO = -i, the i-th argument had an illegal value.
396: // * > 0: if INFO = i, the QR algorithm failed to compute all the
397: // * eigenvalues, and no eigenvectors have been computed;
398: // * elements i+1:N of WR and WI contain eigenvalues which
399: // * have converged.
400: // *
401: // * =====================================================================
402: // *
403: // * .. Parameters ..
404: // * ..
405: // * .. Local Scalars ..
406: // * ..
407: // * .. Local Arrays ..
408: // * ..
409: // * .. External Subroutines ..
410: // * ..
411: // * .. External Functions ..
412: // * ..
413: // * .. Intrinsic Functions ..
414: // INTRINSIC MAX, SQRT;
415: // * ..
416: // * .. Executable Statements ..
417: // *
418: // * Test the input arguments
419: // *
420:
421: #endregion
422:
423:
424: #region Body
425:
426: INFO = 0;
427: LQUERY = (LWORK == - 1);
428: WANTVL = this._lsame.Run(JOBVL, "V");
429: WANTVR = this._lsame.Run(JOBVR, "V");
430: if ((!WANTVL) && (!this._lsame.Run(JOBVL, "N")))
431: {
432: INFO = - 1;
433: }
434: else
435: {
436: if ((!WANTVR) && (!this._lsame.Run(JOBVR, "N")))
437: {
438: INFO = - 2;
439: }
440: else
441: {
442: if (N < 0)
443: {
444: INFO = - 3;
445: }
446: else
447: {
448: if (LDA < Math.Max(1, N))
449: {
450: INFO = - 5;
451: }
452: else
453: {
454: if (LDVL < 1 || (WANTVL && LDVL < N))
455: {
456: INFO = - 9;
457: }
458: else
459: {
460: if (LDVR < 1 || (WANTVR && LDVR < N))
461: {
462: INFO = - 11;
463: }
464: }
465: }
466: }
467: }
468: }
469: // *
470: // * Compute workspace
471: // * (Note: Comments in the code beginning "Workspace:" describe the
472: // * minimal amount of workspace needed at that point in the code,
473: // * as well as the preferred amount for good performance.
474: // * NB refers to the optimal block size for the immediately
475: // * following subroutine, as returned by ILAENV.
476: // * HSWORK refers to the workspace preferred by DHSEQR, as
477: // * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
478: // * the worst case.)
479: // *
480: if (INFO == 0)
481: {
482: if (N == 0)
483: {
484: MINWRK = 1;
485: MAXWRK = 1;
486: }
487: else
488: {
489: MAXWRK = 2 * N + N * this._ilaenv.Run(1, "DGEHRD", " ", N, 1, N, 0);
490: if (WANTVL)
491: {
492: MINWRK = 4 * N;
493: MAXWRK = Math.Max(MAXWRK, 2 * N + (N - 1) * this._ilaenv.Run(1, "DORGHR", " ", N, 1, N, - 1));
494: this._dhseqr.Run("S", "V", N, 1, N, ref A, offset_a
495: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VL, offset_vl, LDVL, ref WORK, offset_work
496: , - 1, ref INFO);
497: HSWORK = (int)WORK[1 + o_work];
498: MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
499: MAXWRK = Math.Max(MAXWRK, 4 * N);
500: }
501: else
502: {
503: if (WANTVR)
504: {
505: MINWRK = 4 * N;
506: MAXWRK = Math.Max(MAXWRK, 2 * N + (N - 1) * this._ilaenv.Run(1, "DORGHR", " ", N, 1, N, - 1));
507: this._dhseqr.Run("S", "V", N, 1, N, ref A, offset_a
508: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, offset_work
509: , - 1, ref INFO);
510: HSWORK = (int)WORK[1 + o_work];
511: MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
512: MAXWRK = Math.Max(MAXWRK, 4 * N);
513: }
514: else
515: {
516: MINWRK = 3 * N;
517: this._dhseqr.Run("E", "N", N, 1, N, ref A, offset_a
518: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, offset_work
519: , - 1, ref INFO);
520: HSWORK = (int)WORK[1 + o_work];
521: MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
522: }
523: }
524: MAXWRK = Math.Max(MAXWRK, MINWRK);
525: }
526: WORK[1 + o_work] = MAXWRK;
527: // *
528: if (LWORK < MINWRK && !LQUERY)
529: {
530: INFO = - 13;
531: }
532: }
533: // *
534: if (INFO != 0)
535: {
536: this._xerbla.Run("DGEEV ", - INFO);
537: return;
538: }
539: else
540: {
541: if (LQUERY)
542: {
543: return;
544: }
545: }
546: // *
547: // * Quick return if possible
548: // *
549: if (N == 0) return;
550: // *
551: // * Get machine constants
552: // *
553: EPS = this._dlamch.Run("P");
554: SMLNUM = this._dlamch.Run("S");
555: BIGNUM = ONE / SMLNUM;
556: this._dlabad.Run(ref SMLNUM, ref BIGNUM);
557: SMLNUM = Math.Sqrt(SMLNUM) / EPS;
558: BIGNUM = ONE / SMLNUM;
559: // *
560: // * Scale A if max element outside range [SMLNUM,BIGNUM]
561: // *
562: ANRM = this._dlange.Run("M", N, N, A, offset_a, LDA, ref DUM, offset_dum);
563: SCALEA = false;
564: if (ANRM > ZERO && ANRM < SMLNUM)
565: {
566: SCALEA = true;
567: CSCALE = SMLNUM;
568: }
569: else
570: {
571: if (ANRM > BIGNUM)
572: {
573: SCALEA = true;
574: CSCALE = BIGNUM;
575: }
576: }
577: if (SCALEA)
578: {
579: this._dlascl.Run("G", 0, 0, ANRM, CSCALE, N
580: , N, ref A, offset_a, LDA, ref IERR);
581: }
582: // *
583: // * Balance the matrix
584: // * (Workspace: need N)
585: // *
586: IBAL = 1;
587: this._dgebal.Run("B", N, ref A, offset_a, LDA, ref ILO, ref IHI
588: , ref WORK, IBAL + o_work, ref IERR);
589: // *
590: // * Reduce to upper Hessenberg form
591: // * (Workspace: need 3*N, prefer 2*N+N*NB)
592: // *
593: ITAU = IBAL + N;
594: IWRK = ITAU + N;
595: this._dgehrd.Run(N, ILO, IHI, ref A, offset_a, LDA, ref WORK, ITAU + o_work
596: , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
597: // *
598: if (WANTVL)
599: {
600: // *
601: // * Want left eigenvectors
602: // * Copy Householder vectors to VL
603: // *
604: FortranLib.Copy(ref SIDE , "L");
605: this._dlacpy.Run("L", N, N, A, offset_a, LDA, ref VL, offset_vl
606: , LDVL);
607: // *
608: // * Generate orthogonal matrix in VL
609: // * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
610: // *
611: this._dorghr.Run(N, ILO, IHI, ref VL, offset_vl, LDVL, WORK, ITAU + o_work
612: , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
613: // *
614: // * Perform QR iteration, accumulating Schur vectors in VL
615: // * (Workspace: need N+1, prefer N+HSWORK (see comments) )
616: // *
617: IWRK = ITAU;
618: this._dhseqr.Run("S", "V", N, ILO, IHI, ref A, offset_a
619: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VL, offset_vl, LDVL, ref WORK, IWRK + o_work
620: , LWORK - IWRK + 1, ref INFO);
621: // *
622: if (WANTVR)
623: {
624: // *
625: // * Want left and right eigenvectors
626: // * Copy Schur vectors to VR
627: // *
628: FortranLib.Copy(ref SIDE , "B");
629: this._dlacpy.Run("F", N, N, VL, offset_vl, LDVL, ref VR, offset_vr
630: , LDVR);
631: }
632: // *
633: }
634: else
635: {
636: if (WANTVR)
637: {
638: // *
639: // * Want right eigenvectors
640: // * Copy Householder vectors to VR
641: // *
642: FortranLib.Copy(ref SIDE , "R");
643: this._dlacpy.Run("L", N, N, A, offset_a, LDA, ref VR, offset_vr
644: , LDVR);
645: // *
646: // * Generate orthogonal matrix in VR
647: // * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
648: // *
649: this._dorghr.Run(N, ILO, IHI, ref VR, offset_vr, LDVR, WORK, ITAU + o_work
650: , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
651: // *
652: // * Perform QR iteration, accumulating Schur vectors in VR
653: // * (Workspace: need N+1, prefer N+HSWORK (see comments) )
654: // *
655: IWRK = ITAU;
656: this._dhseqr.Run("S", "V", N, ILO, IHI, ref A, offset_a
657: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, IWRK + o_work
658: , LWORK - IWRK + 1, ref INFO);
659: // *
660: }
661: else
662: {
663: // *
664: // * Compute eigenvalues only
665: // * (Workspace: need N+1, prefer N+HSWORK (see comments) )
666: // *
667: IWRK = ITAU;
668: this._dhseqr.Run("E", "N", N, ILO, IHI, ref A, offset_a
669: , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, IWRK + o_work
670: , LWORK - IWRK + 1, ref INFO);
671: }
672: }
673: // *
674: // * If INFO > 0 from DHSEQR, then quit
675: // *
676: if (INFO > 0) goto LABEL50;
677: // *
678: if (WANTVL || WANTVR)
679: {
680: // *
681: // * Compute left and/or right eigenvectors
682: // * (Workspace: need 4*N)
683: // *
684: this._dtrevc.Run(SIDE, "B", ref SELECT, offset_select, N, A, offset_a, LDA
685: , ref VL, offset_vl, LDVL, ref VR, offset_vr, LDVR, N, ref NOUT
686: , ref WORK, IWRK + o_work, ref IERR);
687: }
688: // *
689: if (WANTVL)
690: {
691: // *
692: // * Undo balancing of left eigenvectors
693: // * (Workspace: need N)
694: // *
695: this._dgebak.Run("B", "L", N, ILO, IHI, WORK, IBAL + o_work
696: , N, ref VL, offset_vl, LDVL, ref IERR);
697: // *
698: // * Normalize left eigenvectors and make largest component real
699: // *
700: for (I = 1; I <= N; I++)
701: {
702: if (WI[I + o_wi] == ZERO)
703: {
704: SCL = ONE / this._dnrm2.Run(N, VL, 1+I * LDVL + o_vl, 1);
705: this._dscal.Run(N, SCL, ref VL, 1+I * LDVL + o_vl, 1);
706: }
707: else
708: {
709: if (WI[I + o_wi] > ZERO)
710: {
711: SCL = ONE / this._dlapy2.Run(this._dnrm2.Run(N, VL, 1+I * LDVL + o_vl, 1), this._dnrm2.Run(N, VL, 1+(I + 1) * LDVL + o_vl, 1));
712: this._dscal.Run(N, SCL, ref VL, 1+I * LDVL + o_vl, 1);
713: this._dscal.Run(N, SCL, ref VL, 1+(I + 1) * LDVL + o_vl, 1);
714: for (K = 1; K <= N; K++)
715: {
716: WORK[IWRK + K - 1 + o_work] = Math.Pow(VL[K+I * LDVL + o_vl],2) + Math.Pow(VL[K+(I + 1) * LDVL + o_vl],2);
717: }
718: K = this._idamax.Run(N, WORK, IWRK + o_work, 1);
719: this._dlartg.Run(VL[K+I * LDVL + o_vl], VL[K+(I + 1) * LDVL + o_vl], ref CS, ref SN, ref R);
720: this._drot.Run(N, ref VL, 1+I * LDVL + o_vl, 1, ref VL, 1+(I + 1) * LDVL + o_vl, 1, CS
721: , SN);
722: VL[K+(I + 1) * LDVL + o_vl] = ZERO;
723: }
724: }
725: }
726: }
727: // *
728: if (WANTVR)
729: {
730: // *
731: // * Undo balancing of right eigenvectors
732: // * (Workspace: need N)
733: // *
734: this._dgebak.Run("B", "R", N, ILO, IHI, WORK, IBAL + o_work
735: , N, ref VR, offset_vr, LDVR, ref IERR);
736: // *
737: // * Normalize right eigenvectors and make largest component real
738: // *
739: for (I = 1; I <= N; I++)
740: {
741: if (WI[I + o_wi] == ZERO)
742: {
743: SCL = ONE / this._dnrm2.Run(N, VR, 1+I * LDVR + o_vr, 1);
744: this._dscal.Run(N, SCL, ref VR, 1+I * LDVR + o_vr, 1);
745: }
746: else
747: {
748: if (WI[I + o_wi] > ZERO)
749: {
750: SCL = ONE / this._dlapy2.Run(this._dnrm2.Run(N, VR, 1+I * LDVR + o_vr, 1), this._dnrm2.Run(N, VR, 1+(I + 1) * LDVR + o_vr, 1));
751: this._dscal.Run(N, SCL, ref VR, 1+I * LDVR + o_vr, 1);
752: this._dscal.Run(N, SCL, ref VR, 1+(I + 1) * LDVR + o_vr, 1);
753: for (K = 1; K <= N; K++)
754: {
755: WORK[IWRK + K - 1 + o_work] = Math.Pow(VR[K+I * LDVR + o_vr],2) + Math.Pow(VR[K+(I + 1) * LDVR + o_vr],2);
756: }
757: K = this._idamax.Run(N, WORK, IWRK + o_work, 1);
758: this._dlartg.Run(VR[K+I * LDVR + o_vr], VR[K+(I + 1) * LDVR + o_vr], ref CS, ref SN, ref R);
759: this._drot.Run(N, ref VR, 1+I * LDVR + o_vr, 1, ref VR, 1+(I + 1) * LDVR + o_vr, 1, CS
760: , SN);
761: VR[K+(I + 1) * LDVR + o_vr] = ZERO;
762: }
763: }
764: }
765: }
766: // *
767: // * Undo scaling if necessary
768: // *
769: LABEL50:;
770: if (SCALEA)
771: {
772: this._dlascl.Run("G", 0, 0, CSCALE, ANRM, N - INFO
773: , 1, ref WR, INFO + 1 + o_wr, Math.Max(N - INFO, 1), ref IERR);
774: this._dlascl.Run("G", 0, 0, CSCALE, ANRM, N - INFO
775: , 1, ref WI, INFO + 1 + o_wi, Math.Max(N - INFO, 1), ref IERR);
776: if (INFO > 0)
777: {
778: this._dlascl.Run("G", 0, 0, CSCALE, ANRM, ILO - 1
779: , 1, ref WR, offset_wr, N, ref IERR);
780: this._dlascl.Run("G", 0, 0, CSCALE, ANRM, ILO - 1
781: , 1, ref WI, offset_wi, N, ref IERR);
782: }
783: }
784: // *
785: WORK[1 + o_work] = MAXWRK;
786: return;
787: // *
788: // * End of DGEEV
789: // *
790:
791: #endregion
792:
793: }
794: }
795: }