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: /// DLAED7 computes the updated eigensystem of a diagonal
27: /// matrix after modification by a rank-one symmetric matrix. This
28: /// routine is used only for the eigenproblem which requires all
29: /// eigenvalues and optionally eigenvectors of a dense symmetric matrix
30: /// that has been reduced to tridiagonal form. DLAED1 handles
31: /// the case in which all eigenvalues and eigenvectors of a symmetric
32: /// tridiagonal matrix are desired.
33: ///
34: /// T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
35: ///
36: /// where Z = Q'u, u is a vector of length N with ones in the
37: /// CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
38: ///
39: /// The eigenvectors of the original matrix are stored in Q, and the
40: /// eigenvalues are in D. The algorithm consists of three stages:
41: ///
42: /// The first stage consists of deflating the size of the problem
43: /// when there are multiple eigenvalues or if there is a zero in
44: /// the Z vector. For each such occurence the dimension of the
45: /// secular equation problem is reduced by one. This stage is
46: /// performed by the routine DLAED8.
47: ///
48: /// The second stage consists of calculating the updated
49: /// eigenvalues. This is done by finding the roots of the secular
50: /// equation via the routine DLAED4 (as called by DLAED9).
51: /// This routine also calculates the eigenvectors of the current
52: /// problem.
53: ///
54: /// The final stage consists of computing the updated eigenvectors
55: /// directly using the updated eigenvalues. The eigenvectors for
56: /// the current problem are multiplied with the eigenvectors from
57: /// the overall problem.
58: ///
59: ///</summary>
60: public class DLAED7
61: {
62:
63:
64: #region Dependencies
65:
66: DGEMM _dgemm; DLAED8 _dlaed8; DLAED9 _dlaed9; DLAEDA _dlaeda; DLAMRG _dlamrg; XERBLA _xerbla;
67:
68: #endregion
69:
70:
71: #region Fields
72:
73: const double ONE = 1.0E0; const double ZERO = 0.0E0; int COLTYP = 0; int CURR = 0; int I = 0; int IDLMDA = 0;
74: int INDX = 0;int INDXC = 0; int INDXP = 0; int IQ2 = 0; int IS = 0; int IW = 0; int IZ = 0; int K = 0; int LDQ2 = 0;
75: int N1 = 0;int N2 = 0; int PTR = 0;
76:
77: #endregion
78:
79: public DLAED7(DGEMM dgemm, DLAED8 dlaed8, DLAED9 dlaed9, DLAEDA dlaeda, DLAMRG dlamrg, XERBLA xerbla)
80: {
81:
82:
83: #region Set Dependencies
84:
85: this._dgemm = dgemm; this._dlaed8 = dlaed8; this._dlaed9 = dlaed9; this._dlaeda = dlaeda; this._dlamrg = dlamrg;
86: this._xerbla = xerbla;
87:
88: #endregion
89:
90: }
91:
92: public DLAED7()
93: {
94:
95:
96: #region Dependencies (Initialization)
97:
98: LSAME lsame = new LSAME();
99: XERBLA xerbla = new XERBLA();
100: IDAMAX idamax = new IDAMAX();
101: DLAMC3 dlamc3 = new DLAMC3();
102: DLAPY2 dlapy2 = new DLAPY2();
103: DCOPY dcopy = new DCOPY();
104: DLAMRG dlamrg = new DLAMRG();
105: DROT drot = new DROT();
106: DSCAL dscal = new DSCAL();
107: DNRM2 dnrm2 = new DNRM2();
108: DLAED5 dlaed5 = new DLAED5();
109: DGEMM dgemm = new DGEMM(lsame, xerbla);
110: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
111: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
112: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
113: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
114: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
115: DLACPY dlacpy = new DLACPY(lsame);
116: DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
117: DLAED6 dlaed6 = new DLAED6(dlamch);
118: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
119: DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
120: DGEMV dgemv = new DGEMV(lsame, xerbla);
121: DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
122:
123: #endregion
124:
125:
126: #region Set Dependencies
127:
128: this._dgemm = dgemm; this._dlaed8 = dlaed8; this._dlaed9 = dlaed9; this._dlaeda = dlaeda; this._dlamrg = dlamrg;
129: this._xerbla = xerbla;
130:
131: #endregion
132:
133: }
134: /// <summary>
135: /// Purpose
136: /// =======
137: ///
138: /// DLAED7 computes the updated eigensystem of a diagonal
139: /// matrix after modification by a rank-one symmetric matrix. This
140: /// routine is used only for the eigenproblem which requires all
141: /// eigenvalues and optionally eigenvectors of a dense symmetric matrix
142: /// that has been reduced to tridiagonal form. DLAED1 handles
143: /// the case in which all eigenvalues and eigenvectors of a symmetric
144: /// tridiagonal matrix are desired.
145: ///
146: /// T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
147: ///
148: /// where Z = Q'u, u is a vector of length N with ones in the
149: /// CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
150: ///
151: /// The eigenvectors of the original matrix are stored in Q, and the
152: /// eigenvalues are in D. The algorithm consists of three stages:
153: ///
154: /// The first stage consists of deflating the size of the problem
155: /// when there are multiple eigenvalues or if there is a zero in
156: /// the Z vector. For each such occurence the dimension of the
157: /// secular equation problem is reduced by one. This stage is
158: /// performed by the routine DLAED8.
159: ///
160: /// The second stage consists of calculating the updated
161: /// eigenvalues. This is done by finding the roots of the secular
162: /// equation via the routine DLAED4 (as called by DLAED9).
163: /// This routine also calculates the eigenvectors of the current
164: /// problem.
165: ///
166: /// The final stage consists of computing the updated eigenvectors
167: /// directly using the updated eigenvalues. The eigenvectors for
168: /// the current problem are multiplied with the eigenvectors from
169: /// the overall problem.
170: ///
171: ///</summary>
172: /// <param name="ICOMPQ">
173: /// (input) INTEGER
174: /// = 0: Compute eigenvalues only.
175: /// = 1: Compute eigenvectors of original dense symmetric matrix
176: /// also. On entry, Q contains the orthogonal matrix used
177: /// to reduce the original matrix to tridiagonal form.
178: ///</param>
179: /// <param name="N">
180: /// (input) INTEGER
181: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
182: ///</param>
183: /// <param name="QSIZ">
184: /// (input) INTEGER
185: /// The dimension of the orthogonal matrix used to reduce
186: /// the full matrix to tridiagonal form. QSIZ .GE. N if ICOMPQ = 1.
187: ///</param>
188: /// <param name="TLVLS">
189: /// (input) INTEGER
190: /// The total number of merging levels in the overall divide and
191: /// conquer tree.
192: ///</param>
193: /// <param name="CURLVL">
194: /// (input) INTEGER
195: /// The current level in the overall merge routine,
196: /// 0 .LE. CURLVL .LE. TLVLS.
197: ///</param>
198: /// <param name="CURPBM">
199: /// (input) INTEGER
200: /// The current problem in the current level in the overall
201: /// merge routine (counting from upper left to lower right).
202: ///</param>
203: /// <param name="D">
204: /// (input/output) DOUBLE PRECISION array, dimension (N)
205: /// On entry, the eigenvalues of the rank-1-perturbed matrix.
206: /// On exit, the eigenvalues of the repaired matrix.
207: ///</param>
208: /// <param name="Q">
209: /// (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
210: /// On entry, the eigenvectors of the rank-1-perturbed matrix.
211: /// On exit, the eigenvectors of the repaired tridiagonal matrix.
212: ///</param>
213: /// <param name="LDQ">
214: /// (input) INTEGER
215: /// The leading dimension of the array Q. LDQ .GE. max(1,N).
216: ///</param>
217: /// <param name="INDXQ">
218: /// (output) INTEGER array, dimension (N)
219: /// The permutation which will reintegrate the subproblem just
220: /// solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
221: /// will be in ascending order.
222: ///</param>
223: /// <param name="RHO">
224: /// (input) DOUBLE PRECISION
225: /// The subdiagonal element used to create the rank-1
226: /// modification.
227: ///</param>
228: /// <param name="CUTPNT">
229: /// (input) INTEGER
230: /// Contains the location of the last eigenvalue in the leading
231: /// sub-matrix. min(1,N) .LE. CUTPNT .LE. N.
232: ///</param>
233: /// <param name="QSTORE">
234: /// (input/output) DOUBLE PRECISION array, dimension (N**2+1)
235: /// Stores eigenvectors of submatrices encountered during
236: /// divide and conquer, packed together. QPTR points to
237: /// beginning of the submatrices.
238: ///</param>
239: /// <param name="QPTR">
240: /// (input/output) INTEGER array, dimension (N+2)
241: /// List of indices pointing to beginning of submatrices stored
242: /// in QSTORE. The submatrices are numbered starting at the
243: /// bottom left of the divide and conquer tree, from left to
244: /// right and bottom to top.
245: ///</param>
246: /// <param name="PRMPTR">
247: /// (input) INTEGER array, dimension (N lg N)
248: /// Contains a list of pointers which indicate where in PERM a
249: /// level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
250: /// indicates the size of the permutation and also the size of
251: /// the full, non-deflated problem.
252: ///</param>
253: /// <param name="PERM">
254: /// (input) INTEGER array, dimension (N lg N)
255: /// Contains the permutations (from deflation and sorting) to be
256: /// applied to each eigenblock.
257: ///</param>
258: /// <param name="GIVPTR">
259: /// (input) INTEGER array, dimension (N lg N)
260: /// Contains a list of pointers which indicate where in GIVCOL a
261: /// level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
262: /// indicates the number of Givens rotations.
263: ///</param>
264: /// <param name="GIVCOL">
265: /// (input) INTEGER array, dimension (2, N lg N)
266: /// Each pair of numbers indicates a pair of columns to take place
267: /// in a Givens rotation.
268: ///</param>
269: /// <param name="GIVNUM">
270: /// (input) DOUBLE PRECISION array, dimension (2, N lg N)
271: /// Each number indicates the S value to be used in the
272: /// corresponding Givens rotation.
273: ///</param>
274: /// <param name="WORK">
275: /// (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N)
276: ///</param>
277: /// <param name="IWORK">
278: /// (workspace) INTEGER array, dimension (4*N)
279: ///</param>
280: /// <param name="INFO">
281: /// (output) INTEGER
282: /// = 0: successful exit.
283: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
284: /// .GT. 0: if INFO = 1, an eigenvalue did not converge
285: ///</param>
286: public void Run(int ICOMPQ, int N, int QSIZ, int TLVLS, int CURLVL, int CURPBM
287: , ref double[] D, int offset_d, ref double[] Q, int offset_q, int LDQ, ref int[] INDXQ, int offset_indxq, ref double RHO, int CUTPNT
288: , ref double[] QSTORE, int offset_qstore, ref int[] QPTR, int offset_qptr, ref int[] PRMPTR, int offset_prmptr, ref int[] PERM, int offset_perm, ref int[] GIVPTR, int offset_givptr, ref int[] GIVCOL, int offset_givcol
289: , ref double[] GIVNUM, int offset_givnum, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
290: {
291:
292: #region Array Index Correction
293:
294: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_indxq = -1 + offset_indxq;
295: int o_qstore = -1 + offset_qstore; int o_qptr = -1 + offset_qptr; int o_prmptr = -1 + offset_prmptr;
296: int o_perm = -1 + offset_perm; int o_givptr = -1 + offset_givptr; int o_givcol = -3 + offset_givcol;
297: int o_givnum = -3 + offset_givnum; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
298:
299: #endregion
300:
301:
302: #region Prolog
303:
304: // *
305: // * -- LAPACK routine (version 3.1) --
306: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
307: // * November 2006
308: // *
309: // * .. Scalar Arguments ..
310: // * ..
311: // * .. Array Arguments ..
312: // * ..
313: // *
314: // * Purpose
315: // * =======
316: // *
317: // * DLAED7 computes the updated eigensystem of a diagonal
318: // * matrix after modification by a rank-one symmetric matrix. This
319: // * routine is used only for the eigenproblem which requires all
320: // * eigenvalues and optionally eigenvectors of a dense symmetric matrix
321: // * that has been reduced to tridiagonal form. DLAED1 handles
322: // * the case in which all eigenvalues and eigenvectors of a symmetric
323: // * tridiagonal matrix are desired.
324: // *
325: // * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
326: // *
327: // * where Z = Q'u, u is a vector of length N with ones in the
328: // * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
329: // *
330: // * The eigenvectors of the original matrix are stored in Q, and the
331: // * eigenvalues are in D. The algorithm consists of three stages:
332: // *
333: // * The first stage consists of deflating the size of the problem
334: // * when there are multiple eigenvalues or if there is a zero in
335: // * the Z vector. For each such occurence the dimension of the
336: // * secular equation problem is reduced by one. This stage is
337: // * performed by the routine DLAED8.
338: // *
339: // * The second stage consists of calculating the updated
340: // * eigenvalues. This is done by finding the roots of the secular
341: // * equation via the routine DLAED4 (as called by DLAED9).
342: // * This routine also calculates the eigenvectors of the current
343: // * problem.
344: // *
345: // * The final stage consists of computing the updated eigenvectors
346: // * directly using the updated eigenvalues. The eigenvectors for
347: // * the current problem are multiplied with the eigenvectors from
348: // * the overall problem.
349: // *
350: // * Arguments
351: // * =========
352: // *
353: // * ICOMPQ (input) INTEGER
354: // * = 0: Compute eigenvalues only.
355: // * = 1: Compute eigenvectors of original dense symmetric matrix
356: // * also. On entry, Q contains the orthogonal matrix used
357: // * to reduce the original matrix to tridiagonal form.
358: // *
359: // * N (input) INTEGER
360: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
361: // *
362: // * QSIZ (input) INTEGER
363: // * The dimension of the orthogonal matrix used to reduce
364: // * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1.
365: // *
366: // * TLVLS (input) INTEGER
367: // * The total number of merging levels in the overall divide and
368: // * conquer tree.
369: // *
370: // * CURLVL (input) INTEGER
371: // * The current level in the overall merge routine,
372: // * 0 <= CURLVL <= TLVLS.
373: // *
374: // * CURPBM (input) INTEGER
375: // * The current problem in the current level in the overall
376: // * merge routine (counting from upper left to lower right).
377: // *
378: // * D (input/output) DOUBLE PRECISION array, dimension (N)
379: // * On entry, the eigenvalues of the rank-1-perturbed matrix.
380: // * On exit, the eigenvalues of the repaired matrix.
381: // *
382: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
383: // * On entry, the eigenvectors of the rank-1-perturbed matrix.
384: // * On exit, the eigenvectors of the repaired tridiagonal matrix.
385: // *
386: // * LDQ (input) INTEGER
387: // * The leading dimension of the array Q. LDQ >= max(1,N).
388: // *
389: // * INDXQ (output) INTEGER array, dimension (N)
390: // * The permutation which will reintegrate the subproblem just
391: // * solved back into sorted order, i.e., D( INDXQ( I = 1, N ) )
392: // * will be in ascending order.
393: // *
394: // * RHO (input) DOUBLE PRECISION
395: // * The subdiagonal element used to create the rank-1
396: // * modification.
397: // *
398: // * CUTPNT (input) INTEGER
399: // * Contains the location of the last eigenvalue in the leading
400: // * sub-matrix. min(1,N) <= CUTPNT <= N.
401: // *
402: // * QSTORE (input/output) DOUBLE PRECISION array, dimension (N**2+1)
403: // * Stores eigenvectors of submatrices encountered during
404: // * divide and conquer, packed together. QPTR points to
405: // * beginning of the submatrices.
406: // *
407: // * QPTR (input/output) INTEGER array, dimension (N+2)
408: // * List of indices pointing to beginning of submatrices stored
409: // * in QSTORE. The submatrices are numbered starting at the
410: // * bottom left of the divide and conquer tree, from left to
411: // * right and bottom to top.
412: // *
413: // * PRMPTR (input) INTEGER array, dimension (N lg N)
414: // * Contains a list of pointers which indicate where in PERM a
415: // * level's permutation is stored. PRMPTR(i+1) - PRMPTR(i)
416: // * indicates the size of the permutation and also the size of
417: // * the full, non-deflated problem.
418: // *
419: // * PERM (input) INTEGER array, dimension (N lg N)
420: // * Contains the permutations (from deflation and sorting) to be
421: // * applied to each eigenblock.
422: // *
423: // * GIVPTR (input) INTEGER array, dimension (N lg N)
424: // * Contains a list of pointers which indicate where in GIVCOL a
425: // * level's Givens rotations are stored. GIVPTR(i+1) - GIVPTR(i)
426: // * indicates the number of Givens rotations.
427: // *
428: // * GIVCOL (input) INTEGER array, dimension (2, N lg N)
429: // * Each pair of numbers indicates a pair of columns to take place
430: // * in a Givens rotation.
431: // *
432: // * GIVNUM (input) DOUBLE PRECISION array, dimension (2, N lg N)
433: // * Each number indicates the S value to be used in the
434: // * corresponding Givens rotation.
435: // *
436: // * WORK (workspace) DOUBLE PRECISION array, dimension (3*N+QSIZ*N)
437: // *
438: // * IWORK (workspace) INTEGER array, dimension (4*N)
439: // *
440: // * INFO (output) INTEGER
441: // * = 0: successful exit.
442: // * < 0: if INFO = -i, the i-th argument had an illegal value.
443: // * > 0: if INFO = 1, an eigenvalue did not converge
444: // *
445: // * Further Details
446: // * ===============
447: // *
448: // * Based on contributions by
449: // * Jeff Rutter, Computer Science Division, University of California
450: // * at Berkeley, USA
451: // *
452: // * =====================================================================
453: // *
454: // * .. Parameters ..
455: // * ..
456: // * .. Local Scalars ..
457: // * ..
458: // * .. External Subroutines ..
459: // * ..
460: // * .. Intrinsic Functions ..
461: // INTRINSIC MAX, MIN;
462: // * ..
463: // * .. Executable Statements ..
464: // *
465: // * Test the input parameters.
466: // *
467:
468: #endregion
469:
470:
471: #region Body
472:
473: INFO = 0;
474: // *
475: if (ICOMPQ < 0 || ICOMPQ > 1)
476: {
477: INFO = - 1;
478: }
479: else
480: {
481: if (N < 0)
482: {
483: INFO = - 2;
484: }
485: else
486: {
487: if (ICOMPQ == 1 && QSIZ < N)
488: {
489: INFO = - 4;
490: }
491: else
492: {
493: if (LDQ < Math.Max(1, N))
494: {
495: INFO = - 9;
496: }
497: else
498: {
499: if (Math.Min(1, N) > CUTPNT || N < CUTPNT)
500: {
501: INFO = - 12;
502: }
503: }
504: }
505: }
506: }
507: if (INFO != 0)
508: {
509: this._xerbla.Run("DLAED7", - INFO);
510: return;
511: }
512: // *
513: // * Quick return if possible
514: // *
515: if (N == 0) return;
516: // *
517: // * The following values are for bookkeeping purposes only. They are
518: // * integer pointers which indicate the portion of the workspace
519: // * used by a particular array in DLAED8 and DLAED9.
520: // *
521: if (ICOMPQ == 1)
522: {
523: LDQ2 = QSIZ;
524: }
525: else
526: {
527: LDQ2 = N;
528: }
529: // *
530: IZ = 1;
531: IDLMDA = IZ + N;
532: IW = IDLMDA + N;
533: IQ2 = IW + N;
534: IS = IQ2 + N * LDQ2;
535: // *
536: INDX = 1;
537: INDXC = INDX + N;
538: COLTYP = INDXC + N;
539: INDXP = COLTYP + N;
540: // *
541: // * Form the z-vector which consists of the last row of Q_1 and the
542: // * first row of Q_2.
543: // *
544: PTR = 1 + (int)Math.Pow(2,TLVLS);
545: for (I = 1; I <= CURLVL - 1; I++)
546: {
547: PTR = PTR + (int)Math.Pow(2, TLVLS - I);
548: }
549: CURR = PTR + CURPBM;
550: this._dlaeda.Run(N, TLVLS, CURLVL, CURPBM, PRMPTR, offset_prmptr, PERM, offset_perm
551: , GIVPTR, offset_givptr, GIVCOL, offset_givcol, GIVNUM, offset_givnum, QSTORE, offset_qstore, QPTR, offset_qptr, ref WORK, IZ + o_work
552: , ref WORK, IZ + N + o_work, ref INFO);
553: // *
554: // * When solving the final problem, we no longer need the stored data,
555: // * so we will overwrite the data from this level onto the previously
556: // * used storage space.
557: // *
558: if (CURLVL == TLVLS)
559: {
560: QPTR[CURR + o_qptr] = 1;
561: PRMPTR[CURR + o_prmptr] = 1;
562: GIVPTR[CURR + o_givptr] = 1;
563: }
564: // *
565: // * Sort and Deflate eigenvalues.
566: // *
567: this._dlaed8.Run(ICOMPQ, ref K, N, QSIZ, ref D, offset_d, ref Q, offset_q
568: , LDQ, ref INDXQ, offset_indxq, ref RHO, CUTPNT, ref WORK, IZ + o_work, ref WORK, IDLMDA + o_work
569: , ref WORK, IQ2 + o_work, LDQ2, ref WORK, IW + o_work, ref PERM, PRMPTR[CURR + o_prmptr] + o_perm, ref GIVPTR[CURR + 1 + o_givptr], ref GIVCOL, 1+(GIVPTR[CURR + o_givptr]) * 2 + o_givcol
570: , ref GIVNUM, 1+(GIVPTR[CURR + o_givptr]) * 2 + o_givnum, ref IWORK, INDXP + o_iwork, ref IWORK, INDX + o_iwork, ref INFO);
571: PRMPTR[CURR + 1 + o_prmptr] = PRMPTR[CURR + o_prmptr] + N;
572: GIVPTR[CURR + 1 + o_givptr] = GIVPTR[CURR + 1 + o_givptr] + GIVPTR[CURR + o_givptr];
573: // *
574: // * Solve Secular Equation.
575: // *
576: if (K != 0)
577: {
578: this._dlaed9.Run(K, 1, K, N, ref D, offset_d, ref WORK, IS + o_work
579: , K, RHO, ref WORK, IDLMDA + o_work, ref WORK, IW + o_work, ref QSTORE, QPTR[CURR + o_qptr] + o_qstore, K
580: , ref INFO);
581: if (INFO != 0) goto LABEL30;
582: if (ICOMPQ == 1)
583: {
584: this._dgemm.Run("N", "N", QSIZ, K, K, ONE
585: , WORK, IQ2 + o_work, LDQ2, QSTORE, QPTR[CURR + o_qptr] + o_qstore, K, ZERO, ref Q, offset_q
586: , LDQ);
587: }
588: QPTR[CURR + 1 + o_qptr] = QPTR[CURR + o_qptr] + (int)Math.Pow(K, 2);
589: // *
590: // * Prepare the INDXQ sorting permutation.
591: // *
592: N1 = K;
593: N2 = N - K;
594: this._dlamrg.Run(N1, N2, D, offset_d, 1, - 1, ref INDXQ, offset_indxq);
595: }
596: else
597: {
598: QPTR[CURR + 1 + o_qptr] = QPTR[CURR + o_qptr];
599: for (I = 1; I <= N; I++)
600: {
601: INDXQ[I + o_indxq] = I;
602: }
603: }
604: // *
605: LABEL30:;
606: return;
607: // *
608: // * End of DLAED7
609: // *
610:
611: #endregion
612:
613: }
614: }
615: }