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: /// DLAED1 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 eigenvectors of a tridiagonal matrix. DLAED7 handles
30: /// the case in which eigenvalues only or eigenvalues and eigenvectors
31: /// of a full symmetric matrix (which was reduced to tridiagonal form)
32: /// 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 DLAED2.
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 DLAED3).
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 DLAED1
61: {
62:
63:
64: #region Dependencies
65:
66: DCOPY _dcopy; DLAED2 _dlaed2; DLAED3 _dlaed3; DLAMRG _dlamrg; XERBLA _xerbla;
67:
68: #endregion
69:
70:
71: #region Fields
72:
73: int COLTYP = 0; int I = 0; int IDLMDA = 0; int INDX = 0; int INDXC = 0; int INDXP = 0; int IQ2 = 0; int IS = 0;
74: int IW = 0;int IZ = 0; int K = 0; int N1 = 0; int N2 = 0; int ZPP1 = 0;
75:
76: #endregion
77:
78: public DLAED1(DCOPY dcopy, DLAED2 dlaed2, DLAED3 dlaed3, DLAMRG dlamrg, XERBLA xerbla)
79: {
80:
81:
82: #region Set Dependencies
83:
84: this._dcopy = dcopy; this._dlaed2 = dlaed2; this._dlaed3 = dlaed3; this._dlamrg = dlamrg; this._xerbla = xerbla;
85:
86: #endregion
87:
88: }
89:
90: public DLAED1()
91: {
92:
93:
94: #region Dependencies (Initialization)
95:
96: DCOPY dcopy = new DCOPY();
97: IDAMAX idamax = new IDAMAX();
98: LSAME lsame = new LSAME();
99: DLAMC3 dlamc3 = new DLAMC3();
100: DLAPY2 dlapy2 = new DLAPY2();
101: DLAMRG dlamrg = new DLAMRG();
102: DROT drot = new DROT();
103: DSCAL dscal = new DSCAL();
104: XERBLA xerbla = new XERBLA();
105: DNRM2 dnrm2 = new DNRM2();
106: DLAED5 dlaed5 = new DLAED5();
107: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
108: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
109: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
110: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
111: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
112: DLACPY dlacpy = new DLACPY(lsame);
113: DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
114: DGEMM dgemm = new DGEMM(lsame, xerbla);
115: DLAED6 dlaed6 = new DLAED6(dlamch);
116: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
117: DLASET dlaset = new DLASET(lsame);
118: DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
119:
120: #endregion
121:
122:
123: #region Set Dependencies
124:
125: this._dcopy = dcopy; this._dlaed2 = dlaed2; this._dlaed3 = dlaed3; this._dlamrg = dlamrg; this._xerbla = xerbla;
126:
127: #endregion
128:
129: }
130: /// <summary>
131: /// Purpose
132: /// =======
133: ///
134: /// DLAED1 computes the updated eigensystem of a diagonal
135: /// matrix after modification by a rank-one symmetric matrix. This
136: /// routine is used only for the eigenproblem which requires all
137: /// eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles
138: /// the case in which eigenvalues only or eigenvalues and eigenvectors
139: /// of a full symmetric matrix (which was reduced to tridiagonal form)
140: /// are desired.
141: ///
142: /// T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
143: ///
144: /// where Z = Q'u, u is a vector of length N with ones in the
145: /// CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
146: ///
147: /// The eigenvectors of the original matrix are stored in Q, and the
148: /// eigenvalues are in D. The algorithm consists of three stages:
149: ///
150: /// The first stage consists of deflating the size of the problem
151: /// when there are multiple eigenvalues or if there is a zero in
152: /// the Z vector. For each such occurence the dimension of the
153: /// secular equation problem is reduced by one. This stage is
154: /// performed by the routine DLAED2.
155: ///
156: /// The second stage consists of calculating the updated
157: /// eigenvalues. This is done by finding the roots of the secular
158: /// equation via the routine DLAED4 (as called by DLAED3).
159: /// This routine also calculates the eigenvectors of the current
160: /// problem.
161: ///
162: /// The final stage consists of computing the updated eigenvectors
163: /// directly using the updated eigenvalues. The eigenvectors for
164: /// the current problem are multiplied with the eigenvectors from
165: /// the overall problem.
166: ///
167: ///</summary>
168: /// <param name="N">
169: /// (input) INTEGER
170: /// The dimension of the symmetric tridiagonal matrix. N .GE. 0.
171: ///</param>
172: /// <param name="D">
173: /// (input/output) DOUBLE PRECISION array, dimension (N)
174: /// On entry, the eigenvalues of the rank-1-perturbed matrix.
175: /// On exit, the eigenvalues of the repaired matrix.
176: ///</param>
177: /// <param name="Q">
178: /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
179: /// On entry, the eigenvectors of the rank-1-perturbed matrix.
180: /// On exit, the eigenvectors of the repaired tridiagonal matrix.
181: ///</param>
182: /// <param name="LDQ">
183: /// (input) INTEGER
184: /// The leading dimension of the array Q. LDQ .GE. max(1,N).
185: ///</param>
186: /// <param name="INDXQ">
187: /// (input/output) INTEGER array, dimension (N)
188: /// On entry, the permutation which separately sorts the two
189: /// subproblems in D into ascending order.
190: /// On exit, the permutation which will reintegrate the
191: /// subproblems back into sorted order,
192: /// i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
193: ///</param>
194: /// <param name="RHO">
195: /// (input) DOUBLE PRECISION
196: /// The subdiagonal entry used to create the rank-1 modification.
197: ///</param>
198: /// <param name="CUTPNT">
199: /// (input) INTEGER
200: /// The location of the last eigenvalue in the leading sub-matrix.
201: /// min(1,N) .LE. CUTPNT .LE. N/2.
202: ///</param>
203: /// <param name="WORK">
204: /// (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)
205: ///</param>
206: /// <param name="IWORK">
207: /// (workspace) INTEGER array, dimension (4*N)
208: ///</param>
209: /// <param name="INFO">
210: /// (output) INTEGER
211: /// = 0: successful exit.
212: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
213: /// .GT. 0: if INFO = 1, an eigenvalue did not converge
214: ///</param>
215: public void Run(int N, ref double[] D, int offset_d, ref double[] Q, int offset_q, int LDQ, ref int[] INDXQ, int offset_indxq, ref double RHO
216: , int CUTPNT, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
217: {
218:
219: #region Array Index Correction
220:
221: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_indxq = -1 + offset_indxq;
222: int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
223:
224: #endregion
225:
226:
227: #region Prolog
228:
229: // *
230: // * -- LAPACK routine (version 3.1) --
231: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
232: // * November 2006
233: // *
234: // * .. Scalar Arguments ..
235: // * ..
236: // * .. Array Arguments ..
237: // * ..
238: // *
239: // * Purpose
240: // * =======
241: // *
242: // * DLAED1 computes the updated eigensystem of a diagonal
243: // * matrix after modification by a rank-one symmetric matrix. This
244: // * routine is used only for the eigenproblem which requires all
245: // * eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles
246: // * the case in which eigenvalues only or eigenvalues and eigenvectors
247: // * of a full symmetric matrix (which was reduced to tridiagonal form)
248: // * are desired.
249: // *
250: // * T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)
251: // *
252: // * where Z = Q'u, u is a vector of length N with ones in the
253: // * CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.
254: // *
255: // * The eigenvectors of the original matrix are stored in Q, and the
256: // * eigenvalues are in D. The algorithm consists of three stages:
257: // *
258: // * The first stage consists of deflating the size of the problem
259: // * when there are multiple eigenvalues or if there is a zero in
260: // * the Z vector. For each such occurence the dimension of the
261: // * secular equation problem is reduced by one. This stage is
262: // * performed by the routine DLAED2.
263: // *
264: // * The second stage consists of calculating the updated
265: // * eigenvalues. This is done by finding the roots of the secular
266: // * equation via the routine DLAED4 (as called by DLAED3).
267: // * This routine also calculates the eigenvectors of the current
268: // * problem.
269: // *
270: // * The final stage consists of computing the updated eigenvectors
271: // * directly using the updated eigenvalues. The eigenvectors for
272: // * the current problem are multiplied with the eigenvectors from
273: // * the overall problem.
274: // *
275: // * Arguments
276: // * =========
277: // *
278: // * N (input) INTEGER
279: // * The dimension of the symmetric tridiagonal matrix. N >= 0.
280: // *
281: // * D (input/output) DOUBLE PRECISION array, dimension (N)
282: // * On entry, the eigenvalues of the rank-1-perturbed matrix.
283: // * On exit, the eigenvalues of the repaired matrix.
284: // *
285: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
286: // * On entry, the eigenvectors of the rank-1-perturbed matrix.
287: // * On exit, the eigenvectors of the repaired tridiagonal matrix.
288: // *
289: // * LDQ (input) INTEGER
290: // * The leading dimension of the array Q. LDQ >= max(1,N).
291: // *
292: // * INDXQ (input/output) INTEGER array, dimension (N)
293: // * On entry, the permutation which separately sorts the two
294: // * subproblems in D into ascending order.
295: // * On exit, the permutation which will reintegrate the
296: // * subproblems back into sorted order,
297: // * i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.
298: // *
299: // * RHO (input) DOUBLE PRECISION
300: // * The subdiagonal entry used to create the rank-1 modification.
301: // *
302: // * CUTPNT (input) INTEGER
303: // * The location of the last eigenvalue in the leading sub-matrix.
304: // * min(1,N) <= CUTPNT <= N/2.
305: // *
306: // * WORK (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)
307: // *
308: // * IWORK (workspace) INTEGER array, dimension (4*N)
309: // *
310: // * INFO (output) INTEGER
311: // * = 0: successful exit.
312: // * < 0: if INFO = -i, the i-th argument had an illegal value.
313: // * > 0: if INFO = 1, an eigenvalue did not converge
314: // *
315: // * Further Details
316: // * ===============
317: // *
318: // * Based on contributions by
319: // * Jeff Rutter, Computer Science Division, University of California
320: // * at Berkeley, USA
321: // * Modified by Francoise Tisseur, University of Tennessee.
322: // *
323: // * =====================================================================
324: // *
325: // * .. Local Scalars ..
326: // * ..
327: // * .. External Subroutines ..
328: // * ..
329: // * .. Intrinsic Functions ..
330: // INTRINSIC MAX, MIN;
331: // * ..
332: // * .. Executable Statements ..
333: // *
334: // * Test the input parameters.
335: // *
336:
337: #endregion
338:
339:
340: #region Body
341:
342: INFO = 0;
343: // *
344: if (N < 0)
345: {
346: INFO = - 1;
347: }
348: else
349: {
350: if (LDQ < Math.Max(1, N))
351: {
352: INFO = - 4;
353: }
354: else
355: {
356: if (Math.Min(1, N / 2) > CUTPNT || (N / 2) < CUTPNT)
357: {
358: INFO = - 7;
359: }
360: }
361: }
362: if (INFO != 0)
363: {
364: this._xerbla.Run("DLAED1", - INFO);
365: return;
366: }
367: // *
368: // * Quick return if possible
369: // *
370: if (N == 0) return;
371: // *
372: // * The following values are integer pointers which indicate
373: // * the portion of the workspace
374: // * used by a particular array in DLAED2 and DLAED3.
375: // *
376: IZ = 1;
377: IDLMDA = IZ + N;
378: IW = IDLMDA + N;
379: IQ2 = IW + N;
380: // *
381: INDX = 1;
382: INDXC = INDX + N;
383: COLTYP = INDXC + N;
384: INDXP = COLTYP + N;
385: // *
386: // *
387: // * Form the z-vector which consists of the last row of Q_1 and the
388: // * first row of Q_2.
389: // *
390: this._dcopy.Run(CUTPNT, Q, CUTPNT+1 * LDQ + o_q, LDQ, ref WORK, IZ + o_work, 1);
391: ZPP1 = CUTPNT + 1;
392: this._dcopy.Run(N - CUTPNT, Q, ZPP1+ZPP1 * LDQ + o_q, LDQ, ref WORK, IZ + CUTPNT + o_work, 1);
393: // *
394: // * Deflate eigenvalues.
395: // *
396: this._dlaed2.Run(ref K, N, CUTPNT, ref D, offset_d, ref Q, offset_q, LDQ
397: , ref INDXQ, offset_indxq, ref RHO, ref WORK, IZ + o_work, ref WORK, IDLMDA + o_work, ref WORK, IW + o_work, ref WORK, IQ2 + o_work
398: , ref IWORK, INDX + o_iwork, ref IWORK, INDXC + o_iwork, ref IWORK, INDXP + o_iwork, ref IWORK, COLTYP + o_iwork, ref INFO);
399: // *
400: if (INFO != 0) goto LABEL20;
401: // *
402: // * Solve Secular Equation.
403: // *
404: if (K != 0)
405: {
406: IS = (IWORK[COLTYP + o_iwork] + IWORK[COLTYP + 1 + o_iwork]) * CUTPNT + (IWORK[COLTYP + 1 + o_iwork] + IWORK[COLTYP + 2 + o_iwork]) * (N - CUTPNT) + IQ2;
407: this._dlaed3.Run(K, N, CUTPNT, ref D, offset_d, ref Q, offset_q, LDQ
408: , RHO, ref WORK, IDLMDA + o_work, WORK, IQ2 + o_work, IWORK, INDXC + o_iwork, IWORK, COLTYP + o_iwork, ref WORK, IW + o_work
409: , ref WORK, IS + o_work, ref INFO);
410: if (INFO != 0) goto LABEL20;
411: // *
412: // * Prepare the INDXQ sorting permutation.
413: // *
414: N1 = K;
415: N2 = N - K;
416: this._dlamrg.Run(N1, N2, D, offset_d, 1, - 1, ref INDXQ, offset_indxq);
417: }
418: else
419: {
420: for (I = 1; I <= N; I++)
421: {
422: INDXQ[I + o_indxq] = I;
423: }
424: }
425: // *
426: LABEL20:;
427: return;
428: // *
429: // * End of DLAED1
430: // *
431:
432: #endregion
433:
434: }
435: }
436: }