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: /// DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
27: /// real symmetric matrix A. If eigenvectors are desired, it uses a
28: /// divide and conquer algorithm.
29: ///
30: /// The divide and conquer algorithm makes very mild assumptions about
31: /// floating point arithmetic. It will work on machines with a guard
32: /// digit in add/subtract, or on those binary machines without guard
33: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
34: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
35: /// without guard digits, but we know of none.
36: ///
37: /// Because of large use of BLAS of level 3, DSYEVD needs N**2 more
38: /// workspace than DSYEVX.
39: ///
40: ///</summary>
41: public class DSYEVD
42: {
43:
44:
45: #region Dependencies
46:
47: LSAME _lsame; DLAMCH _dlamch; DLANSY _dlansy; ILAENV _ilaenv; DLACPY _dlacpy; DLASCL _dlascl; DORMTR _dormtr;
48: DSCAL _dscal;DSTEDC _dstedc; DSTERF _dsterf; DSYTRD _dsytrd; XERBLA _xerbla;
49:
50: #endregion
51:
52:
53: #region Fields
54:
55: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LOWER = false; bool LQUERY = false; bool WANTZ = false;
56: int IINFO = 0;int INDE = 0; int INDTAU = 0; int INDWK2 = 0; int INDWRK = 0; int ISCALE = 0; int LIOPT = 0; int LIWMIN = 0;
57: int LLWORK = 0;int LLWRK2 = 0; int LOPT = 0; int LWMIN = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0;
58: double RMAX = 0;double RMIN = 0; double SAFMIN = 0; double SIGMA = 0; double SMLNUM = 0;
59:
60: #endregion
61:
62: public DSYEVD(LSAME lsame, DLAMCH dlamch, DLANSY dlansy, ILAENV ilaenv, DLACPY dlacpy, DLASCL dlascl, DORMTR dormtr, DSCAL dscal, DSTEDC dstedc, DSTERF dsterf
63: , DSYTRD dsytrd, XERBLA xerbla)
64: {
65:
66:
67: #region Set Dependencies
68:
69: this._lsame = lsame; this._dlamch = dlamch; this._dlansy = dlansy; this._ilaenv = ilaenv; this._dlacpy = dlacpy;
70: this._dlascl = dlascl;this._dormtr = dormtr; this._dscal = dscal; this._dstedc = dstedc; this._dsterf = dsterf;
71: this._dsytrd = dsytrd;this._xerbla = xerbla;
72:
73: #endregion
74:
75: }
76:
77: public DSYEVD()
78: {
79:
80:
81: #region Dependencies (Initialization)
82:
83: LSAME lsame = new LSAME();
84: DLAMC3 dlamc3 = new DLAMC3();
85: DLASSQ dlassq = new DLASSQ();
86: IEEECK ieeeck = new IEEECK();
87: IPARMQ iparmq = new IPARMQ();
88: XERBLA xerbla = new XERBLA();
89: DCOPY dcopy = new DCOPY();
90: DSCAL dscal = new DSCAL();
91: IDAMAX idamax = new IDAMAX();
92: DLAPY2 dlapy2 = new DLAPY2();
93: DLAMRG dlamrg = new DLAMRG();
94: DROT drot = new DROT();
95: DNRM2 dnrm2 = new DNRM2();
96: DLAED5 dlaed5 = new DLAED5();
97: DLAE2 dlae2 = new DLAE2();
98: DLAEV2 dlaev2 = new DLAEV2();
99: DSWAP dswap = new DSWAP();
100: DAXPY daxpy = new DAXPY();
101: DDOT ddot = new DDOT();
102: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
103: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
104: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
105: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
106: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
107: DLANSY dlansy = new DLANSY(dlassq, lsame);
108: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
109: DLACPY dlacpy = new DLACPY(lsame);
110: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
111: DGEMM dgemm = new DGEMM(lsame, xerbla);
112: DTRMM dtrmm = new DTRMM(lsame, xerbla);
113: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
114: DGEMV dgemv = new DGEMV(lsame, xerbla);
115: DTRMV dtrmv = new DTRMV(lsame, xerbla);
116: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
117: DGER dger = new DGER(xerbla);
118: DLARF dlarf = new DLARF(dgemv, dger, lsame);
119: DORM2L dorm2l = new DORM2L(lsame, dlarf, xerbla);
120: DORMQL dormql = new DORMQL(lsame, ilaenv, dlarfb, dlarft, dorm2l, xerbla);
121: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
122: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
123: DORMTR dormtr = new DORMTR(lsame, ilaenv, dormql, dormqr, xerbla);
124: DLANST dlanst = new DLANST(lsame, dlassq);
125: DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
126: DLAED6 dlaed6 = new DLAED6(dlamch);
127: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
128: DLASET dlaset = new DLASET(lsame);
129: DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
130: DLAED1 dlaed1 = new DLAED1(dcopy, dlaed2, dlaed3, dlamrg, xerbla);
131: DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
132: DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
133: DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
134: DLAED7 dlaed7 = new DLAED7(dgemm, dlaed8, dlaed9, dlaeda, dlamrg, xerbla);
135: DLARTG dlartg = new DLARTG(dlamch);
136: DLASR dlasr = new DLASR(lsame, xerbla);
137: DLASRT dlasrt = new DLASRT(lsame, xerbla);
138: DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
139: , dlasrt, dswap, xerbla);
140: DLAED0 dlaed0 = new DLAED0(dcopy, dgemm, dlacpy, dlaed1, dlaed7, dsteqr, xerbla, ilaenv);
141: DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
142: DSTEDC dstedc = new DSTEDC(lsame, ilaenv, dlamch, dlanst, dgemm, dlacpy, dlaed0, dlascl, dlaset, dlasrt
143: , dsteqr, dsterf, dswap, xerbla);
144: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
145: DSYMV dsymv = new DSYMV(lsame, xerbla);
146: DLATRD dlatrd = new DLATRD(daxpy, dgemv, dlarfg, dscal, dsymv, lsame, ddot);
147: DSYR2K dsyr2k = new DSYR2K(lsame, xerbla);
148: DSYR2 dsyr2 = new DSYR2(lsame, xerbla);
149: DSYTD2 dsytd2 = new DSYTD2(daxpy, dlarfg, dsymv, dsyr2, xerbla, lsame, ddot);
150: DSYTRD dsytrd = new DSYTRD(dlatrd, dsyr2k, dsytd2, xerbla, lsame, ilaenv);
151:
152: #endregion
153:
154:
155: #region Set Dependencies
156:
157: this._lsame = lsame; this._dlamch = dlamch; this._dlansy = dlansy; this._ilaenv = ilaenv; this._dlacpy = dlacpy;
158: this._dlascl = dlascl;this._dormtr = dormtr; this._dscal = dscal; this._dstedc = dstedc; this._dsterf = dsterf;
159: this._dsytrd = dsytrd;this._xerbla = xerbla;
160:
161: #endregion
162:
163: }
164: /// <summary>
165: /// Purpose
166: /// =======
167: ///
168: /// DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
169: /// real symmetric matrix A. If eigenvectors are desired, it uses a
170: /// divide and conquer algorithm.
171: ///
172: /// The divide and conquer algorithm makes very mild assumptions about
173: /// floating point arithmetic. It will work on machines with a guard
174: /// digit in add/subtract, or on those binary machines without guard
175: /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
176: /// Cray-2. It could conceivably fail on hexadecimal or decimal machines
177: /// without guard digits, but we know of none.
178: ///
179: /// Because of large use of BLAS of level 3, DSYEVD needs N**2 more
180: /// workspace than DSYEVX.
181: ///
182: ///</summary>
183: /// <param name="JOBZ">
184: /// (input) CHARACTER*1
185: /// = 'N': Compute eigenvalues only;
186: /// = 'V': Compute eigenvalues and eigenvectors.
187: ///</param>
188: /// <param name="UPLO">
189: /// (input) CHARACTER*1
190: /// = 'U': Upper triangle of A is stored;
191: /// = 'L': Lower triangle of A is stored.
192: ///</param>
193: /// <param name="N">
194: /// (input) INTEGER
195: /// The order of the matrix A. N .GE. 0.
196: ///</param>
197: /// <param name="A">
198: /// (input/output) DOUBLE PRECISION array, dimension (LDA, N)
199: /// On entry, the symmetric matrix A. If UPLO = 'U', the
200: /// leading N-by-N upper triangular part of A contains the
201: /// upper triangular part of the matrix A. If UPLO = 'L',
202: /// the leading N-by-N lower triangular part of A contains
203: /// the lower triangular part of the matrix A.
204: /// On exit, if JOBZ = 'V', then if INFO = 0, A contains the
205: /// orthonormal eigenvectors of the matrix A.
206: /// If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
207: /// or the upper triangle (if UPLO='U') of A, including the
208: /// diagonal, is destroyed.
209: ///</param>
210: /// <param name="LDA">
211: /// (input) INTEGER
212: /// The leading dimension of the array A. LDA .GE. max(1,N).
213: ///</param>
214: /// <param name="W">
215: /// (output) DOUBLE PRECISION array, dimension (N)
216: /// If INFO = 0, the eigenvalues in ascending order.
217: ///</param>
218: /// <param name="WORK">
219: /// (workspace/output) DOUBLE PRECISION array,
220: /// dimension (LWORK)
221: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
222: ///</param>
223: /// <param name="LWORK">
224: /// (input) INTEGER
225: /// The dimension of the array WORK.
226: /// If N .LE. 1, LWORK must be at least 1.
227: /// If JOBZ = 'N' and N .GT. 1, LWORK must be at least 2*N+1.
228: /// If JOBZ = 'V' and N .GT. 1, LWORK must be at least
229: /// 1 + 6*N + 2*N**2.
230: ///
231: /// If LWORK = -1, then a workspace query is assumed; the routine
232: /// only calculates the optimal sizes of the WORK and IWORK
233: /// arrays, returns these values as the first entries of the WORK
234: /// and IWORK arrays, and no error message related to LWORK or
235: /// LIWORK is issued by XERBLA.
236: ///</param>
237: /// <param name="IWORK">
238: /// (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
239: /// On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
240: ///</param>
241: /// <param name="LIWORK">
242: /// (input) INTEGER
243: /// The dimension of the array IWORK.
244: /// If N .LE. 1, LIWORK must be at least 1.
245: /// If JOBZ = 'N' and N .GT. 1, LIWORK must be at least 1.
246: /// If JOBZ = 'V' and N .GT. 1, LIWORK must be at least 3 + 5*N.
247: ///
248: /// If LIWORK = -1, then a workspace query is assumed; the
249: /// routine only calculates the optimal sizes of the WORK and
250: /// IWORK arrays, returns these values as the first entries of
251: /// the WORK and IWORK arrays, and no error message related to
252: /// LWORK or LIWORK is issued by XERBLA.
253: ///</param>
254: /// <param name="INFO">
255: /// (output) INTEGER
256: /// = 0: successful exit
257: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
258: /// .GT. 0: if INFO = i and JOBZ = 'N', then the algorithm failed
259: /// to converge; i off-diagonal elements of an intermediate
260: /// tridiagonal form did not converge to zero;
261: /// if INFO = i and JOBZ = 'V', then the algorithm failed
262: /// to compute an eigenvalue while working on the submatrix
263: /// lying in rows and columns INFO/(N+1) through
264: /// mod(INFO,N+1).
265: ///</param>
266: public void Run(string JOBZ, string UPLO, int N, ref double[] A, int offset_a, int LDA, ref double[] W, int offset_w
267: , ref double[] WORK, int offset_work, int LWORK, ref int[] IWORK, int offset_iwork, int LIWORK, ref int INFO)
268: {
269:
270: #region Array Index Correction
271:
272: int o_a = -1 - LDA + offset_a; int o_w = -1 + offset_w; int o_work = -1 + offset_work;
273: int o_iwork = -1 + offset_iwork;
274:
275: #endregion
276:
277:
278: #region Strings
279:
280: JOBZ = JOBZ.Substring(0, 1); UPLO = UPLO.Substring(0, 1);
281:
282: #endregion
283:
284:
285: #region Prolog
286:
287: // *
288: // * -- LAPACK driver routine (version 3.1) --
289: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
290: // * November 2006
291: // *
292: // * .. Scalar Arguments ..
293: // * ..
294: // * .. Array Arguments ..
295: // * ..
296: // *
297: // * Purpose
298: // * =======
299: // *
300: // * DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
301: // * real symmetric matrix A. If eigenvectors are desired, it uses a
302: // * divide and conquer algorithm.
303: // *
304: // * The divide and conquer algorithm makes very mild assumptions about
305: // * floating point arithmetic. It will work on machines with a guard
306: // * digit in add/subtract, or on those binary machines without guard
307: // * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
308: // * Cray-2. It could conceivably fail on hexadecimal or decimal machines
309: // * without guard digits, but we know of none.
310: // *
311: // * Because of large use of BLAS of level 3, DSYEVD needs N**2 more
312: // * workspace than DSYEVX.
313: // *
314: // * Arguments
315: // * =========
316: // *
317: // * JOBZ (input) CHARACTER*1
318: // * = 'N': Compute eigenvalues only;
319: // * = 'V': Compute eigenvalues and eigenvectors.
320: // *
321: // * UPLO (input) CHARACTER*1
322: // * = 'U': Upper triangle of A is stored;
323: // * = 'L': Lower triangle of A is stored.
324: // *
325: // * N (input) INTEGER
326: // * The order of the matrix A. N >= 0.
327: // *
328: // * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
329: // * On entry, the symmetric matrix A. If UPLO = 'U', the
330: // * leading N-by-N upper triangular part of A contains the
331: // * upper triangular part of the matrix A. If UPLO = 'L',
332: // * the leading N-by-N lower triangular part of A contains
333: // * the lower triangular part of the matrix A.
334: // * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
335: // * orthonormal eigenvectors of the matrix A.
336: // * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
337: // * or the upper triangle (if UPLO='U') of A, including the
338: // * diagonal, is destroyed.
339: // *
340: // * LDA (input) INTEGER
341: // * The leading dimension of the array A. LDA >= max(1,N).
342: // *
343: // * W (output) DOUBLE PRECISION array, dimension (N)
344: // * If INFO = 0, the eigenvalues in ascending order.
345: // *
346: // * WORK (workspace/output) DOUBLE PRECISION array,
347: // * dimension (LWORK)
348: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
349: // *
350: // * LWORK (input) INTEGER
351: // * The dimension of the array WORK.
352: // * If N <= 1, LWORK must be at least 1.
353: // * If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
354: // * If JOBZ = 'V' and N > 1, LWORK must be at least
355: // * 1 + 6*N + 2*N**2.
356: // *
357: // * If LWORK = -1, then a workspace query is assumed; the routine
358: // * only calculates the optimal sizes of the WORK and IWORK
359: // * arrays, returns these values as the first entries of the WORK
360: // * and IWORK arrays, and no error message related to LWORK or
361: // * LIWORK is issued by XERBLA.
362: // *
363: // * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
364: // * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
365: // *
366: // * LIWORK (input) INTEGER
367: // * The dimension of the array IWORK.
368: // * If N <= 1, LIWORK must be at least 1.
369: // * If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
370: // * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
371: // *
372: // * If LIWORK = -1, then a workspace query is assumed; the
373: // * routine only calculates the optimal sizes of the WORK and
374: // * IWORK arrays, returns these values as the first entries of
375: // * the WORK and IWORK arrays, and no error message related to
376: // * LWORK or LIWORK is issued by XERBLA.
377: // *
378: // * INFO (output) INTEGER
379: // * = 0: successful exit
380: // * < 0: if INFO = -i, the i-th argument had an illegal value
381: // * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
382: // * to converge; i off-diagonal elements of an intermediate
383: // * tridiagonal form did not converge to zero;
384: // * if INFO = i and JOBZ = 'V', then the algorithm failed
385: // * to compute an eigenvalue while working on the submatrix
386: // * lying in rows and columns INFO/(N+1) through
387: // * mod(INFO,N+1).
388: // *
389: // * Further Details
390: // * ===============
391: // *
392: // * Based on contributions by
393: // * Jeff Rutter, Computer Science Division, University of California
394: // * at Berkeley, USA
395: // * Modified by Francoise Tisseur, University of Tennessee.
396: // *
397: // * Modified description of INFO. Sven, 16 Feb 05.
398: // * =====================================================================
399: // *
400: // * .. Parameters ..
401: // * ..
402: // * .. Local Scalars ..
403: // *
404: // * ..
405: // * .. External Functions ..
406: // * ..
407: // * .. External Subroutines ..
408: // * ..
409: // * .. Intrinsic Functions ..
410: // INTRINSIC MAX, SQRT;
411: // * ..
412: // * .. Executable Statements ..
413: // *
414: // * Test the input parameters.
415: // *
416:
417: #endregion
418:
419:
420: #region Body
421:
422: WANTZ = this._lsame.Run(JOBZ, "V");
423: LOWER = this._lsame.Run(UPLO, "L");
424: LQUERY = (LWORK == - 1 || LIWORK == - 1);
425: // *
426: INFO = 0;
427: if (!(WANTZ || this._lsame.Run(JOBZ, "N")))
428: {
429: INFO = - 1;
430: }
431: else
432: {
433: if (!(LOWER || this._lsame.Run(UPLO, "U")))
434: {
435: INFO = - 2;
436: }
437: else
438: {
439: if (N < 0)
440: {
441: INFO = - 3;
442: }
443: else
444: {
445: if (LDA < Math.Max(1, N))
446: {
447: INFO = - 5;
448: }
449: }
450: }
451: }
452: // *
453: if (INFO == 0)
454: {
455: if (N <= 1)
456: {
457: LIWMIN = 1;
458: LWMIN = 1;
459: LOPT = LWMIN;
460: LIOPT = LIWMIN;
461: }
462: else
463: {
464: if (WANTZ)
465: {
466: LIWMIN = 3 + 5 * N;
467: LWMIN = 1 + 6 * N + 2 * (int)Math.Pow(N, 2);
468: }
469: else
470: {
471: LIWMIN = 1;
472: LWMIN = 2 * N + 1;
473: }
474: LOPT = Math.Max(LWMIN, 2 * N + this._ilaenv.Run(1, "DSYTRD", UPLO, N, - 1, - 1, - 1));
475: LIOPT = LIWMIN;
476: }
477: WORK[1 + o_work] = LOPT;
478: IWORK[1 + o_iwork] = LIOPT;
479: // *
480: if (LWORK < LWMIN && !LQUERY)
481: {
482: INFO = - 8;
483: }
484: else
485: {
486: if (LIWORK < LIWMIN && !LQUERY)
487: {
488: INFO = - 10;
489: }
490: }
491: }
492: // *
493: if (INFO != 0)
494: {
495: this._xerbla.Run("DSYEVD", - INFO);
496: return;
497: }
498: else
499: {
500: if (LQUERY)
501: {
502: return;
503: }
504: }
505: // *
506: // * Quick return if possible
507: // *
508: if (N == 0) return;
509: // *
510: if (N == 1)
511: {
512: W[1 + o_w] = A[1+1 * LDA + o_a];
513: if (WANTZ) A[1+1 * LDA + o_a] = ONE;
514: return;
515: }
516: // *
517: // * Get machine constants.
518: // *
519: SAFMIN = this._dlamch.Run("Safe minimum");
520: EPS = this._dlamch.Run("Precision");
521: SMLNUM = SAFMIN / EPS;
522: BIGNUM = ONE / SMLNUM;
523: RMIN = Math.Sqrt(SMLNUM);
524: RMAX = Math.Sqrt(BIGNUM);
525: // *
526: // * Scale matrix to allowable range, if necessary.
527: // *
528: ANRM = this._dlansy.Run("M", UPLO, N, A, offset_a, LDA, ref WORK, offset_work);
529: ISCALE = 0;
530: if (ANRM > ZERO && ANRM < RMIN)
531: {
532: ISCALE = 1;
533: SIGMA = RMIN / ANRM;
534: }
535: else
536: {
537: if (ANRM > RMAX)
538: {
539: ISCALE = 1;
540: SIGMA = RMAX / ANRM;
541: }
542: }
543: if (ISCALE == 1)
544: {
545: this._dlascl.Run(UPLO, 0, 0, ONE, SIGMA, N
546: , N, ref A, offset_a, LDA, ref INFO);
547: }
548: // *
549: // * Call DSYTRD to reduce symmetric matrix to tridiagonal form.
550: // *
551: INDE = 1;
552: INDTAU = INDE + N;
553: INDWRK = INDTAU + N;
554: LLWORK = LWORK - INDWRK + 1;
555: INDWK2 = INDWRK + N * N;
556: LLWRK2 = LWORK - INDWK2 + 1;
557: // *
558: this._dsytrd.Run(UPLO, N, ref A, offset_a, LDA, ref W, offset_w, ref WORK, INDE + o_work
559: , ref WORK, INDTAU + o_work, ref WORK, INDWRK + o_work, LLWORK, ref IINFO);
560: LOPT = 2 * N + (int)WORK[INDWRK + o_work];
561: // *
562: // * For eigenvalues only, call DSTERF. For eigenvectors, first call
563: // * DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
564: // * tridiagonal matrix, then call DORMTR to multiply it by the
565: // * Householder transformations stored in A.
566: // *
567: if (!WANTZ)
568: {
569: this._dsterf.Run(N, ref W, offset_w, ref WORK, INDE + o_work, ref INFO);
570: }
571: else
572: {
573: this._dstedc.Run("I", N, ref W, offset_w, ref WORK, INDE + o_work, ref WORK, INDWRK + o_work, N
574: , ref WORK, INDWK2 + o_work, LLWRK2, ref IWORK, offset_iwork, LIWORK, ref INFO);
575: this._dormtr.Run("L", UPLO, "N", N, N, ref A, offset_a
576: , LDA, WORK, INDTAU + o_work, ref WORK, INDWRK + o_work, N, ref WORK, INDWK2 + o_work, LLWRK2
577: , ref IINFO);
578: this._dlacpy.Run("A", N, N, WORK, INDWRK + o_work, N, ref A, offset_a
579: , LDA);
580: LOPT = (int)Math.Max(LOPT, 1 + 6 * N + 2 * Math.Pow(N, 2));
581: }
582: // *
583: // * If matrix was scaled, then rescale eigenvalues appropriately.
584: // *
585: if (ISCALE == 1) this._dscal.Run(N, ONE / SIGMA, ref W, offset_w, 1);
586: // *
587: WORK[1 + o_work] = LOPT;
588: IWORK[1 + o_iwork] = LIOPT;
589: // *
590: return;
591: // *
592: // * End of DSYEVD
593: // *
594:
595: #endregion
596:
597: }
598: }
599: }