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: /// DLAED3 finds the roots of the secular equation, as defined by the
27: /// values in D, W, and RHO, between 1 and K. It makes the
28: /// appropriate calls to DLAED4 and then updates the eigenvectors by
29: /// multiplying the matrix of eigenvectors of the pair of eigensystems
30: /// being combined by the matrix of eigenvectors of the K-by-K system
31: /// which is solved here.
32: ///
33: /// This code makes very mild assumptions about floating point
34: /// arithmetic. It will work on machines with a guard digit in
35: /// add/subtract, or on those binary machines without guard digits
36: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
37: /// It could conceivably fail on hexadecimal or decimal machines
38: /// without guard digits, but we know of none.
39: ///
40: ///</summary>
41: public class DLAED3
42: {
43:
44:
45: #region Dependencies
46:
47: DLAMC3 _dlamc3; DNRM2 _dnrm2; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLAED4 _dlaed4; DLASET _dlaset; XERBLA _xerbla;
48:
49: #endregion
50:
51:
52: #region Fields
53:
54: const double ONE = 1.0E0; const double ZERO = 0.0E0; int I = 0; int II = 0; int IQ2 = 0; int J = 0; int N12 = 0;
55: int N2 = 0;int N23 = 0; double TEMP = 0;
56:
57: #endregion
58:
59: public DLAED3(DLAMC3 dlamc3, DNRM2 dnrm2, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLAED4 dlaed4, DLASET dlaset, XERBLA xerbla)
60: {
61:
62:
63: #region Set Dependencies
64:
65: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
66: this._dlaed4 = dlaed4;this._dlaset = dlaset; this._xerbla = xerbla;
67:
68: #endregion
69:
70: }
71:
72: public DLAED3()
73: {
74:
75:
76: #region Dependencies (Initialization)
77:
78: DLAMC3 dlamc3 = new DLAMC3();
79: DNRM2 dnrm2 = new DNRM2();
80: DCOPY dcopy = new DCOPY();
81: LSAME lsame = new LSAME();
82: XERBLA xerbla = new XERBLA();
83: DLAED5 dlaed5 = new DLAED5();
84: DGEMM dgemm = new DGEMM(lsame, xerbla);
85: DLACPY dlacpy = new DLACPY(lsame);
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91: DLAED6 dlaed6 = new DLAED6(dlamch);
92: DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
93: DLASET dlaset = new DLASET(lsame);
94:
95: #endregion
96:
97:
98: #region Set Dependencies
99:
100: this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy;
101: this._dlaed4 = dlaed4;this._dlaset = dlaset; this._xerbla = xerbla;
102:
103: #endregion
104:
105: }
106: /// <summary>
107: /// Purpose
108: /// =======
109: ///
110: /// DLAED3 finds the roots of the secular equation, as defined by the
111: /// values in D, W, and RHO, between 1 and K. It makes the
112: /// appropriate calls to DLAED4 and then updates the eigenvectors by
113: /// multiplying the matrix of eigenvectors of the pair of eigensystems
114: /// being combined by the matrix of eigenvectors of the K-by-K system
115: /// which is solved here.
116: ///
117: /// This code makes very mild assumptions about floating point
118: /// arithmetic. It will work on machines with a guard digit in
119: /// add/subtract, or on those binary machines without guard digits
120: /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
121: /// It could conceivably fail on hexadecimal or decimal machines
122: /// without guard digits, but we know of none.
123: ///
124: ///</summary>
125: /// <param name="K">
126: /// (input) INTEGER
127: /// The number of terms in the rational function to be solved by
128: /// DLAED4. K .GE. 0.
129: ///</param>
130: /// <param name="N">
131: /// (input) INTEGER
132: /// The number of rows and columns in the Q matrix.
133: /// N .GE. K (deflation may result in N.GT.K).
134: ///</param>
135: /// <param name="N1">
136: /// (input) INTEGER
137: /// The location of the last eigenvalue in the leading submatrix.
138: /// min(1,N) .LE. N1 .LE. N/2.
139: ///</param>
140: /// <param name="D">
141: /// (output) DOUBLE PRECISION array, dimension (N)
142: /// D(I) contains the updated eigenvalues for
143: /// 1 .LE. I .LE. K.
144: ///</param>
145: /// <param name="Q">
146: /// (output) DOUBLE PRECISION array, dimension (LDQ,N)
147: /// Initially the first K columns are used as workspace.
148: /// On output the columns 1 to K contain
149: /// the updated eigenvectors.
150: ///</param>
151: /// <param name="LDQ">
152: /// (input) INTEGER
153: /// The leading dimension of the array Q. LDQ .GE. max(1,N).
154: ///</param>
155: /// <param name="RHO">
156: /// (input) DOUBLE PRECISION
157: /// The value of the parameter in the rank one update equation.
158: /// RHO .GE. 0 required.
159: ///</param>
160: /// <param name="DLAMDA">
161: /// (input/output) DOUBLE PRECISION array, dimension (K)
162: /// The first K elements of this array contain the old roots
163: /// of the deflated updating problem. These are the poles
164: /// of the secular equation. May be changed on output by
165: /// having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
166: /// Cray-2, or Cray C-90, as described above.
167: ///</param>
168: /// <param name="Q2">
169: /// (input) DOUBLE PRECISION array, dimension (LDQ2, N)
170: /// The first K columns of this matrix contain the non-deflated
171: /// eigenvectors for the split problem.
172: ///</param>
173: /// <param name="INDX">
174: /// (input) INTEGER array, dimension (N)
175: /// The permutation used to arrange the columns of the deflated
176: /// Q matrix into three groups (see DLAED2).
177: /// The rows of the eigenvectors found by DLAED4 must be likewise
178: /// permuted before the matrix multiply can take place.
179: ///</param>
180: /// <param name="CTOT">
181: /// (input) INTEGER array, dimension (4)
182: /// A count of the total number of the various types of columns
183: /// in Q, as described in INDX. The fourth column type is any
184: /// column which has been deflated.
185: ///</param>
186: /// <param name="W">
187: /// (input/output) DOUBLE PRECISION array, dimension (K)
188: /// The first K elements of this array contain the components
189: /// of the deflation-adjusted updating vector. Destroyed on
190: /// output.
191: ///</param>
192: /// <param name="S">
193: /// (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K
194: /// Will contain the eigenvectors of the repaired matrix which
195: /// will be multiplied by the previously accumulated eigenvectors
196: /// to update the system.
197: ///</param>
198: /// <param name="INFO">
199: /// (output) INTEGER
200: /// = 0: successful exit.
201: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
202: /// .GT. 0: if INFO = 1, an eigenvalue did not converge
203: ///</param>
204: public void Run(int K, int N, int N1, ref double[] D, int offset_d, ref double[] Q, int offset_q, int LDQ
205: , double RHO, ref double[] DLAMDA, int offset_dlamda, double[] Q2, int offset_q2, int[] INDX, int offset_indx, int[] CTOT, int offset_ctot, ref double[] W, int offset_w
206: , ref double[] S, int offset_s, ref int INFO)
207: {
208:
209: #region Array Index Correction
210:
211: int o_d = -1 + offset_d; int o_q = -1 - LDQ + offset_q; int o_dlamda = -1 + offset_dlamda;
212: int o_q2 = -1 + offset_q2; int o_indx = -1 + offset_indx; int o_ctot = -1 + offset_ctot; int o_w = -1 + offset_w;
213: int o_s = -1 + offset_s;
214:
215: #endregion
216:
217:
218: #region Prolog
219:
220: // *
221: // * -- LAPACK routine (version 3.1) --
222: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
223: // * November 2006
224: // *
225: // * .. Scalar Arguments ..
226: // * ..
227: // * .. Array Arguments ..
228: // * ..
229: // *
230: // * Purpose
231: // * =======
232: // *
233: // * DLAED3 finds the roots of the secular equation, as defined by the
234: // * values in D, W, and RHO, between 1 and K. It makes the
235: // * appropriate calls to DLAED4 and then updates the eigenvectors by
236: // * multiplying the matrix of eigenvectors of the pair of eigensystems
237: // * being combined by the matrix of eigenvectors of the K-by-K system
238: // * which is solved here.
239: // *
240: // * This code makes very mild assumptions about floating point
241: // * arithmetic. It will work on machines with a guard digit in
242: // * add/subtract, or on those binary machines without guard digits
243: // * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
244: // * It could conceivably fail on hexadecimal or decimal machines
245: // * without guard digits, but we know of none.
246: // *
247: // * Arguments
248: // * =========
249: // *
250: // * K (input) INTEGER
251: // * The number of terms in the rational function to be solved by
252: // * DLAED4. K >= 0.
253: // *
254: // * N (input) INTEGER
255: // * The number of rows and columns in the Q matrix.
256: // * N >= K (deflation may result in N>K).
257: // *
258: // * N1 (input) INTEGER
259: // * The location of the last eigenvalue in the leading submatrix.
260: // * min(1,N) <= N1 <= N/2.
261: // *
262: // * D (output) DOUBLE PRECISION array, dimension (N)
263: // * D(I) contains the updated eigenvalues for
264: // * 1 <= I <= K.
265: // *
266: // * Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
267: // * Initially the first K columns are used as workspace.
268: // * On output the columns 1 to K contain
269: // * the updated eigenvectors.
270: // *
271: // * LDQ (input) INTEGER
272: // * The leading dimension of the array Q. LDQ >= max(1,N).
273: // *
274: // * RHO (input) DOUBLE PRECISION
275: // * The value of the parameter in the rank one update equation.
276: // * RHO >= 0 required.
277: // *
278: // * DLAMDA (input/output) DOUBLE PRECISION array, dimension (K)
279: // * The first K elements of this array contain the old roots
280: // * of the deflated updating problem. These are the poles
281: // * of the secular equation. May be changed on output by
282: // * having lowest order bit set to zero on Cray X-MP, Cray Y-MP,
283: // * Cray-2, or Cray C-90, as described above.
284: // *
285: // * Q2 (input) DOUBLE PRECISION array, dimension (LDQ2, N)
286: // * The first K columns of this matrix contain the non-deflated
287: // * eigenvectors for the split problem.
288: // *
289: // * INDX (input) INTEGER array, dimension (N)
290: // * The permutation used to arrange the columns of the deflated
291: // * Q matrix into three groups (see DLAED2).
292: // * The rows of the eigenvectors found by DLAED4 must be likewise
293: // * permuted before the matrix multiply can take place.
294: // *
295: // * CTOT (input) INTEGER array, dimension (4)
296: // * A count of the total number of the various types of columns
297: // * in Q, as described in INDX. The fourth column type is any
298: // * column which has been deflated.
299: // *
300: // * W (input/output) DOUBLE PRECISION array, dimension (K)
301: // * The first K elements of this array contain the components
302: // * of the deflation-adjusted updating vector. Destroyed on
303: // * output.
304: // *
305: // * S (workspace) DOUBLE PRECISION array, dimension (N1 + 1)*K
306: // * Will contain the eigenvectors of the repaired matrix which
307: // * will be multiplied by the previously accumulated eigenvectors
308: // * to update the system.
309: // *
310: // * LDS (input) INTEGER
311: // * The leading dimension of S. LDS >= max(1,K).
312: // *
313: // * INFO (output) INTEGER
314: // * = 0: successful exit.
315: // * < 0: if INFO = -i, the i-th argument had an illegal value.
316: // * > 0: if INFO = 1, an eigenvalue did not converge
317: // *
318: // * Further Details
319: // * ===============
320: // *
321: // * Based on contributions by
322: // * Jeff Rutter, Computer Science Division, University of California
323: // * at Berkeley, USA
324: // * Modified by Francoise Tisseur, University of Tennessee.
325: // *
326: // * =====================================================================
327: // *
328: // * .. Parameters ..
329: // * ..
330: // * .. Local Scalars ..
331: // * ..
332: // * .. External Functions ..
333: // * ..
334: // * .. External Subroutines ..
335: // * ..
336: // * .. Intrinsic Functions ..
337: // INTRINSIC MAX, SIGN, SQRT;
338: // * ..
339: // * .. Executable Statements ..
340: // *
341: // * Test the input parameters.
342: // *
343:
344: #endregion
345:
346:
347: #region Body
348:
349: INFO = 0;
350: // *
351: if (K < 0)
352: {
353: INFO = - 1;
354: }
355: else
356: {
357: if (N < K)
358: {
359: INFO = - 2;
360: }
361: else
362: {
363: if (LDQ < Math.Max(1, N))
364: {
365: INFO = - 6;
366: }
367: }
368: }
369: if (INFO != 0)
370: {
371: this._xerbla.Run("DLAED3", - INFO);
372: return;
373: }
374: // *
375: // * Quick return if possible
376: // *
377: if (K == 0) return;
378: // *
379: // * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
380: // * be computed with high relative accuracy (barring over/underflow).
381: // * This is a problem on machines without a guard digit in
382: // * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
383: // * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
384: // * which on any of these machines zeros out the bottommost
385: // * bit of DLAMDA(I) if it is 1; this makes the subsequent
386: // * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
387: // * occurs. On binary machines with a guard digit (almost all
388: // * machines) it does not change DLAMDA(I) at all. On hexadecimal
389: // * and decimal machines with a guard digit, it slightly
390: // * changes the bottommost bits of DLAMDA(I). It does not account
391: // * for hexadecimal or decimal machines without guard digits
392: // * (we know of none). We use a subroutine call to compute
393: // * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
394: // * this code.
395: // *
396: for (I = 1; I <= K; I++)
397: {
398: DLAMDA[I + o_dlamda] = this._dlamc3.Run(DLAMDA[I + o_dlamda], DLAMDA[I + o_dlamda]) - DLAMDA[I + o_dlamda];
399: }
400: // *
401: for (J = 1; J <= K; J++)
402: {
403: this._dlaed4.Run(K, J, DLAMDA, offset_dlamda, W, offset_w, ref Q, 1+J * LDQ + o_q, RHO
404: , ref D[J + o_d], ref INFO);
405: // *
406: // * If the zero finder fails, the computation is terminated.
407: // *
408: if (INFO != 0) goto LABEL120;
409: }
410: // *
411: if (K == 1) goto LABEL110;
412: if (K == 2)
413: {
414: for (J = 1; J <= K; J++)
415: {
416: W[1 + o_w] = Q[1+J * LDQ + o_q];
417: W[2 + o_w] = Q[2+J * LDQ + o_q];
418: II = INDX[1 + o_indx];
419: Q[1+J * LDQ + o_q] = W[II + o_w];
420: II = INDX[2 + o_indx];
421: Q[2+J * LDQ + o_q] = W[II + o_w];
422: }
423: goto LABEL110;
424: }
425: // *
426: // * Compute updated W.
427: // *
428: this._dcopy.Run(K, W, offset_w, 1, ref S, offset_s, 1);
429: // *
430: // * Initialize W(I) = Q(I,I)
431: // *
432: this._dcopy.Run(K, Q, offset_q, LDQ + 1, ref W, offset_w, 1);
433: for (J = 1; J <= K; J++)
434: {
435: for (I = 1; I <= J - 1; I++)
436: {
437: W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
438: }
439: for (I = J + 1; I <= K; I++)
440: {
441: W[I + o_w] = W[I + o_w] * (Q[I+J * LDQ + o_q] / (DLAMDA[I + o_dlamda] - DLAMDA[J + o_dlamda]));
442: }
443: }
444: for (I = 1; I <= K; I++)
445: {
446: W[I + o_w] = FortranLib.Sign(Math.Sqrt( - W[I + o_w]),S[I + o_s]);
447: }
448: // *
449: // * Compute eigenvectors of the modified rank-1 modification.
450: // *
451: for (J = 1; J <= K; J++)
452: {
453: for (I = 1; I <= K; I++)
454: {
455: S[I + o_s] = W[I + o_w] / Q[I+J * LDQ + o_q];
456: }
457: TEMP = this._dnrm2.Run(K, S, offset_s, 1);
458: for (I = 1; I <= K; I++)
459: {
460: II = INDX[I + o_indx];
461: Q[I+J * LDQ + o_q] = S[II + o_s] / TEMP;
462: }
463: }
464: // *
465: // * Compute the updated eigenvectors.
466: // *
467: LABEL110:;
468: // *
469: N2 = N - N1;
470: N12 = CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
471: N23 = CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
472: // *
473: this._dlacpy.Run("A", N23, K, Q, CTOT[1 + o_ctot] + 1+1 * LDQ + o_q, LDQ, ref S, offset_s
474: , N23);
475: IQ2 = N1 * N12 + 1;
476: if (N23 != 0)
477: {
478: this._dgemm.Run("N", "N", N2, K, N23, ONE
479: , Q2, IQ2 + o_q2, N2, S, offset_s, N23, ZERO, ref Q, N1 + 1+1 * LDQ + o_q
480: , LDQ);
481: }
482: else
483: {
484: this._dlaset.Run("A", N2, K, ZERO, ZERO, ref Q, N1 + 1+1 * LDQ + o_q
485: , LDQ);
486: }
487: // *
488: this._dlacpy.Run("A", N12, K, Q, offset_q, LDQ, ref S, offset_s
489: , N12);
490: if (N12 != 0)
491: {
492: this._dgemm.Run("N", "N", N1, K, N12, ONE
493: , Q2, offset_q2, N1, S, offset_s, N12, ZERO, ref Q, offset_q
494: , LDQ);
495: }
496: else
497: {
498: this._dlaset.Run("A", N1, K, ZERO, ZERO, ref Q, 1+1 * LDQ + o_q
499: , LDQ);
500: }
501: // *
502: // *
503: LABEL120:;
504: return;
505: // *
506: // * End of DLAED3
507: // *
508:
509: #endregion
510:
511: }
512: }
513: }