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 routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLAED0 computes all eigenvalues and corresponding eigenvectors of a
27: /// symmetric tridiagonal matrix using the divide and conquer method.
28: ///
29: ///</summary>
30: public class DLAED0
31: {
32:
33:
34: #region Dependencies
35:
36: DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLAED1 _dlaed1; DLAED7 _dlaed7; DSTEQR _dsteqr; XERBLA _xerbla;
37: ILAENV _ilaenv;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int CURLVL = 0; int CURPRB = 0;
45: int CURR = 0;int I = 0; int IGIVCL = 0; int IGIVNM = 0; int IGIVPT = 0; int INDXQ = 0; int IPERM = 0; int IPRMPT = 0;
46: int IQ = 0;int IQPTR = 0; int IWREM = 0; int J = 0; int K = 0; int LGN = 0; int MATSIZ = 0; int MSD2 = 0; int SMLSIZ = 0;
47: int SMM1 = 0;int SPM1 = 0; int SPM2 = 0; int SUBMAT = 0; int SUBPBS = 0; int TLVLS = 0; double TEMP = 0;
48:
49: #endregion
50:
51: public DLAED0(DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLAED1 dlaed1, DLAED7 dlaed7, DSTEQR dsteqr, XERBLA xerbla, ILAENV ilaenv)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; this._dlaed1 = dlaed1; this._dlaed7 = dlaed7;
58: this._dsteqr = dsteqr;this._xerbla = xerbla; this._ilaenv = ilaenv;
59:
60: #endregion
61:
62: }
63:
64: public DLAED0()
65: {
66:
67:
68: #region Dependencies (Initialization)
69:
70: DCOPY dcopy = new DCOPY();
71: LSAME lsame = new LSAME();
72: XERBLA xerbla = new XERBLA();
73: IDAMAX idamax = new IDAMAX();
74: DLAMC3 dlamc3 = new DLAMC3();
75: DLAPY2 dlapy2 = new DLAPY2();
76: DLAMRG dlamrg = new DLAMRG();
77: DROT drot = new DROT();
78: DSCAL dscal = new DSCAL();
79: DNRM2 dnrm2 = new DNRM2();
80: DLAED5 dlaed5 = new DLAED5();
81: DLASSQ dlassq = new DLASSQ();
82: DLAE2 dlae2 = new DLAE2();
83: DLAEV2 dlaev2 = new DLAEV2();
84: DSWAP dswap = new DSWAP();
85: IEEECK ieeeck = new IEEECK();
86: IPARMQ iparmq = new IPARMQ();
87: DGEMM dgemm = new DGEMM(lsame, xerbla);
88: DLACPY dlacpy = new DLACPY(lsame);
89: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
90: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
91: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
92: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
93: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
94: DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
95: DLAED6 dlaed6 = new DLAED6(dlamch);
96: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
97: DLASET dlaset = new DLASET(lsame);
98: DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
99: DLAED1 dlaed1 = new DLAED1(dcopy, dlaed2, dlaed3, dlamrg, xerbla);
100: DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
101: DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
102: DGEMV dgemv = new DGEMV(lsame, xerbla);
103: DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
104: DLAED7 dlaed7 = new DLAED7(dgemm, dlaed8, dlaed9, dlaeda, dlamrg, xerbla);
105: DLANST dlanst = new DLANST(lsame, dlassq);
106: DLARTG dlartg = new DLARTG(dlamch);
107: DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
108: DLASR dlasr = new DLASR(lsame, xerbla);
109: DLASRT dlasrt = new DLASRT(lsame, xerbla);
110: DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
111: , dlasrt, dswap, xerbla);
112: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
113:
114: #endregion
115:
116:
117: #region Set Dependencies
118:
119: this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; this._dlaed1 = dlaed1; this._dlaed7 = dlaed7;
120: this._dsteqr = dsteqr;this._xerbla = xerbla; this._ilaenv = ilaenv;
121:
122: #endregion
123:
124: }
125: /// <summary>
126: /// Purpose
127: /// =======
128: ///
129: /// DLAED0 computes all eigenvalues and corresponding eigenvectors of a
130: /// symmetric tridiagonal matrix using the divide and conquer method.
131: ///
132: ///</summary>
133: /// <param name="ICOMPQ">
134: /// (input) INTEGER
135: /// = 0: Compute eigenvalues only.
136: /// = 1: Compute eigenvectors of original dense symmetric matrix
137: /// also. On entry, Q contains the orthogonal matrix used
138: /// to reduce the original matrix to tridiagonal form.
139: /// = 2: Compute eigenvalues and eigenvectors of tridiagonal
140: /// matrix.
141: ///</param>
142: /// <param name="QSIZ">
143: /// (input) INTEGER
144: /// The dimension of the orthogonal matrix used to reduce
145: /// the full matrix to tridiagonal form. QSIZ .GE. N if ICOMPQ = 1.
146: ///</param>
147: /// <param name="N">
148: /// (input) INTEGER
149: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
150: ///</param>
151: /// <param name="D">
152: /// (input/output) DOUBLE PRECISION array, dimension (N)
153: /// On entry, the main diagonal of the tridiagonal matrix.
154: /// On exit, its eigenvalues.
155: ///</param>
156: /// <param name="E">
157: /// (input) DOUBLE PRECISION array, dimension (N-1)
158: /// The off-diagonal elements of the tridiagonal matrix.
159: /// On exit, E has been destroyed.
160: ///</param>
161: /// <param name="Q">
162: /// (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
163: /// On entry, Q must contain an N-by-N orthogonal matrix.
164: /// If ICOMPQ = 0 Q is not referenced.
165: /// If ICOMPQ = 1 On entry, Q is a subset of the columns of the
166: /// orthogonal matrix used to reduce the full
167: /// matrix to tridiagonal form corresponding to
168: /// the subset of the full matrix which is being
169: /// decomposed at this time.
170: /// If ICOMPQ = 2 On entry, Q will be the identity matrix.
171: /// On exit, Q contains the eigenvectors of the
172: /// tridiagonal matrix.
173: ///</param>
174: /// <param name="LDQ">
175: /// (input) INTEGER
176: /// The leading dimension of the array Q. If eigenvectors are
177: /// desired, then LDQ .GE. max(1,N). In any case, LDQ .GE. 1.
178: ///</param>
179: /// <param name="QSTORE">
180: /// (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
181: /// Referenced only when ICOMPQ = 1. Used to store parts of
182: /// the eigenvector matrix when the updating matrix multiplies
183: /// take place.
184: ///</param>
185: /// <param name="LDQS">
186: /// (input) INTEGER
187: /// The leading dimension of the array QSTORE. If ICOMPQ = 1,
188: /// then LDQS .GE. max(1,N). In any case, LDQS .GE. 1.
189: ///</param>
190: /// <param name="WORK">
191: /// (workspace) DOUBLE PRECISION array,
192: /// If ICOMPQ = 0 or 1, the dimension of WORK must be at least
193: /// 1 + 3*N + 2*N*lg N + 2*N**2
194: /// ( lg( N ) = smallest integer k
195: /// such that 2^k .GE. N )
196: /// If ICOMPQ = 2, the dimension of WORK must be at least
197: /// 4*N + N**2.
198: ///</param>
199: /// <param name="IWORK">
200: /// (workspace) INTEGER array,
201: /// If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
202: /// 6 + 6*N + 5*N*lg N.
203: /// ( lg( N ) = smallest integer k
204: /// such that 2^k .GE. N )
205: /// If ICOMPQ = 2, the dimension of IWORK must be at least
206: /// 3 + 5*N.
207: ///</param>
208: /// <param name="INFO">
209: /// (output) INTEGER
210: /// = 0: successful exit.
211: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
212: /// .GT. 0: The algorithm failed to compute an eigenvalue while
213: /// working on the submatrix lying in rows and columns
214: /// INFO/(N+1) through mod(INFO,N+1).
215: ///</param>
216: public void Run(int ICOMPQ, int QSIZ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Q, int offset_q
217: , int LDQ, ref double[] QSTORE, int offset_qstore, int LDQS, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
218: {
219:
220: #region Array Index Correction
221:
222: int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_q = -1 - LDQ + offset_q;
223: int o_qstore = -1 - LDQS + offset_qstore; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
224:
225: #endregion
226:
227:
228: #region Prolog
229:
230: // *
231: // * -- LAPACK routine (version 3.1) --
232: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
233: // * November 2006
234: // *
235: // * .. Scalar Arguments ..
236: // * ..
237: // * .. Array Arguments ..
238: // * ..
239: // *
240: // * Purpose
241: // * =======
242: // *
243: // * DLAED0 computes all eigenvalues and corresponding eigenvectors of a
244: // * symmetric tridiagonal matrix using the divide and conquer method.
245: // *
246: // * Arguments
247: // * =========
248: // *
249: // * ICOMPQ (input) INTEGER
250: // * = 0: Compute eigenvalues only.
251: // * = 1: Compute eigenvectors of original dense symmetric matrix
252: // * also. On entry, Q contains the orthogonal matrix used
253: // * to reduce the original matrix to tridiagonal form.
254: // * = 2: Compute eigenvalues and eigenvectors of tridiagonal
255: // * matrix.
256: // *
257: // * QSIZ (input) INTEGER
258: // * The dimension of the orthogonal matrix used to reduce
259: // * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
260: // *
261: // * N (input) INTEGER
262: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
263: // *
264: // * D (input/output) DOUBLE PRECISION array, dimension (N)
265: // * On entry, the main diagonal of the tridiagonal matrix.
266: // * On exit, its eigenvalues.
267: // *
268: // * E (input) DOUBLE PRECISION array, dimension (N-1)
269: // * The off-diagonal elements of the tridiagonal matrix.
270: // * On exit, E has been destroyed.
271: // *
272: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
273: // * On entry, Q must contain an N-by-N orthogonal matrix.
274: // * If ICOMPQ = 0 Q is not referenced.
275: // * If ICOMPQ = 1 On entry, Q is a subset of the columns of the
276: // * orthogonal matrix used to reduce the full
277: // * matrix to tridiagonal form corresponding to
278: // * the subset of the full matrix which is being
279: // * decomposed at this time.
280: // * If ICOMPQ = 2 On entry, Q will be the identity matrix.
281: // * On exit, Q contains the eigenvectors of the
282: // * tridiagonal matrix.
283: // *
284: // * LDQ (input) INTEGER
285: // * The leading dimension of the array Q. If eigenvectors are
286: // * desired, then LDQ >= max(1,N). In any case, LDQ >= 1.
287: // *
288: // * QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
289: // * Referenced only when ICOMPQ = 1. Used to store parts of
290: // * the eigenvector matrix when the updating matrix multiplies
291: // * take place.
292: // *
293: // * LDQS (input) INTEGER
294: // * The leading dimension of the array QSTORE. If ICOMPQ = 1,
295: // * then LDQS >= max(1,N). In any case, LDQS >= 1.
296: // *
297: // * WORK (workspace) DOUBLE PRECISION array,
298: // * If ICOMPQ = 0 or 1, the dimension of WORK must be at least
299: // * 1 + 3*N + 2*N*lg N + 2*N**2
300: // * ( lg( N ) = smallest integer k
301: // * such that 2^k >= N )
302: // * If ICOMPQ = 2, the dimension of WORK must be at least
303: // * 4*N + N**2.
304: // *
305: // * IWORK (workspace) INTEGER array,
306: // * If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
307: // * 6 + 6*N + 5*N*lg N.
308: // * ( lg( N ) = smallest integer k
309: // * such that 2^k >= N )
310: // * If ICOMPQ = 2, the dimension of IWORK must be at least
311: // * 3 + 5*N.
312: // *
313: // * INFO (output) INTEGER
314: // * = 0: successful exit.
315: // * < 0: if INFO = -i, the i-th argument had an illegal value.
316: // * > 0: The algorithm failed to compute an eigenvalue while
317: // * working on the submatrix lying in rows and columns
318: // * INFO/(N+1) through mod(INFO,N+1).
319: // *
320: // * Further Details
321: // * ===============
322: // *
323: // * Based on contributions by
324: // * Jeff Rutter, Computer Science Division, University of California
325: // * at Berkeley, USA
326: // *
327: // * =====================================================================
328: // *
329: // * .. Parameters ..
330: // * ..
331: // * .. Local Scalars ..
332: // * ..
333: // * .. External Subroutines ..
334: // * ..
335: // * .. External Functions ..
336: // * ..
337: // * .. Intrinsic Functions ..
338: // INTRINSIC ABS, DBLE, INT, LOG, MAX;
339: // * ..
340: // * .. Executable Statements ..
341: // *
342: // * Test the input parameters.
343: // *
344:
345: #endregion
346:
347:
348: #region Body
349:
350: INFO = 0;
351: // *
352: if (ICOMPQ < 0 || ICOMPQ > 2)
353: {
354: INFO = - 1;
355: }
356: else
357: {
358: if ((ICOMPQ == 1) && (QSIZ < Math.Max(0, N)))
359: {
360: INFO = - 2;
361: }
362: else
363: {
364: if (N < 0)
365: {
366: INFO = - 3;
367: }
368: else
369: {
370: if (LDQ < Math.Max(1, N))
371: {
372: INFO = - 7;
373: }
374: else
375: {
376: if (LDQS < Math.Max(1, N))
377: {
378: INFO = - 9;
379: }
380: }
381: }
382: }
383: }
384: if (INFO != 0)
385: {
386: this._xerbla.Run("DLAED0", - INFO);
387: return;
388: }
389: // *
390: // * Quick return if possible
391: // *
392: if (N == 0) return;
393: // *
394: SMLSIZ = this._ilaenv.Run(9, "DLAED0", " ", 0, 0, 0, 0);
395: // *
396: // * Determine the size and placement of the submatrices, and save in
397: // * the leading elements of IWORK.
398: // *
399: IWORK[1 + o_iwork] = N;
400: SUBPBS = 1;
401: TLVLS = 0;
402: LABEL10:;
403: if (IWORK[SUBPBS + o_iwork] > SMLSIZ)
404: {
405: for (J = SUBPBS; J >= 1; J += - 1)
406: {
407: IWORK[2 * J + o_iwork] = (IWORK[J + o_iwork] + 1) / 2;
408: IWORK[2 * J - 1 + o_iwork] = IWORK[J + o_iwork] / 2;
409: }
410: TLVLS = TLVLS + 1;
411: SUBPBS = 2 * SUBPBS;
412: goto LABEL10;
413: }
414: for (J = 2; J <= SUBPBS; J++)
415: {
416: IWORK[J + o_iwork] = IWORK[J + o_iwork] + IWORK[J - 1 + o_iwork];
417: }
418: // *
419: // * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
420: // * using rank-1 modifications (cuts).
421: // *
422: SPM1 = SUBPBS - 1;
423: for (I = 1; I <= SPM1; I++)
424: {
425: SUBMAT = IWORK[I + o_iwork] + 1;
426: SMM1 = SUBMAT - 1;
427: D[SMM1 + o_d] = D[SMM1 + o_d] - Math.Abs(E[SMM1 + o_e]);
428: D[SUBMAT + o_d] = D[SUBMAT + o_d] - Math.Abs(E[SMM1 + o_e]);
429: }
430: // *
431: INDXQ = 4 * N + 3;
432: if (ICOMPQ != 2)
433: {
434: // *
435: // * Set up workspaces for eigenvalues only/accumulate new vectors
436: // * routine
437: // *
438: TEMP = Math.Log(Convert.ToDouble(N)) / Math.Log(TWO);
439: LGN = Convert.ToInt32(Math.Truncate(TEMP));
440: if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
441: if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
442: IPRMPT = INDXQ + N + 1;
443: IPERM = IPRMPT + N * LGN;
444: IQPTR = IPERM + N * LGN;
445: IGIVPT = IQPTR + N + 2;
446: IGIVCL = IGIVPT + N * LGN;
447: // *
448: IGIVNM = 1;
449: IQ = IGIVNM + 2 * N * LGN;
450: IWREM = IQ + (int)Math.Pow(N, 2) + 1;
451: // *
452: // * Initialize pointers
453: // *
454: for (I = 0; I <= SUBPBS; I++)
455: {
456: IWORK[IPRMPT + I + o_iwork] = 1;
457: IWORK[IGIVPT + I + o_iwork] = 1;
458: }
459: IWORK[IQPTR + o_iwork] = 1;
460: }
461: // *
462: // * Solve each submatrix eigenproblem at the bottom of the divide and
463: // * conquer tree.
464: // *
465: CURR = 0;
466: for (I = 0; I <= SPM1; I++)
467: {
468: if (I == 0)
469: {
470: SUBMAT = 1;
471: MATSIZ = IWORK[1 + o_iwork];
472: }
473: else
474: {
475: SUBMAT = IWORK[I + o_iwork] + 1;
476: MATSIZ = IWORK[I + 1 + o_iwork] - IWORK[I + o_iwork];
477: }
478: if (ICOMPQ == 2)
479: {
480: this._dsteqr.Run("I", MATSIZ, ref D, SUBMAT + o_d, ref E, SUBMAT + o_e, ref Q, SUBMAT+SUBMAT * LDQ + o_q, LDQ
481: , ref WORK, offset_work, ref INFO);
482: if (INFO != 0) goto LABEL130;
483: }
484: else
485: {
486: this._dsteqr.Run("I", MATSIZ, ref D, SUBMAT + o_d, ref E, SUBMAT + o_e, ref WORK, IQ - 1 + IWORK[IQPTR + CURR + o_iwork] + o_work, MATSIZ
487: , ref WORK, offset_work, ref INFO);
488: if (INFO != 0) goto LABEL130;
489: if (ICOMPQ == 1)
490: {
491: this._dgemm.Run("N", "N", QSIZ, MATSIZ, MATSIZ, ONE
492: , Q, 1+SUBMAT * LDQ + o_q, LDQ, WORK, IQ - 1 + IWORK[IQPTR + CURR + o_iwork] + o_work, MATSIZ, ZERO, ref QSTORE, 1+SUBMAT * LDQS + o_qstore
493: , LDQS);
494: }
495: IWORK[IQPTR + CURR + 1 + o_iwork] = IWORK[IQPTR + CURR + o_iwork] + (int)Math.Pow(MATSIZ, 2);
496: CURR = CURR + 1;
497: }
498: K = 1;
499: for (J = SUBMAT; J <= IWORK[I + 1 + o_iwork]; J++)
500: {
501: IWORK[INDXQ + J + o_iwork] = K;
502: K = K + 1;
503: }
504: }
505: // *
506: // * Successively merge eigensystems of adjacent submatrices
507: // * into eigensystem for the corresponding larger matrix.
508: // *
509: // * while ( SUBPBS > 1 )
510: // *
511: CURLVL = 1;
512: LABEL80:;
513: if (SUBPBS > 1)
514: {
515: SPM2 = SUBPBS - 2;
516: for (I = 0; I <= SPM2; I += 2)
517: {
518: if (I == 0)
519: {
520: SUBMAT = 1;
521: MATSIZ = IWORK[2 + o_iwork];
522: MSD2 = IWORK[1 + o_iwork];
523: CURPRB = 0;
524: }
525: else
526: {
527: SUBMAT = IWORK[I + o_iwork] + 1;
528: MATSIZ = IWORK[I + 2 + o_iwork] - IWORK[I + o_iwork];
529: MSD2 = MATSIZ / 2;
530: CURPRB = CURPRB + 1;
531: }
532: // *
533: // * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
534: // * into an eigensystem of size MATSIZ.
535: // * DLAED1 is used only for the full eigensystem of a tridiagonal
536: // * matrix.
537: // * DLAED7 handles the cases in which eigenvalues only or eigenvalues
538: // * and eigenvectors of a full symmetric matrix (which was reduced to
539: // * tridiagonal form) are desired.
540: // *
541: if (ICOMPQ == 2)
542: {
543: this._dlaed1.Run(MATSIZ, ref D, SUBMAT + o_d, ref Q, SUBMAT+SUBMAT * LDQ + o_q, LDQ, ref IWORK, INDXQ + SUBMAT + o_iwork, ref E[SUBMAT + MSD2 - 1 + o_e]
544: , MSD2, ref WORK, offset_work, ref IWORK, SUBPBS + 1 + o_iwork, ref INFO);
545: }
546: else
547: {
548: this._dlaed7.Run(ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB
549: , ref D, SUBMAT + o_d, ref QSTORE, 1+SUBMAT * LDQS + o_qstore, LDQS, ref IWORK, INDXQ + SUBMAT + o_iwork, ref E[SUBMAT + MSD2 - 1 + o_e], MSD2
550: , ref WORK, IQ + o_work, ref IWORK, IQPTR + o_iwork, ref IWORK, IPRMPT + o_iwork, ref IWORK, IPERM + o_iwork, ref IWORK, IGIVPT + o_iwork, ref IWORK, IGIVCL + o_iwork
551: , ref WORK, IGIVNM + o_work, ref WORK, IWREM + o_work, ref IWORK, SUBPBS + 1 + o_iwork, ref INFO);
552: }
553: if (INFO != 0) goto LABEL130;
554: IWORK[I / 2 + 1 + o_iwork] = IWORK[I + 2 + o_iwork];
555: }
556: SUBPBS = SUBPBS / 2;
557: CURLVL = CURLVL + 1;
558: goto LABEL80;
559: }
560: // *
561: // * end while
562: // *
563: // * Re-merge the eigenvalues/vectors which were deflated at the final
564: // * merge step.
565: // *
566: if (ICOMPQ == 1)
567: {
568: for (I = 1; I <= N; I++)
569: {
570: J = IWORK[INDXQ + I + o_iwork];
571: WORK[I + o_work] = D[J + o_d];
572: this._dcopy.Run(QSIZ, QSTORE, 1+J * LDQS + o_qstore, 1, ref Q, 1+I * LDQ + o_q, 1);
573: }
574: this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
575: }
576: else
577: {
578: if (ICOMPQ == 2)
579: {
580: for (I = 1; I <= N; I++)
581: {
582: J = IWORK[INDXQ + I + o_iwork];
583: WORK[I + o_work] = D[J + o_d];
584: this._dcopy.Run(N, Q, 1+J * LDQ + o_q, 1, ref WORK, N * I + 1 + o_work, 1);
585: }
586: this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
587: this._dlacpy.Run("A", N, N, WORK, N + 1 + o_work, N, ref Q, offset_q
588: , LDQ);
589: }
590: else
591: {
592: for (I = 1; I <= N; I++)
593: {
594: J = IWORK[INDXQ + I + o_iwork];
595: WORK[I + o_work] = D[J + o_d];
596: }
597: this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
598: }
599: }
600: goto LABEL140;
601: // *
602: LABEL130:;
603: INFO = SUBMAT * (N + 1) + SUBMAT + MATSIZ - 1;
604: // *
605: LABEL140:;
606: return;
607: // *
608: // * End of DLAED0
609: // *
610:
611: #endregion
612:
613: }
614: }
615: }