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: /// DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
27: /// an orthogonal similarity transformation: Q' * A * Q = H .
28: ///
29: ///</summary>
30: public class DGEHD2
31: {
32:
33:
34: #region Dependencies
35:
36: DLARF _dlarf; DLARFG _dlarfg; XERBLA _xerbla;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: const double ONE = 1.0E+0; int I = 0; double AII = 0;
44:
45: #endregion
46:
47: public DGEHD2(DLARF dlarf, DLARFG dlarfg, XERBLA xerbla)
48: {
49:
50:
51: #region Set Dependencies
52:
53: this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla;
54:
55: #endregion
56:
57: }
58:
59: public DGEHD2()
60: {
61:
62:
63: #region Dependencies (Initialization)
64:
65: LSAME lsame = new LSAME();
66: XERBLA xerbla = new XERBLA();
67: DLAMC3 dlamc3 = new DLAMC3();
68: DLAPY2 dlapy2 = new DLAPY2();
69: DNRM2 dnrm2 = new DNRM2();
70: DSCAL dscal = new DSCAL();
71: DGEMV dgemv = new DGEMV(lsame, xerbla);
72: DGER dger = new DGER(xerbla);
73: DLARF dlarf = new DLARF(dgemv, dger, lsame);
74: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
75: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
76: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
77: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
78: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
79: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
80:
81: #endregion
82:
83:
84: #region Set Dependencies
85:
86: this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla;
87:
88: #endregion
89:
90: }
91: /// <summary>
92: /// Purpose
93: /// =======
94: ///
95: /// DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
96: /// an orthogonal similarity transformation: Q' * A * Q = H .
97: ///
98: ///</summary>
99: /// <param name="N">
100: /// (input) INTEGER
101: /// The order of the matrix A. N .GE. 0.
102: ///</param>
103: /// <param name="ILO">
104: /// (input) INTEGER
105: ///</param>
106: /// <param name="IHI">
107: /// (input) INTEGER
108: /// It is assumed that A is already upper triangular in rows
109: /// and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
110: /// set by a previous call to DGEBAL; otherwise they should be
111: /// set to 1 and N respectively. See Further Details.
112: /// 1 .LE. ILO .LE. IHI .LE. max(1,N).
113: ///</param>
114: /// <param name="A">
115: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
116: /// On entry, the n by n general matrix to be reduced.
117: /// On exit, the upper triangle and the first subdiagonal of A
118: /// are overwritten with the upper Hessenberg matrix H, and the
119: /// elements below the first subdiagonal, with the array TAU,
120: /// represent the orthogonal matrix Q as a product of elementary
121: /// reflectors. See Further Details.
122: ///</param>
123: /// <param name="LDA">
124: /// (input) INTEGER
125: /// The leading dimension of the array A. LDA .GE. max(1,N).
126: ///</param>
127: /// <param name="TAU">
128: /// (output) DOUBLE PRECISION array, dimension (N-1)
129: /// The scalar factors of the elementary reflectors (see Further
130: /// Details).
131: ///</param>
132: /// <param name="WORK">
133: /// (workspace) DOUBLE PRECISION array, dimension (N)
134: ///</param>
135: /// <param name="INFO">
136: /// (output) INTEGER
137: /// = 0: successful exit.
138: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
139: ///</param>
140: public void Run(int N, int ILO, int IHI, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
141: , ref double[] WORK, int offset_work, ref int INFO)
142: {
143:
144: #region Array Index Correction
145:
146: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
147:
148: #endregion
149:
150:
151: #region Prolog
152:
153: // *
154: // * -- LAPACK routine (version 3.1) --
155: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
156: // * November 2006
157: // *
158: // * .. Scalar Arguments ..
159: // * ..
160: // * .. Array Arguments ..
161: // * ..
162: // *
163: // * Purpose
164: // * =======
165: // *
166: // * DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
167: // * an orthogonal similarity transformation: Q' * A * Q = H .
168: // *
169: // * Arguments
170: // * =========
171: // *
172: // * N (input) INTEGER
173: // * The order of the matrix A. N >= 0.
174: // *
175: // * ILO (input) INTEGER
176: // * IHI (input) INTEGER
177: // * It is assumed that A is already upper triangular in rows
178: // * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
179: // * set by a previous call to DGEBAL; otherwise they should be
180: // * set to 1 and N respectively. See Further Details.
181: // * 1 <= ILO <= IHI <= max(1,N).
182: // *
183: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
184: // * On entry, the n by n general matrix to be reduced.
185: // * On exit, the upper triangle and the first subdiagonal of A
186: // * are overwritten with the upper Hessenberg matrix H, and the
187: // * elements below the first subdiagonal, with the array TAU,
188: // * represent the orthogonal matrix Q as a product of elementary
189: // * reflectors. See Further Details.
190: // *
191: // * LDA (input) INTEGER
192: // * The leading dimension of the array A. LDA >= max(1,N).
193: // *
194: // * TAU (output) DOUBLE PRECISION array, dimension (N-1)
195: // * The scalar factors of the elementary reflectors (see Further
196: // * Details).
197: // *
198: // * WORK (workspace) DOUBLE PRECISION array, dimension (N)
199: // *
200: // * INFO (output) INTEGER
201: // * = 0: successful exit.
202: // * < 0: if INFO = -i, the i-th argument had an illegal value.
203: // *
204: // * Further Details
205: // * ===============
206: // *
207: // * The matrix Q is represented as a product of (ihi-ilo) elementary
208: // * reflectors
209: // *
210: // * Q = H(ilo) H(ilo+1) . . . H(ihi-1).
211: // *
212: // * Each H(i) has the form
213: // *
214: // * H(i) = I - tau * v * v'
215: // *
216: // * where tau is a real scalar, and v is a real vector with
217: // * v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
218: // * exit in A(i+2:ihi,i), and tau in TAU(i).
219: // *
220: // * The contents of A are illustrated by the following example, with
221: // * n = 7, ilo = 2 and ihi = 6:
222: // *
223: // * on entry, on exit,
224: // *
225: // * ( a a a a a a a ) ( a a h h h h a )
226: // * ( a a a a a a ) ( a h h h h a )
227: // * ( a a a a a a ) ( h h h h h h )
228: // * ( a a a a a a ) ( v2 h h h h h )
229: // * ( a a a a a a ) ( v2 v3 h h h h )
230: // * ( a a a a a a ) ( v2 v3 v4 h h h )
231: // * ( a ) ( a )
232: // *
233: // * where a denotes an element of the original matrix A, h denotes a
234: // * modified element of the upper Hessenberg matrix H, and vi denotes an
235: // * element of the vector defining H(i).
236: // *
237: // * =====================================================================
238: // *
239: // * .. Parameters ..
240: // * ..
241: // * .. Local Scalars ..
242: // * ..
243: // * .. External Subroutines ..
244: // * ..
245: // * .. Intrinsic Functions ..
246: // INTRINSIC MAX, MIN;
247: // * ..
248: // * .. Executable Statements ..
249: // *
250: // * Test the input parameters
251: // *
252:
253: #endregion
254:
255:
256: #region Body
257:
258: INFO = 0;
259: if (N < 0)
260: {
261: INFO = - 1;
262: }
263: else
264: {
265: if (ILO < 1 || ILO > Math.Max(1, N))
266: {
267: INFO = - 2;
268: }
269: else
270: {
271: if (IHI < Math.Min(ILO, N) || IHI > N)
272: {
273: INFO = - 3;
274: }
275: else
276: {
277: if (LDA < Math.Max(1, N))
278: {
279: INFO = - 5;
280: }
281: }
282: }
283: }
284: if (INFO != 0)
285: {
286: this._xerbla.Run("DGEHD2", - INFO);
287: return;
288: }
289: // *
290: for (I = ILO; I <= IHI - 1; I++)
291: {
292: // *
293: // * Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
294: // *
295: this._dlarfg.Run(IHI - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, N)+I * LDA + o_a, 1, ref TAU[I + o_tau]);
296: AII = A[I + 1+I * LDA + o_a];
297: A[I + 1+I * LDA + o_a] = ONE;
298: // *
299: // * Apply H(i) to A(1:ihi,i+1:ihi) from the right
300: // *
301: this._dlarf.Run("Right", IHI, IHI - I, A, I + 1+I * LDA + o_a, 1, TAU[I + o_tau]
302: , ref A, 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
303: // *
304: // * Apply H(i) to A(i+1:ihi,i+1:n) from the left
305: // *
306: this._dlarf.Run("Left", IHI - I, N - I, A, I + 1+I * LDA + o_a, 1, TAU[I + o_tau]
307: , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
308: // *
309: A[I + 1+I * LDA + o_a] = AII;
310: }
311: // *
312: return;
313: // *
314: // * End of DGEHD2
315: // *
316:
317: #endregion
318:
319: }
320: }
321: }