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