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: /// DGEHRD 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 DGEHRD
31: {
32:
33:
34: #region Dependencies
35:
36: DAXPY _daxpy; DGEHD2 _dgehd2; DGEMM _dgemm; DLAHR2 _dlahr2; DLARFB _dlarfb; DTRMM _dtrmm; XERBLA _xerbla; ILAENV _ilaenv;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: const int NBMAX = 64; const int LDT = NBMAX + 1; const double ZERO = 0.0E+0; const double ONE = 1.0E+0;
44: bool LQUERY = false;int I = 0; int IB = 0; int IINFO = 0; int IWS = 0; int J = 0; int LDWORK = 0; int LWKOPT = 0;
45: int NB = 0;int NBMIN = 0; int NH = 0; int NX = 0; double EI = 0; double[] T = new double[LDT * NBMAX]; int offset_t = 0;
46:
47: #endregion
48:
49: public DGEHRD(DAXPY daxpy, DGEHD2 dgehd2, DGEMM dgemm, DLAHR2 dlahr2, DLARFB dlarfb, DTRMM dtrmm, XERBLA xerbla, ILAENV ilaenv)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._daxpy = daxpy; this._dgehd2 = dgehd2; this._dgemm = dgemm; this._dlahr2 = dlahr2; this._dlarfb = dlarfb;
56: this._dtrmm = dtrmm;this._xerbla = xerbla; this._ilaenv = ilaenv;
57:
58: #endregion
59:
60: }
61:
62: public DGEHRD()
63: {
64:
65:
66: #region Dependencies (Initialization)
67:
68: DAXPY daxpy = new DAXPY();
69: LSAME lsame = new LSAME();
70: XERBLA xerbla = new XERBLA();
71: DLAMC3 dlamc3 = new DLAMC3();
72: DLAPY2 dlapy2 = new DLAPY2();
73: DNRM2 dnrm2 = new DNRM2();
74: DSCAL dscal = new DSCAL();
75: DCOPY dcopy = new DCOPY();
76: IEEECK ieeeck = new IEEECK();
77: IPARMQ iparmq = new IPARMQ();
78: DGEMV dgemv = new DGEMV(lsame, xerbla);
79: DGER dger = new DGER(xerbla);
80: DLARF dlarf = new DLARF(dgemv, dger, lsame);
81: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
82: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
83: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
84: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
85: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
86: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
87: DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
88: DGEMM dgemm = new DGEMM(lsame, xerbla);
89: DLACPY dlacpy = new DLACPY(lsame);
90: DTRMM dtrmm = new DTRMM(lsame, xerbla);
91: DTRMV dtrmv = new DTRMV(lsame, xerbla);
92: DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
93: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
94: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
95:
96: #endregion
97:
98:
99: #region Set Dependencies
100:
101: this._daxpy = daxpy; this._dgehd2 = dgehd2; this._dgemm = dgemm; this._dlahr2 = dlahr2; this._dlarfb = dlarfb;
102: this._dtrmm = dtrmm;this._xerbla = xerbla; this._ilaenv = ilaenv;
103:
104: #endregion
105:
106: }
107: /// <summary>
108: /// Purpose
109: /// =======
110: ///
111: /// DGEHRD reduces a real general matrix A to upper Hessenberg form H by
112: /// an orthogonal similarity transformation: Q' * A * Q = H .
113: ///
114: ///</summary>
115: /// <param name="N">
116: /// (input) INTEGER
117: /// The order of the matrix A. N .GE. 0.
118: ///</param>
119: /// <param name="ILO">
120: /// (input) INTEGER
121: ///</param>
122: /// <param name="IHI">
123: /// (input) INTEGER
124: /// It is assumed that A is already upper triangular in rows
125: /// and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
126: /// set by a previous call to DGEBAL; otherwise they should be
127: /// set to 1 and N respectively. See Further Details.
128: /// 1 .LE. ILO .LE. IHI .LE. N, if N .GT. 0; ILO=1 and IHI=0, if N=0.
129: ///</param>
130: /// <param name="A">
131: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
132: /// On entry, the N-by-N general matrix to be reduced.
133: /// On exit, the upper triangle and the first subdiagonal of A
134: /// are overwritten with the upper Hessenberg matrix H, and the
135: /// elements below the first subdiagonal, with the array TAU,
136: /// represent the orthogonal matrix Q as a product of elementary
137: /// reflectors. See Further Details.
138: ///</param>
139: /// <param name="LDA">
140: /// (input) INTEGER
141: /// The leading dimension of the array A. LDA .GE. max(1,N).
142: ///</param>
143: /// <param name="TAU">
144: /// (output) DOUBLE PRECISION array, dimension (N-1)
145: /// The scalar factors of the elementary reflectors (see Further
146: /// Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
147: /// zero.
148: ///</param>
149: /// <param name="WORK">
150: /// (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
151: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
152: ///</param>
153: /// <param name="LWORK">
154: /// (input) INTEGER
155: /// The length of the array WORK. LWORK .GE. max(1,N).
156: /// For optimum performance LWORK .GE. N*NB, where NB is the
157: /// optimal blocksize.
158: ///
159: /// If LWORK = -1, then a workspace query is assumed; the routine
160: /// only calculates the optimal size of the WORK array, returns
161: /// this value as the first entry of the WORK array, and no error
162: /// message related to LWORK is issued by XERBLA.
163: ///</param>
164: /// <param name="INFO">
165: /// (output) INTEGER
166: /// = 0: successful exit
167: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
168: ///</param>
169: public void Run(int N, int ILO, int IHI, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
170: , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
171: {
172:
173: #region Array Index Correction
174:
175: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
176:
177: #endregion
178:
179:
180: #region Prolog
181:
182: // *
183: // * -- LAPACK routine (version 3.1) --
184: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
185: // * November 2006
186: // *
187: // * .. Scalar Arguments ..
188: // * ..
189: // * .. Array Arguments ..
190: // * ..
191: // *
192: // * Purpose
193: // * =======
194: // *
195: // * DGEHRD reduces a real general matrix A to upper Hessenberg form H by
196: // * an orthogonal similarity transformation: Q' * A * Q = H .
197: // *
198: // * Arguments
199: // * =========
200: // *
201: // * N (input) INTEGER
202: // * The order of the matrix A. N >= 0.
203: // *
204: // * ILO (input) INTEGER
205: // * IHI (input) INTEGER
206: // * It is assumed that A is already upper triangular in rows
207: // * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
208: // * set by a previous call to DGEBAL; otherwise they should be
209: // * set to 1 and N respectively. See Further Details.
210: // * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
211: // *
212: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
213: // * On entry, the N-by-N general matrix to be reduced.
214: // * On exit, the upper triangle and the first subdiagonal of A
215: // * are overwritten with the upper Hessenberg matrix H, and the
216: // * elements below the first subdiagonal, with the array TAU,
217: // * represent the orthogonal matrix Q as a product of elementary
218: // * reflectors. See Further Details.
219: // *
220: // * LDA (input) INTEGER
221: // * The leading dimension of the array A. LDA >= max(1,N).
222: // *
223: // * TAU (output) DOUBLE PRECISION array, dimension (N-1)
224: // * The scalar factors of the elementary reflectors (see Further
225: // * Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
226: // * zero.
227: // *
228: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
229: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
230: // *
231: // * LWORK (input) INTEGER
232: // * The length of the array WORK. LWORK >= max(1,N).
233: // * For optimum performance LWORK >= N*NB, where NB is the
234: // * optimal blocksize.
235: // *
236: // * If LWORK = -1, then a workspace query is assumed; the routine
237: // * only calculates the optimal size of the WORK array, returns
238: // * this value as the first entry of the WORK array, and no error
239: // * message related to LWORK is issued by XERBLA.
240: // *
241: // * INFO (output) INTEGER
242: // * = 0: successful exit
243: // * < 0: if INFO = -i, the i-th argument had an illegal value.
244: // *
245: // * Further Details
246: // * ===============
247: // *
248: // * The matrix Q is represented as a product of (ihi-ilo) elementary
249: // * reflectors
250: // *
251: // * Q = H(ilo) H(ilo+1) . . . H(ihi-1).
252: // *
253: // * Each H(i) has the form
254: // *
255: // * H(i) = I - tau * v * v'
256: // *
257: // * where tau is a real scalar, and v is a real vector with
258: // * v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
259: // * exit in A(i+2:ihi,i), and tau in TAU(i).
260: // *
261: // * The contents of A are illustrated by the following example, with
262: // * n = 7, ilo = 2 and ihi = 6:
263: // *
264: // * on entry, on exit,
265: // *
266: // * ( a a a a a a a ) ( a a h h h h a )
267: // * ( a a a a a a ) ( a h h h h a )
268: // * ( a a a a a a ) ( h h h h h h )
269: // * ( a a a a a a ) ( v2 h h h h h )
270: // * ( a a a a a a ) ( v2 v3 h h h h )
271: // * ( a a a a a a ) ( v2 v3 v4 h h h )
272: // * ( a ) ( a )
273: // *
274: // * where a denotes an element of the original matrix A, h denotes a
275: // * modified element of the upper Hessenberg matrix H, and vi denotes an
276: // * element of the vector defining H(i).
277: // *
278: // * This file is a slight modification of LAPACK-3.0's DGEHRD
279: // * subroutine incorporating improvements proposed by Quintana-Orti and
280: // * Van de Geijn (2005).
281: // *
282: // * =====================================================================
283: // *
284: // * .. Parameters ..
285: // * ..
286: // * .. Local Scalars ..
287: // * ..
288: // * .. Local Arrays ..
289: // * ..
290: // * .. External Subroutines ..
291: // * ..
292: // * .. Intrinsic Functions ..
293: // INTRINSIC MAX, MIN;
294: // * ..
295: // * .. External Functions ..
296: // * ..
297: // * .. Executable Statements ..
298: // *
299: // * Test the input parameters
300: // *
301:
302: #endregion
303:
304:
305: #region Body
306:
307: INFO = 0;
308: NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DGEHRD", " ", N, ILO, IHI, - 1));
309: LWKOPT = N * NB;
310: WORK[1 + o_work] = LWKOPT;
311: LQUERY = (LWORK == - 1);
312: if (N < 0)
313: {
314: INFO = - 1;
315: }
316: else
317: {
318: if (ILO < 1 || ILO > Math.Max(1, N))
319: {
320: INFO = - 2;
321: }
322: else
323: {
324: if (IHI < Math.Min(ILO, N) || IHI > N)
325: {
326: INFO = - 3;
327: }
328: else
329: {
330: if (LDA < Math.Max(1, N))
331: {
332: INFO = - 5;
333: }
334: else
335: {
336: if (LWORK < Math.Max(1, N) && !LQUERY)
337: {
338: INFO = - 8;
339: }
340: }
341: }
342: }
343: }
344: if (INFO != 0)
345: {
346: this._xerbla.Run("DGEHRD", - INFO);
347: return;
348: }
349: else
350: {
351: if (LQUERY)
352: {
353: return;
354: }
355: }
356: // *
357: // * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
358: // *
359: for (I = 1; I <= ILO - 1; I++)
360: {
361: TAU[I + o_tau] = ZERO;
362: }
363: for (I = Math.Max(1, IHI); I <= N - 1; I++)
364: {
365: TAU[I + o_tau] = ZERO;
366: }
367: // *
368: // * Quick return if possible
369: // *
370: NH = IHI - ILO + 1;
371: if (NH <= 1)
372: {
373: WORK[1 + o_work] = 1;
374: return;
375: }
376: // *
377: // * Determine the block size
378: // *
379: NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DGEHRD", " ", N, ILO, IHI, - 1));
380: NBMIN = 2;
381: IWS = 1;
382: if (NB > 1 && NB < NH)
383: {
384: // *
385: // * Determine when to cross over from blocked to unblocked code
386: // * (last block is always handled by unblocked code)
387: // *
388: NX = Math.Max(NB, this._ilaenv.Run(3, "DGEHRD", " ", N, ILO, IHI, - 1));
389: if (NX < NH)
390: {
391: // *
392: // * Determine if workspace is large enough for blocked code
393: // *
394: IWS = N * NB;
395: if (LWORK < IWS)
396: {
397: // *
398: // * Not enough workspace to use optimal NB: determine the
399: // * minimum value of NB, and reduce NB or force use of
400: // * unblocked code
401: // *
402: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGEHRD", " ", N, ILO, IHI, - 1));
403: if (LWORK >= N * NBMIN)
404: {
405: NB = LWORK / N;
406: }
407: else
408: {
409: NB = 1;
410: }
411: }
412: }
413: }
414: LDWORK = N;
415: // *
416: if (NB < NBMIN || NB >= NH)
417: {
418: // *
419: // * Use unblocked code below
420: // *
421: I = ILO;
422: // *
423: }
424: else
425: {
426: // *
427: // * Use blocked code
428: // *
429: for (I = ILO; (NB >= 0) ? (I <= IHI - 1 - NX) : (I >= IHI - 1 - NX); I += NB)
430: {
431: IB = Math.Min(NB, IHI - I);
432: // *
433: // * Reduce columns i:i+ib-1 to Hessenberg form, returning the
434: // * matrices V and T of the block reflector H = I - V*T*V'
435: // * which performs the reduction, and also the matrix Y = A*V*T
436: // *
437: this._dlahr2.Run(IHI, I, IB, ref A, 1+I * LDA + o_a, LDA, ref TAU, I + o_tau
438: , ref T, offset_t, LDT, ref WORK, offset_work, LDWORK);
439: // *
440: // * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
441: // * right, computing A := A - Y * V'. V(i+ib,ib-1) must be set
442: // * to 1
443: // *
444: EI = A[I + IB+(I + IB - 1) * LDA + o_a];
445: A[I + IB+(I + IB - 1) * LDA + o_a] = ONE;
446: this._dgemm.Run("No transpose", "Transpose", IHI, IHI - I - IB + 1, IB, - ONE
447: , WORK, offset_work, LDWORK, A, I + IB+I * LDA + o_a, LDA, ONE, ref A, 1+(I + IB) * LDA + o_a
448: , LDA);
449: A[I + IB+(I + IB - 1) * LDA + o_a] = EI;
450: // *
451: // * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
452: // * right
453: // *
454: this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", I, IB - 1
455: , ONE, A, I + 1+I * LDA + o_a, LDA, ref WORK, offset_work, LDWORK);
456: for (J = 0; J <= IB - 2; J++)
457: {
458: this._daxpy.Run(I, - ONE, WORK, LDWORK * J + 1 + o_work, 1, ref A, 1+(I + J + 1) * LDA + o_a, 1);
459: }
460: // *
461: // * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
462: // * left
463: // *
464: this._dlarfb.Run("Left", "Transpose", "Forward", "Columnwise", IHI - I, N - I - IB + 1
465: , IB, A, I + 1+I * LDA + o_a, LDA, T, offset_t, LDT, ref A, I + 1+(I + IB) * LDA + o_a
466: , LDA, ref WORK, offset_work, LDWORK);
467: }
468: }
469: // *
470: // * Use unblocked code to reduce the rest of the matrix
471: // *
472: this._dgehd2.Run(N, I, IHI, ref A, offset_a, LDA, ref TAU, offset_tau
473: , ref WORK, offset_work, ref IINFO);
474: WORK[1 + o_work] = IWS;
475: // *
476: return;
477: // *
478: // * End of DGEHRD
479: // *
480:
481: #endregion
482:
483: }
484: }
485: }