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 auxiliary routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
27: /// matrix A so that elements below the k-th subdiagonal are zero. The
28: /// reduction is performed by an orthogonal similarity transformation
29: /// Q' * A * Q. The routine returns the matrices V and T which determine
30: /// Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
31: ///
32: /// This is an auxiliary routine called by DGEHRD.
33: ///
34: ///</summary>
35: public class DLAHR2
36: {
37:
38:
39: #region Dependencies
40:
41: DAXPY _daxpy; DCOPY _dcopy; DGEMM _dgemm; DGEMV _dgemv; DLACPY _dlacpy; DLARFG _dlarfg; DSCAL _dscal; DTRMM _dtrmm;
42: DTRMV _dtrmv;
43:
44: #endregion
45:
46:
47: #region Fields
48:
49: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; double EI = 0;
50:
51: #endregion
52:
53: public DLAHR2(DAXPY daxpy, DCOPY dcopy, DGEMM dgemm, DGEMV dgemv, DLACPY dlacpy, DLARFG dlarfg, DSCAL dscal, DTRMM dtrmm, DTRMV dtrmv)
54: {
55:
56:
57: #region Set Dependencies
58:
59: this._daxpy = daxpy; this._dcopy = dcopy; this._dgemm = dgemm; this._dgemv = dgemv; this._dlacpy = dlacpy;
60: this._dlarfg = dlarfg;this._dscal = dscal; this._dtrmm = dtrmm; this._dtrmv = dtrmv;
61:
62: #endregion
63:
64: }
65:
66: public DLAHR2()
67: {
68:
69:
70: #region Dependencies (Initialization)
71:
72: DAXPY daxpy = new DAXPY();
73: DCOPY dcopy = new DCOPY();
74: LSAME lsame = new LSAME();
75: XERBLA xerbla = new XERBLA();
76: DLAMC3 dlamc3 = new DLAMC3();
77: DLAPY2 dlapy2 = new DLAPY2();
78: DNRM2 dnrm2 = new DNRM2();
79: DSCAL dscal = new DSCAL();
80: DGEMM dgemm = new DGEMM(lsame, xerbla);
81: DGEMV dgemv = new DGEMV(lsame, xerbla);
82: DLACPY dlacpy = new DLACPY(lsame);
83: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
84: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
85: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
86: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
87: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
88: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
89: DTRMM dtrmm = new DTRMM(lsame, xerbla);
90: DTRMV dtrmv = new DTRMV(lsame, xerbla);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._daxpy = daxpy; this._dcopy = dcopy; this._dgemm = dgemm; this._dgemv = dgemv; this._dlacpy = dlacpy;
98: this._dlarfg = dlarfg;this._dscal = dscal; this._dtrmm = dtrmm; this._dtrmv = dtrmv;
99:
100: #endregion
101:
102: }
103: /// <summary>
104: /// Purpose
105: /// =======
106: ///
107: /// DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
108: /// matrix A so that elements below the k-th subdiagonal are zero. The
109: /// reduction is performed by an orthogonal similarity transformation
110: /// Q' * A * Q. The routine returns the matrices V and T which determine
111: /// Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
112: ///
113: /// This is an auxiliary routine called by DGEHRD.
114: ///
115: ///</summary>
116: /// <param name="N">
117: /// (input) INTEGER
118: /// The order of the matrix A.
119: ///</param>
120: /// <param name="K">
121: /// (input) INTEGER
122: /// The offset for the reduction. Elements below the k-th
123: /// subdiagonal in the first NB columns are reduced to zero.
124: /// K .LT. N.
125: ///</param>
126: /// <param name="NB">
127: /// (input) INTEGER
128: /// The number of columns to be reduced.
129: ///</param>
130: /// <param name="A">
131: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)
132: /// On entry, the n-by-(n-k+1) general matrix A.
133: /// On exit, the elements on and above the k-th subdiagonal in
134: /// the first NB columns are overwritten with the corresponding
135: /// elements of the reduced matrix; the elements below the k-th
136: /// subdiagonal, with the array TAU, represent the matrix Q as a
137: /// product of elementary reflectors. The other columns of A are
138: /// unchanged. See Further Details.
139: ///</param>
140: /// <param name="LDA">
141: /// (input) INTEGER
142: /// The leading dimension of the array A. LDA .GE. max(1,N).
143: ///</param>
144: /// <param name="TAU">
145: /// (output) DOUBLE PRECISION array, dimension (NB)
146: /// The scalar factors of the elementary reflectors. See Further
147: /// Details.
148: ///</param>
149: /// <param name="T">
150: /// (output) DOUBLE PRECISION array, dimension (LDT,NB)
151: /// The upper triangular matrix T.
152: ///</param>
153: /// <param name="LDT">
154: /// (input) INTEGER
155: /// The leading dimension of the array T. LDT .GE. NB.
156: ///</param>
157: /// <param name="Y">
158: /// (output) DOUBLE PRECISION array, dimension (LDY,NB)
159: /// The n-by-nb matrix Y.
160: ///</param>
161: /// <param name="LDY">
162: /// (input) INTEGER
163: /// The leading dimension of the array Y. LDY .GE. N.
164: ///</param>
165: public void Run(int N, int K, int NB, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
166: , ref double[] T, int offset_t, int LDT, ref double[] Y, int offset_y, int LDY)
167: {
168:
169: #region Array Index Correction
170:
171: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_t = -1 - LDT + offset_t;
172: int o_y = -1 - LDY + offset_y;
173:
174: #endregion
175:
176:
177: #region Prolog
178:
179: // *
180: // * -- LAPACK auxiliary routine (version 3.1) --
181: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
182: // * November 2006
183: // *
184: // * .. Scalar Arguments ..
185: // * ..
186: // * .. Array Arguments ..
187: // * ..
188: // *
189: // * Purpose
190: // * =======
191: // *
192: // * DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)
193: // * matrix A so that elements below the k-th subdiagonal are zero. The
194: // * reduction is performed by an orthogonal similarity transformation
195: // * Q' * A * Q. The routine returns the matrices V and T which determine
196: // * Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
197: // *
198: // * This is an auxiliary routine called by DGEHRD.
199: // *
200: // * Arguments
201: // * =========
202: // *
203: // * N (input) INTEGER
204: // * The order of the matrix A.
205: // *
206: // * K (input) INTEGER
207: // * The offset for the reduction. Elements below the k-th
208: // * subdiagonal in the first NB columns are reduced to zero.
209: // * K < N.
210: // *
211: // * NB (input) INTEGER
212: // * The number of columns to be reduced.
213: // *
214: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)
215: // * On entry, the n-by-(n-k+1) general matrix A.
216: // * On exit, the elements on and above the k-th subdiagonal in
217: // * the first NB columns are overwritten with the corresponding
218: // * elements of the reduced matrix; the elements below the k-th
219: // * subdiagonal, with the array TAU, represent the matrix Q as a
220: // * product of elementary reflectors. The other columns of A are
221: // * unchanged. See Further Details.
222: // *
223: // * LDA (input) INTEGER
224: // * The leading dimension of the array A. LDA >= max(1,N).
225: // *
226: // * TAU (output) DOUBLE PRECISION array, dimension (NB)
227: // * The scalar factors of the elementary reflectors. See Further
228: // * Details.
229: // *
230: // * T (output) DOUBLE PRECISION array, dimension (LDT,NB)
231: // * The upper triangular matrix T.
232: // *
233: // * LDT (input) INTEGER
234: // * The leading dimension of the array T. LDT >= NB.
235: // *
236: // * Y (output) DOUBLE PRECISION array, dimension (LDY,NB)
237: // * The n-by-nb matrix Y.
238: // *
239: // * LDY (input) INTEGER
240: // * The leading dimension of the array Y. LDY >= N.
241: // *
242: // * Further Details
243: // * ===============
244: // *
245: // * The matrix Q is represented as a product of nb elementary reflectors
246: // *
247: // * Q = H(1) H(2) . . . H(nb).
248: // *
249: // * Each H(i) has the form
250: // *
251: // * H(i) = I - tau * v * v'
252: // *
253: // * where tau is a real scalar, and v is a real vector with
254: // * v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
255: // * A(i+k+1:n,i), and tau in TAU(i).
256: // *
257: // * The elements of the vectors v together form the (n-k+1)-by-nb matrix
258: // * V which is needed, with T and Y, to apply the transformation to the
259: // * unreduced part of the matrix, using an update of the form:
260: // * A := (I - V*T*V') * (A - Y*V').
261: // *
262: // * The contents of A on exit are illustrated by the following example
263: // * with n = 7, k = 3 and nb = 2:
264: // *
265: // * ( a a a a a )
266: // * ( a a a a a )
267: // * ( a a a a a )
268: // * ( h h a a a )
269: // * ( v1 h a a a )
270: // * ( v1 v2 a a a )
271: // * ( v1 v2 a a a )
272: // *
273: // * where a denotes an element of the original matrix A, h denotes a
274: // * modified element of the upper Hessenberg matrix H, and vi denotes an
275: // * element of the vector defining H(i).
276: // *
277: // * This file is a slight modification of LAPACK-3.0's DLAHRD
278: // * incorporating improvements proposed by Quintana-Orti and Van de
279: // * Gejin. Note that the entries of A(1:K,2:NB) differ from those
280: // * returned by the original LAPACK routine. This function is
281: // * not backward compatible with LAPACK3.0.
282: // *
283: // * =====================================================================
284: // *
285: // * .. Parameters ..
286: // * ..
287: // * .. Local Scalars ..
288: // * ..
289: // * .. External Subroutines ..
290: // * ..
291: // * .. Intrinsic Functions ..
292: // INTRINSIC MIN;
293: // * ..
294: // * .. Executable Statements ..
295: // *
296: // * Quick return if possible
297: // *
298:
299: #endregion
300:
301:
302: #region Body
303:
304: if (N <= 1) return;
305: // *
306: for (I = 1; I <= NB; I++)
307: {
308: if (I > 1)
309: {
310: // *
311: // * Update A(K+1:N,I)
312: // *
313: // * Update I-th column of A - Y * V'
314: // *
315: this._dgemv.Run("NO TRANSPOSE", N - K, I - 1, - ONE, Y, K + 1+1 * LDY + o_y, LDY
316: , A, K + I - 1+1 * LDA + o_a, LDA, ONE, ref A, K + 1+I * LDA + o_a, 1);
317: // *
318: // * Apply I - V * T' * V' to this column (call it b) from the
319: // * left, using the last column of T as workspace
320: // *
321: // * Let V = ( V1 ) and b = ( b1 ) (first I-1 rows)
322: // * ( V2 ) ( b2 )
323: // *
324: // * where V1 is unit lower triangular
325: // *
326: // * w := V1' * b1
327: // *
328: this._dcopy.Run(I - 1, A, K + 1+I * LDA + o_a, 1, ref T, 1+NB * LDT + o_t, 1);
329: this._dtrmv.Run("Lower", "Transpose", "UNIT", I - 1, A, K + 1+1 * LDA + o_a, LDA
330: , ref T, 1+NB * LDT + o_t, 1);
331: // *
332: // * w := w + V2'*b2
333: // *
334: this._dgemv.Run("Transpose", N - K - I + 1, I - 1, ONE, A, K + I+1 * LDA + o_a, LDA
335: , A, K + I+I * LDA + o_a, 1, ONE, ref T, 1+NB * LDT + o_t, 1);
336: // *
337: // * w := T'*w
338: // *
339: this._dtrmv.Run("Upper", "Transpose", "NON-UNIT", I - 1, T, offset_t, LDT
340: , ref T, 1+NB * LDT + o_t, 1);
341: // *
342: // * b2 := b2 - V2*w
343: // *
344: this._dgemv.Run("NO TRANSPOSE", N - K - I + 1, I - 1, - ONE, A, K + I+1 * LDA + o_a, LDA
345: , T, 1+NB * LDT + o_t, 1, ONE, ref A, K + I+I * LDA + o_a, 1);
346: // *
347: // * b1 := b1 - V1*w
348: // *
349: this._dtrmv.Run("Lower", "NO TRANSPOSE", "UNIT", I - 1, A, K + 1+1 * LDA + o_a, LDA
350: , ref T, 1+NB * LDT + o_t, 1);
351: this._daxpy.Run(I - 1, - ONE, T, 1+NB * LDT + o_t, 1, ref A, K + 1+I * LDA + o_a, 1);
352: // *
353: A[K + I - 1+(I - 1) * LDA + o_a] = EI;
354: }
355: // *
356: // * Generate the elementary reflector H(I) to annihilate
357: // * A(K+I+1:N,I)
358: // *
359: this._dlarfg.Run(N - K - I + 1, ref A[K + I+I * LDA + o_a], ref A, Math.Min(K + I + 1, N)+I * LDA + o_a, 1, ref TAU[I + o_tau]);
360: EI = A[K + I+I * LDA + o_a];
361: A[K + I+I * LDA + o_a] = ONE;
362: // *
363: // * Compute Y(K+1:N,I)
364: // *
365: this._dgemv.Run("NO TRANSPOSE", N - K, N - K - I + 1, ONE, A, K + 1+(I + 1) * LDA + o_a, LDA
366: , A, K + I+I * LDA + o_a, 1, ZERO, ref Y, K + 1+I * LDY + o_y, 1);
367: this._dgemv.Run("Transpose", N - K - I + 1, I - 1, ONE, A, K + I+1 * LDA + o_a, LDA
368: , A, K + I+I * LDA + o_a, 1, ZERO, ref T, 1+I * LDT + o_t, 1);
369: this._dgemv.Run("NO TRANSPOSE", N - K, I - 1, - ONE, Y, K + 1+1 * LDY + o_y, LDY
370: , T, 1+I * LDT + o_t, 1, ONE, ref Y, K + 1+I * LDY + o_y, 1);
371: this._dscal.Run(N - K, TAU[I + o_tau], ref Y, K + 1+I * LDY + o_y, 1);
372: // *
373: // * Compute T(1:I,I)
374: // *
375: this._dscal.Run(I - 1, - TAU[I + o_tau], ref T, 1+I * LDT + o_t, 1);
376: this._dtrmv.Run("Upper", "No Transpose", "NON-UNIT", I - 1, T, offset_t, LDT
377: , ref T, 1+I * LDT + o_t, 1);
378: T[I+I * LDT + o_t] = TAU[I + o_tau];
379: // *
380: }
381: A[K + NB+NB * LDA + o_a] = EI;
382: // *
383: // * Compute Y(1:K,1:NB)
384: // *
385: this._dlacpy.Run("ALL", K, NB, A, 1+2 * LDA + o_a, LDA, ref Y, offset_y
386: , LDY);
387: this._dtrmm.Run("RIGHT", "Lower", "NO TRANSPOSE", "UNIT", K, NB
388: , ONE, A, K + 1+1 * LDA + o_a, LDA, ref Y, offset_y, LDY);
389: if (N > K + NB)
390: {
391: this._dgemm.Run("NO TRANSPOSE", "NO TRANSPOSE", K, NB, N - K - NB, ONE
392: , A, 1+(2 + NB) * LDA + o_a, LDA, A, K + 1 + NB+1 * LDA + o_a, LDA, ONE, ref Y, offset_y
393: , LDY);
394: }
395: this._dtrmm.Run("RIGHT", "Upper", "NO TRANSPOSE", "NON-UNIT", K, NB
396: , ONE, T, offset_t, LDT, ref Y, offset_y, LDY);
397: // *
398: return;
399: // *
400: // * End of DLAHR2
401: // *
402:
403: #endregion
404:
405: }
406: }
407: }