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: /// DLAED8 merges the two sets of eigenvalues together into a single
27: /// sorted set. Then it tries to deflate the size of the problem.
28: /// There are two ways in which deflation can occur: when two or more
29: /// eigenvalues are close together or if there is a tiny element in the
30: /// Z vector. For each such occurrence the order of the related secular
31: /// equation problem is reduced by one.
32: ///
33: ///</summary>
34: public class DLAED8
35: {
36:
37:
38: #region Dependencies
39:
40: IDAMAX _idamax; DLAMCH _dlamch; DLAPY2 _dlapy2; DCOPY _dcopy; DLACPY _dlacpy; DLAMRG _dlamrg; DROT _drot; DSCAL _dscal;
41: XERBLA _xerbla;
42:
43: #endregion
44:
45:
46: #region Fields
47:
48: const double MONE = - 1.0E0; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0;
49: const double EIGHT = 8.0E0;int I = 0; int IMAX = 0; int J = 0; int JLAM = 0; int JMAX = 0; int JP = 0; int K2 = 0;
50: int N1 = 0;int N1P1 = 0; int N2 = 0; double C = 0; double EPS = 0; double S = 0; double T = 0; double TAU = 0;
51: double TOL = 0;
52:
53: #endregion
54:
55: public DLAED8(IDAMAX idamax, DLAMCH dlamch, DLAPY2 dlapy2, DCOPY dcopy, DLACPY dlacpy, DLAMRG dlamrg, DROT drot, DSCAL dscal, XERBLA xerbla)
56: {
57:
58:
59: #region Set Dependencies
60:
61: this._idamax = idamax; this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy;
62: this._dlamrg = dlamrg;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla;
63:
64: #endregion
65:
66: }
67:
68: public DLAED8()
69: {
70:
71:
72: #region Dependencies (Initialization)
73:
74: IDAMAX idamax = new IDAMAX();
75: LSAME lsame = new LSAME();
76: DLAMC3 dlamc3 = new DLAMC3();
77: DLAPY2 dlapy2 = new DLAPY2();
78: DCOPY dcopy = new DCOPY();
79: DLAMRG dlamrg = new DLAMRG();
80: DROT drot = new DROT();
81: DSCAL dscal = new DSCAL();
82: XERBLA xerbla = new XERBLA();
83: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
84: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
85: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
86: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
87: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
88: DLACPY dlacpy = new DLACPY(lsame);
89:
90: #endregion
91:
92:
93: #region Set Dependencies
94:
95: this._idamax = idamax; this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy;
96: this._dlamrg = dlamrg;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla;
97:
98: #endregion
99:
100: }
101: /// <summary>
102: /// Purpose
103: /// =======
104: ///
105: /// DLAED8 merges the two sets of eigenvalues together into a single
106: /// sorted set. Then it tries to deflate the size of the problem.
107: /// There are two ways in which deflation can occur: when two or more
108: /// eigenvalues are close together or if there is a tiny element in the
109: /// Z vector. For each such occurrence the order of the related secular
110: /// equation problem is reduced by one.
111: ///
112: ///</summary>
113: /// <param name="ICOMPQ">
114: /// (input) INTEGER
115: /// = 0: Compute eigenvalues only.
116: /// = 1: Compute eigenvectors of original dense symmetric matrix
117: /// also. On entry, Q contains the orthogonal matrix used
118: /// to reduce the original matrix to tridiagonal form.
119: ///</param>
120: /// <param name="K">
121: /// (output) INTEGER
122: /// The number of non-deflated eigenvalues, and the order of the
123: /// related secular equation.
124: ///</param>
125: /// <param name="N">
126: /// (input) INTEGER
127: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
128: ///</param>
129: /// <param name="QSIZ">
130: /// (input) INTEGER
131: /// The dimension of the orthogonal matrix used to reduce
132: /// the full matrix to tridiagonal form. QSIZ .GE. N if ICOMPQ = 1.
133: ///</param>
134: /// <param name="D">
135: /// (input/output) DOUBLE PRECISION array, dimension (N)
136: /// On entry, the eigenvalues of the two submatrices to be
137: /// combined. On exit, the trailing (N-K) updated eigenvalues
138: /// (those which were deflated) sorted into increasing order.
139: ///</param>
140: /// <param name="Q">
141: /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
142: /// If ICOMPQ = 0, Q is not referenced. Otherwise,
143: /// on entry, Q contains the eigenvectors of the partially solved
144: /// system which has been previously updated in matrix
145: /// multiplies with other partially solved eigensystems.
146: /// On exit, Q contains the trailing (N-K) updated eigenvectors
147: /// (those which were deflated) in its last N-K columns.
148: ///</param>
149: /// <param name="LDQ">
150: /// (input) INTEGER
151: /// The leading dimension of the array Q. LDQ .GE. max(1,N).
152: ///</param>
153: /// <param name="INDXQ">
154: /// (input) INTEGER array, dimension (N)
155: /// The permutation which separately sorts the two sub-problems
156: /// in D into ascending order. Note that elements in the second
157: /// half of this permutation must first have CUTPNT added to
158: /// their values in order to be accurate.
159: ///</param>
160: /// <param name="RHO">
161: /// (input/output) DOUBLE PRECISION
162: /// On entry, the off-diagonal element associated with the rank-1
163: /// cut which originally split the two submatrices which are now
164: /// being recombined.
165: /// On exit, RHO has been modified to the value required by
166: /// DLAED3.
167: ///</param>
168: /// <param name="CUTPNT">
169: /// (input) INTEGER
170: /// The location of the last eigenvalue in the leading
171: /// sub-matrix. min(1,N) .LE. CUTPNT .LE. N.
172: ///</param>
173: /// <param name="Z">
174: /// (input) DOUBLE PRECISION array, dimension (N)
175: /// On entry, Z contains the updating vector (the last row of
176: /// the first sub-eigenvector matrix and the first row of the
177: /// second sub-eigenvector matrix).
178: /// On exit, the contents of Z are destroyed by the updating
179: /// process.
180: ///</param>
181: /// <param name="DLAMDA">
182: /// (output) DOUBLE PRECISION array, dimension (N)
183: /// A copy of the first K eigenvalues which will be used by
184: /// DLAED3 to form the secular equation.
185: ///</param>
186: /// <param name="Q2">
187: /// (output) DOUBLE PRECISION array, dimension (LDQ2,N)
188: /// If ICOMPQ = 0, Q2 is not referenced. Otherwise,
189: /// a copy of the first K eigenvectors which will be used by
190: /// DLAED7 in a matrix multiply (DGEMM) to update the new
191: /// eigenvectors.
192: ///</param>
193: /// <param name="LDQ2">
194: /// (input) INTEGER
195: /// The leading dimension of the array Q2. LDQ2 .GE. max(1,N).
196: ///</param>
197: /// <param name="W">
198: /// (output) DOUBLE PRECISION array, dimension (N)
199: /// The first k values of the final deflation-altered z-vector and
200: /// will be passed to DLAED3.
201: ///</param>
202: /// <param name="PERM">
203: /// (output) INTEGER array, dimension (N)
204: /// The permutations (from deflation and sorting) to be applied
205: /// to each eigenblock.
206: ///</param>
207: /// <param name="GIVPTR">
208: /// (output) INTEGER
209: /// The number of Givens rotations which took place in this
210: /// subproblem.
211: ///</param>
212: /// <param name="GIVCOL">
213: /// (output) INTEGER array, dimension (2, N)
214: /// Each pair of numbers indicates a pair of columns to take place
215: /// in a Givens rotation.
216: ///</param>
217: /// <param name="GIVNUM">
218: /// (output) DOUBLE PRECISION array, dimension (2, N)
219: /// Each number indicates the S value to be used in the
220: /// corresponding Givens rotation.
221: ///</param>
222: /// <param name="INDXP">
223: /// (workspace) INTEGER array, dimension (N)
224: /// The permutation used to place deflated values of D at the end
225: /// of the array. INDXP(1:K) points to the nondeflated D-values
226: /// and INDXP(K+1:N) points to the deflated eigenvalues.
227: ///</param>
228: /// <param name="INDX">
229: /// (workspace) INTEGER array, dimension (N)
230: /// The permutation used to sort the contents of D into ascending
231: /// order.
232: ///</param>
233: /// <param name="INFO">
234: /// (output) INTEGER
235: /// = 0: successful exit.
236: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
237: ///</param>
238: public void Run(int ICOMPQ, ref int K, int N, int QSIZ, ref double[] D, int offset_d, ref double[] Q, int offset_q
239: , int LDQ, ref int[] INDXQ, int offset_indxq, ref double RHO, int CUTPNT, ref double[] Z, int offset_z, ref double[] DLAMDA, int offset_dlamda
240: , ref double[] Q2, int offset_q2, int LDQ2, ref double[] W, int offset_w, ref int[] PERM, int offset_perm, ref int GIVPTR, ref int[] GIVCOL, int offset_givcol
241: , ref double[] GIVNUM, int offset_givnum, ref int[] INDXP, int offset_indxp, ref int[] INDX, int offset_indx, ref int INFO)
242: {
243:
244: #region Array Index Correction
245:
246: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_indxq = -1 + offset_indxq; int o_z = -1 + offset_z;
247: int o_dlamda = -1 + offset_dlamda; int o_q2 = -1 - LDQ2 + offset_q2; int o_w = -1 + offset_w;
248: int o_perm = -1 + offset_perm; int o_givcol = -3 + offset_givcol; int o_givnum = -3 + offset_givnum;
249: int o_indxp = -1 + offset_indxp; int o_indx = -1 + offset_indx;
250:
251: #endregion
252:
253:
254: #region Prolog
255:
256: // *
257: // * -- LAPACK routine (version 3.1) --
258: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
259: // * November 2006
260: // *
261: // * .. Scalar Arguments ..
262: // * ..
263: // * .. Array Arguments ..
264: // * ..
265: // *
266: // * Purpose
267: // * =======
268: // *
269: // * DLAED8 merges the two sets of eigenvalues together into a single
270: // * sorted set. Then it tries to deflate the size of the problem.
271: // * There are two ways in which deflation can occur: when two or more
272: // * eigenvalues are close together or if there is a tiny element in the
273: // * Z vector. For each such occurrence the order of the related secular
274: // * equation problem is reduced by one.
275: // *
276: // * Arguments
277: // * =========
278: // *
279: // * ICOMPQ (input) INTEGER
280: // * = 0: Compute eigenvalues only.
281: // * = 1: Compute eigenvectors of original dense symmetric matrix
282: // * also. On entry, Q contains the orthogonal matrix used
283: // * to reduce the original matrix to tridiagonal form.
284: // *
285: // * K (output) INTEGER
286: // * The number of non-deflated eigenvalues, and the order of the
287: // * related secular equation.
288: // *
289: // * N (input) INTEGER
290: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
291: // *
292: // * QSIZ (input) INTEGER
293: // * The dimension of the orthogonal matrix used to reduce
294: // * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
295: // *
296: // * D (input/output) DOUBLE PRECISION array, dimension (N)
297: // * On entry, the eigenvalues of the two submatrices to be
298: // * combined. On exit, the trailing (N-K) updated eigenvalues
299: // * (those which were deflated) sorted into increasing order.
300: // *
301: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
302: // * If ICOMPQ = 0, Q is not referenced. Otherwise,
303: // * on entry, Q contains the eigenvectors of the partially solved
304: // * system which has been previously updated in matrix
305: // * multiplies with other partially solved eigensystems.
306: // * On exit, Q contains the trailing (N-K) updated eigenvectors
307: // * (those which were deflated) in its last N-K columns.
308: // *
309: // * LDQ (input) INTEGER
310: // * The leading dimension of the array Q. LDQ >= max(1,N).
311: // *
312: // * INDXQ (input) INTEGER array, dimension (N)
313: // * The permutation which separately sorts the two sub-problems
314: // * in D into ascending order. Note that elements in the second
315: // * half of this permutation must first have CUTPNT added to
316: // * their values in order to be accurate.
317: // *
318: // * RHO (input/output) DOUBLE PRECISION
319: // * On entry, the off-diagonal element associated with the rank-1
320: // * cut which originally split the two submatrices which are now
321: // * being recombined.
322: // * On exit, RHO has been modified to the value required by
323: // * DLAED3.
324: // *
325: // * CUTPNT (input) INTEGER
326: // * The location of the last eigenvalue in the leading
327: // * sub-matrix. min(1,N) <= CUTPNT <= N.
328: // *
329: // * Z (input) DOUBLE PRECISION array, dimension (N)
330: // * On entry, Z contains the updating vector (the last row of
331: // * the first sub-eigenvector matrix and the first row of the
332: // * second sub-eigenvector matrix).
333: // * On exit, the contents of Z are destroyed by the updating
334: // * process.
335: // *
336: // * DLAMDA (output) DOUBLE PRECISION array, dimension (N)
337: // * A copy of the first K eigenvalues which will be used by
338: // * DLAED3 to form the secular equation.
339: // *
340: // * Q2 (output) DOUBLE PRECISION array, dimension (LDQ2,N)
341: // * If ICOMPQ = 0, Q2 is not referenced. Otherwise,
342: // * a copy of the first K eigenvectors which will be used by
343: // * DLAED7 in a matrix multiply (DGEMM) to update the new
344: // * eigenvectors.
345: // *
346: // * LDQ2 (input) INTEGER
347: // * The leading dimension of the array Q2. LDQ2 >= max(1,N).
348: // *
349: // * W (output) DOUBLE PRECISION array, dimension (N)
350: // * The first k values of the final deflation-altered z-vector and
351: // * will be passed to DLAED3.
352: // *
353: // * PERM (output) INTEGER array, dimension (N)
354: // * The permutations (from deflation and sorting) to be applied
355: // * to each eigenblock.
356: // *
357: // * GIVPTR (output) INTEGER
358: // * The number of Givens rotations which took place in this
359: // * subproblem.
360: // *
361: // * GIVCOL (output) INTEGER array, dimension (2, N)
362: // * Each pair of numbers indicates a pair of columns to take place
363: // * in a Givens rotation.
364: // *
365: // * GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
366: // * Each number indicates the S value to be used in the
367: // * corresponding Givens rotation.
368: // *
369: // * INDXP (workspace) INTEGER array, dimension (N)
370: // * The permutation used to place deflated values of D at the end
371: // * of the array. INDXP(1:K) points to the nondeflated D-values
372: // * and INDXP(K+1:N) points to the deflated eigenvalues.
373: // *
374: // * INDX (workspace) INTEGER array, dimension (N)
375: // * The permutation used to sort the contents of D into ascending
376: // * order.
377: // *
378: // * INFO (output) INTEGER
379: // * = 0: successful exit.
380: // * < 0: if INFO = -i, the i-th argument had an illegal value.
381: // *
382: // * Further Details
383: // * ===============
384: // *
385: // * Based on contributions by
386: // * Jeff Rutter, Computer Science Division, University of California
387: // * at Berkeley, USA
388: // *
389: // * =====================================================================
390: // *
391: // * .. Parameters ..
392: // * ..
393: // * .. Local Scalars ..
394: // *
395: // * ..
396: // * .. External Functions ..
397: // * ..
398: // * .. External Subroutines ..
399: // * ..
400: // * .. Intrinsic Functions ..
401: // INTRINSIC ABS, MAX, MIN, SQRT;
402: // * ..
403: // * .. Executable Statements ..
404: // *
405: // * Test the input parameters.
406: // *
407:
408: #endregion
409:
410:
411: #region Body
412:
413: INFO = 0;
414: // *
415: if (ICOMPQ < 0 || ICOMPQ > 1)
416: {
417: INFO = - 1;
418: }
419: else
420: {
421: if (N < 0)
422: {
423: INFO = - 3;
424: }
425: else
426: {
427: if (ICOMPQ == 1 && QSIZ < N)
428: {
429: INFO = - 4;
430: }
431: else
432: {
433: if (LDQ < Math.Max(1, N))
434: {
435: INFO = - 7;
436: }
437: else
438: {
439: if (CUTPNT < Math.Min(1, N) || CUTPNT > N)
440: {
441: INFO = - 10;
442: }
443: else
444: {
445: if (LDQ2 < Math.Max(1, N))
446: {
447: INFO = - 14;
448: }
449: }
450: }
451: }
452: }
453: }
454: if (INFO != 0)
455: {
456: this._xerbla.Run("DLAED8", - INFO);
457: return;
458: }
459: // *
460: // * Quick return if possible
461: // *
462: if (N == 0) return;
463: // *
464: N1 = CUTPNT;
465: N2 = N - N1;
466: N1P1 = N1 + 1;
467: // *
468: if (RHO < ZERO)
469: {
470: this._dscal.Run(N2, MONE, ref Z, N1P1 + o_z, 1);
471: }
472: // *
473: // * Normalize z so that norm(z) = 1
474: // *
475: T = ONE / Math.Sqrt(TWO);
476: for (J = 1; J <= N; J++)
477: {
478: INDX[J + o_indx] = J;
479: }
480: this._dscal.Run(N, T, ref Z, offset_z, 1);
481: RHO = Math.Abs(TWO * RHO);
482: // *
483: // * Sort the eigenvalues into increasing order
484: // *
485: for (I = CUTPNT + 1; I <= N; I++)
486: {
487: INDXQ[I + o_indxq] = INDXQ[I + o_indxq] + CUTPNT;
488: }
489: for (I = 1; I <= N; I++)
490: {
491: DLAMDA[I + o_dlamda] = D[INDXQ[I + o_indxq] + o_d];
492: W[I + o_w] = Z[INDXQ[I + o_indxq] + o_z];
493: }
494: I = 1;
495: J = CUTPNT + 1;
496: this._dlamrg.Run(N1, N2, DLAMDA, offset_dlamda, 1, 1, ref INDX, offset_indx);
497: for (I = 1; I <= N; I++)
498: {
499: D[I + o_d] = DLAMDA[INDX[I + o_indx] + o_dlamda];
500: Z[I + o_z] = W[INDX[I + o_indx] + o_w];
501: }
502: // *
503: // * Calculate the allowable deflation tolerence
504: // *
505: IMAX = this._idamax.Run(N, Z, offset_z, 1);
506: JMAX = this._idamax.Run(N, D, offset_d, 1);
507: EPS = this._dlamch.Run("Epsilon");
508: TOL = EIGHT * EPS * Math.Abs(D[JMAX + o_d]);
509: // *
510: // * If the rank-1 modifier is small enough, no more needs to be done
511: // * except to reorganize Q so that its columns correspond with the
512: // * elements in D.
513: // *
514: if (RHO * Math.Abs(Z[IMAX + o_z]) <= TOL)
515: {
516: K = 0;
517: if (ICOMPQ == 0)
518: {
519: for (J = 1; J <= N; J++)
520: {
521: PERM[J + o_perm] = INDXQ[INDX[J + o_indx] + o_indxq];
522: }
523: }
524: else
525: {
526: for (J = 1; J <= N; J++)
527: {
528: PERM[J + o_perm] = INDXQ[INDX[J + o_indx] + o_indxq];
529: this._dcopy.Run(QSIZ, Q, 1+(PERM[J + o_perm]) * LDQ + o_q, 1, ref Q2, 1+J * LDQ2 + o_q2, 1);
530: }
531: this._dlacpy.Run("A", QSIZ, N, Q2, 1+1 * LDQ2 + o_q2, LDQ2, ref Q, 1+1 * LDQ + o_q
532: , LDQ);
533: }
534: return;
535: }
536: // *
537: // * If there are multiple eigenvalues then the problem deflates. Here
538: // * the number of equal eigenvalues are found. As each equal
539: // * eigenvalue is found, an elementary reflector is computed to rotate
540: // * the corresponding eigensubspace so that the corresponding
541: // * components of Z are zero in this new basis.
542: // *
543: K = 0;
544: GIVPTR = 0;
545: K2 = N + 1;
546: for (J = 1; J <= N; J++)
547: {
548: if (RHO * Math.Abs(Z[J + o_z]) <= TOL)
549: {
550: // *
551: // * Deflate due to small z component.
552: // *
553: K2 = K2 - 1;
554: INDXP[K2 + o_indxp] = J;
555: if (J == N) goto LABEL110;
556: }
557: else
558: {
559: JLAM = J;
560: goto LABEL80;
561: }
562: }
563: LABEL80:;
564: J = J + 1;
565: if (J > N) goto LABEL100;
566: if (RHO * Math.Abs(Z[J + o_z]) <= TOL)
567: {
568: // *
569: // * Deflate due to small z component.
570: // *
571: K2 = K2 - 1;
572: INDXP[K2 + o_indxp] = J;
573: }
574: else
575: {
576: // *
577: // * Check if eigenvalues are close enough to allow deflation.
578: // *
579: S = Z[JLAM + o_z];
580: C = Z[J + o_z];
581: // *
582: // * Find sqrt(a**2+b**2) without overflow or
583: // * destructive underflow.
584: // *
585: TAU = this._dlapy2.Run(C, S);
586: T = D[J + o_d] - D[JLAM + o_d];
587: C = C / TAU;
588: S = - S / TAU;
589: if (Math.Abs(T * C * S) <= TOL)
590: {
591: // *
592: // * Deflation is possible.
593: // *
594: Z[J + o_z] = TAU;
595: Z[JLAM + o_z] = ZERO;
596: // *
597: // * Record the appropriate Givens rotation
598: // *
599: GIVPTR = GIVPTR + 1;
600: GIVCOL[1+GIVPTR * 2 + o_givcol] = INDXQ[INDX[JLAM + o_indx] + o_indxq];
601: GIVCOL[2+GIVPTR * 2 + o_givcol] = INDXQ[INDX[J + o_indx] + o_indxq];
602: GIVNUM[1+GIVPTR * 2 + o_givnum] = C;
603: GIVNUM[2+GIVPTR * 2 + o_givnum] = S;
604: if (ICOMPQ == 1)
605: {
606: this._drot.Run(QSIZ, ref Q, 1+(INDXQ[INDX[JLAM + o_indx] + o_indxq]) * LDQ + o_q, 1, ref Q, 1+(INDXQ[INDX[J + o_indx] + o_indxq]) * LDQ + o_q, 1, C
607: , S);
608: }
609: T = D[JLAM + o_d] * C * C + D[J + o_d] * S * S;
610: D[J + o_d] = D[JLAM + o_d] * S * S + D[J + o_d] * C * C;
611: D[JLAM + o_d] = T;
612: K2 = K2 - 1;
613: I = 1;
614: LABEL90:;
615: if (K2 + I <= N)
616: {
617: if (D[JLAM + o_d] < D[INDXP[K2 + I + o_indxp] + o_d])
618: {
619: INDXP[K2 + I - 1 + o_indxp] = INDXP[K2 + I + o_indxp];
620: INDXP[K2 + I + o_indxp] = JLAM;
621: I = I + 1;
622: goto LABEL90;
623: }
624: else
625: {
626: INDXP[K2 + I - 1 + o_indxp] = JLAM;
627: }
628: }
629: else
630: {
631: INDXP[K2 + I - 1 + o_indxp] = JLAM;
632: }
633: JLAM = J;
634: }
635: else
636: {
637: K = K + 1;
638: W[K + o_w] = Z[JLAM + o_z];
639: DLAMDA[K + o_dlamda] = D[JLAM + o_d];
640: INDXP[K + o_indxp] = JLAM;
641: JLAM = J;
642: }
643: }
644: goto LABEL80;
645: LABEL100:;
646: // *
647: // * Record the last eigenvalue.
648: // *
649: K = K + 1;
650: W[K + o_w] = Z[JLAM + o_z];
651: DLAMDA[K + o_dlamda] = D[JLAM + o_d];
652: INDXP[K + o_indxp] = JLAM;
653: // *
654: LABEL110:;
655: // *
656: // * Sort the eigenvalues and corresponding eigenvectors into DLAMDA
657: // * and Q2 respectively. The eigenvalues/vectors which were not
658: // * deflated go into the first K slots of DLAMDA and Q2 respectively,
659: // * while those which were deflated go into the last N - K slots.
660: // *
661: if (ICOMPQ == 0)
662: {
663: for (J = 1; J <= N; J++)
664: {
665: JP = INDXP[J + o_indxp];
666: DLAMDA[J + o_dlamda] = D[JP + o_d];
667: PERM[J + o_perm] = INDXQ[INDX[JP + o_indx] + o_indxq];
668: }
669: }
670: else
671: {
672: for (J = 1; J <= N; J++)
673: {
674: JP = INDXP[J + o_indxp];
675: DLAMDA[J + o_dlamda] = D[JP + o_d];
676: PERM[J + o_perm] = INDXQ[INDX[JP + o_indx] + o_indxq];
677: this._dcopy.Run(QSIZ, Q, 1+(PERM[J + o_perm]) * LDQ + o_q, 1, ref Q2, 1+J * LDQ2 + o_q2, 1);
678: }
679: }
680: // *
681: // * The deflated eigenvalues and their corresponding vectors go back
682: // * into the last N - K slots of D and Q respectively.
683: // *
684: if (K < N)
685: {
686: if (ICOMPQ == 0)
687: {
688: this._dcopy.Run(N - K, DLAMDA, K + 1 + o_dlamda, 1, ref D, K + 1 + o_d, 1);
689: }
690: else
691: {
692: this._dcopy.Run(N - K, DLAMDA, K + 1 + o_dlamda, 1, ref D, K + 1 + o_d, 1);
693: this._dlacpy.Run("A", QSIZ, N - K, Q2, 1+(K + 1) * LDQ2 + o_q2, LDQ2, ref Q, 1+(K + 1) * LDQ + o_q
694: , LDQ);
695: }
696: }
697: // *
698: return;
699: // *
700: // * End of DLAED8
701: // *
702:
703: #endregion
704:
705: }
706: }
707: }