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: /// DLARFT forms the triangular factor T of a real block reflector H
27: /// of order n, which is defined as a product of k elementary reflectors.
28: ///
29: /// If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
30: ///
31: /// If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
32: ///
33: /// If STOREV = 'C', the vector which defines the elementary reflector
34: /// H(i) is stored in the i-th column of the array V, and
35: ///
36: /// H = I - V * T * V'
37: ///
38: /// If STOREV = 'R', the vector which defines the elementary reflector
39: /// H(i) is stored in the i-th row of the array V, and
40: ///
41: /// H = I - V' * T * V
42: ///
43: ///</summary>
44: public class DLARFT
45: {
46:
47:
48: #region Dependencies
49:
50: DGEMV _dgemv; DTRMV _dtrmv; LSAME _lsame;
51:
52: #endregion
53:
54:
55: #region Fields
56:
57: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int I = 0; int J = 0; double VII = 0;
58:
59: #endregion
60:
61: public DLARFT(DGEMV dgemv, DTRMV dtrmv, LSAME lsame)
62: {
63:
64:
65: #region Set Dependencies
66:
67: this._dgemv = dgemv; this._dtrmv = dtrmv; this._lsame = lsame;
68:
69: #endregion
70:
71: }
72:
73: public DLARFT()
74: {
75:
76:
77: #region Dependencies (Initialization)
78:
79: LSAME lsame = new LSAME();
80: XERBLA xerbla = new XERBLA();
81: DGEMV dgemv = new DGEMV(lsame, xerbla);
82: DTRMV dtrmv = new DTRMV(lsame, xerbla);
83:
84: #endregion
85:
86:
87: #region Set Dependencies
88:
89: this._dgemv = dgemv; this._dtrmv = dtrmv; this._lsame = lsame;
90:
91: #endregion
92:
93: }
94: /// <summary>
95: /// Purpose
96: /// =======
97: ///
98: /// DLARFT forms the triangular factor T of a real block reflector H
99: /// of order n, which is defined as a product of k elementary reflectors.
100: ///
101: /// If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
102: ///
103: /// If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
104: ///
105: /// If STOREV = 'C', the vector which defines the elementary reflector
106: /// H(i) is stored in the i-th column of the array V, and
107: ///
108: /// H = I - V * T * V'
109: ///
110: /// If STOREV = 'R', the vector which defines the elementary reflector
111: /// H(i) is stored in the i-th row of the array V, and
112: ///
113: /// H = I - V' * T * V
114: ///
115: ///</summary>
116: /// <param name="DIRECT">
117: /// (input) CHARACTER*1
118: /// Specifies the order in which the elementary reflectors are
119: /// multiplied to form the block reflector:
120: /// = 'F': H = H(1) H(2) . . . H(k) (Forward)
121: /// = 'B': H = H(k) . . . H(2) H(1) (Backward)
122: ///</param>
123: /// <param name="STOREV">
124: /// (input) CHARACTER*1
125: /// Specifies how the vectors which define the elementary
126: /// reflectors are stored (see also Further Details):
127: /// = 'C': columnwise
128: /// = 'R': rowwise
129: ///</param>
130: /// <param name="N">
131: /// (input) INTEGER
132: /// The order of the block reflector H. N .GE. 0.
133: ///</param>
134: /// <param name="K">
135: /// (input) INTEGER
136: /// The order of the triangular factor T (= the number of
137: /// elementary reflectors). K .GE. 1.
138: ///</param>
139: /// <param name="V">
140: /// (input/output) DOUBLE PRECISION array, dimension
141: /// (LDV,K) if STOREV = 'C'
142: /// (LDV,N) if STOREV = 'R'
143: /// The matrix V. See further details.
144: ///</param>
145: /// <param name="LDV">
146: /// (input) INTEGER
147: /// The leading dimension of the array V.
148: /// If STOREV = 'C', LDV .GE. max(1,N); if STOREV = 'R', LDV .GE. K.
149: ///</param>
150: /// <param name="TAU">
151: /// (input) DOUBLE PRECISION array, dimension (K)
152: /// TAU(i) must contain the scalar factor of the elementary
153: /// reflector H(i).
154: ///</param>
155: /// <param name="T">
156: /// (output) DOUBLE PRECISION array, dimension (LDT,K)
157: /// The k by k triangular factor T of the block reflector.
158: /// If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
159: /// lower triangular. The rest of the array is not used.
160: ///</param>
161: /// <param name="LDT">
162: /// (input) INTEGER
163: /// The leading dimension of the array T. LDT .GE. K.
164: ///</param>
165: public void Run(string DIRECT, string STOREV, int N, int K, ref double[] V, int offset_v, int LDV
166: , double[] TAU, int offset_tau, ref double[] T, int offset_t, int LDT)
167: {
168:
169: #region Array Index Correction
170:
171: int o_v = -1 - LDV + offset_v; int o_tau = -1 + offset_tau; int o_t = -1 - LDT + offset_t;
172:
173: #endregion
174:
175:
176: #region Strings
177:
178: DIRECT = DIRECT.Substring(0, 1); STOREV = STOREV.Substring(0, 1);
179:
180: #endregion
181:
182:
183: #region Prolog
184:
185: // *
186: // * -- LAPACK auxiliary routine (version 3.1) --
187: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
188: // * November 2006
189: // *
190: // * .. Scalar Arguments ..
191: // * ..
192: // * .. Array Arguments ..
193: // * ..
194: // *
195: // * Purpose
196: // * =======
197: // *
198: // * DLARFT forms the triangular factor T of a real block reflector H
199: // * of order n, which is defined as a product of k elementary reflectors.
200: // *
201: // * If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
202: // *
203: // * If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
204: // *
205: // * If STOREV = 'C', the vector which defines the elementary reflector
206: // * H(i) is stored in the i-th column of the array V, and
207: // *
208: // * H = I - V * T * V'
209: // *
210: // * If STOREV = 'R', the vector which defines the elementary reflector
211: // * H(i) is stored in the i-th row of the array V, and
212: // *
213: // * H = I - V' * T * V
214: // *
215: // * Arguments
216: // * =========
217: // *
218: // * DIRECT (input) CHARACTER*1
219: // * Specifies the order in which the elementary reflectors are
220: // * multiplied to form the block reflector:
221: // * = 'F': H = H(1) H(2) . . . H(k) (Forward)
222: // * = 'B': H = H(k) . . . H(2) H(1) (Backward)
223: // *
224: // * STOREV (input) CHARACTER*1
225: // * Specifies how the vectors which define the elementary
226: // * reflectors are stored (see also Further Details):
227: // * = 'C': columnwise
228: // * = 'R': rowwise
229: // *
230: // * N (input) INTEGER
231: // * The order of the block reflector H. N >= 0.
232: // *
233: // * K (input) INTEGER
234: // * The order of the triangular factor T (= the number of
235: // * elementary reflectors). K >= 1.
236: // *
237: // * V (input/output) DOUBLE PRECISION array, dimension
238: // * (LDV,K) if STOREV = 'C'
239: // * (LDV,N) if STOREV = 'R'
240: // * The matrix V. See further details.
241: // *
242: // * LDV (input) INTEGER
243: // * The leading dimension of the array V.
244: // * If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
245: // *
246: // * TAU (input) DOUBLE PRECISION array, dimension (K)
247: // * TAU(i) must contain the scalar factor of the elementary
248: // * reflector H(i).
249: // *
250: // * T (output) DOUBLE PRECISION array, dimension (LDT,K)
251: // * The k by k triangular factor T of the block reflector.
252: // * If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
253: // * lower triangular. The rest of the array is not used.
254: // *
255: // * LDT (input) INTEGER
256: // * The leading dimension of the array T. LDT >= K.
257: // *
258: // * Further Details
259: // * ===============
260: // *
261: // * The shape of the matrix V and the storage of the vectors which define
262: // * the H(i) is best illustrated by the following example with n = 5 and
263: // * k = 3. The elements equal to 1 are not stored; the corresponding
264: // * array elements are modified but restored on exit. The rest of the
265: // * array is not used.
266: // *
267: // * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
268: // *
269: // * V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
270: // * ( v1 1 ) ( 1 v2 v2 v2 )
271: // * ( v1 v2 1 ) ( 1 v3 v3 )
272: // * ( v1 v2 v3 )
273: // * ( v1 v2 v3 )
274: // *
275: // * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
276: // *
277: // * V = ( v1 v2 v3 ) V = ( v1 v1 1 )
278: // * ( v1 v2 v3 ) ( v2 v2 v2 1 )
279: // * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
280: // * ( 1 v3 )
281: // * ( 1 )
282: // *
283: // * =====================================================================
284: // *
285: // * .. Parameters ..
286: // * ..
287: // * .. Local Scalars ..
288: // * ..
289: // * .. External Subroutines ..
290: // * ..
291: // * .. External Functions ..
292: // * ..
293: // * .. Executable Statements ..
294: // *
295: // * Quick return if possible
296: // *
297:
298: #endregion
299:
300:
301: #region Body
302:
303: if (N == 0) return;
304: // *
305: if (this._lsame.Run(DIRECT, "F"))
306: {
307: for (I = 1; I <= K; I++)
308: {
309: if (TAU[I + o_tau] == ZERO)
310: {
311: // *
312: // * H(i) = I
313: // *
314: for (J = 1; J <= I; J++)
315: {
316: T[J+I * LDT + o_t] = ZERO;
317: }
318: }
319: else
320: {
321: // *
322: // * general case
323: // *
324: VII = V[I+I * LDV + o_v];
325: V[I+I * LDV + o_v] = ONE;
326: if (this._lsame.Run(STOREV, "C"))
327: {
328: // *
329: // * T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i)
330: // *
331: this._dgemv.Run("Transpose", N - I + 1, I - 1, - TAU[I + o_tau], V, I+1 * LDV + o_v, LDV
332: , V, I+I * LDV + o_v, 1, ZERO, ref T, 1+I * LDT + o_t, 1);
333: }
334: else
335: {
336: // *
337: // * T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)'
338: // *
339: this._dgemv.Run("No transpose", I - 1, N - I + 1, - TAU[I + o_tau], V, 1+I * LDV + o_v, LDV
340: , V, I+I * LDV + o_v, LDV, ZERO, ref T, 1+I * LDT + o_t, 1);
341: }
342: V[I+I * LDV + o_v] = VII;
343: // *
344: // * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
345: // *
346: this._dtrmv.Run("Upper", "No transpose", "Non-unit", I - 1, T, offset_t, LDT
347: , ref T, 1+I * LDT + o_t, 1);
348: T[I+I * LDT + o_t] = TAU[I + o_tau];
349: }
350: }
351: }
352: else
353: {
354: for (I = K; I >= 1; I += - 1)
355: {
356: if (TAU[I + o_tau] == ZERO)
357: {
358: // *
359: // * H(i) = I
360: // *
361: for (J = I; J <= K; J++)
362: {
363: T[J+I * LDT + o_t] = ZERO;
364: }
365: }
366: else
367: {
368: // *
369: // * general case
370: // *
371: if (I < K)
372: {
373: if (this._lsame.Run(STOREV, "C"))
374: {
375: VII = V[N - K + I+I * LDV + o_v];
376: V[N - K + I+I * LDV + o_v] = ONE;
377: // *
378: // * T(i+1:k,i) :=
379: // * - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i)
380: // *
381: this._dgemv.Run("Transpose", N - K + I, K - I, - TAU[I + o_tau], V, 1+(I + 1) * LDV + o_v, LDV
382: , V, 1+I * LDV + o_v, 1, ZERO, ref T, I + 1+I * LDT + o_t, 1);
383: V[N - K + I+I * LDV + o_v] = VII;
384: }
385: else
386: {
387: VII = V[I+(N - K + I) * LDV + o_v];
388: V[I+(N - K + I) * LDV + o_v] = ONE;
389: // *
390: // * T(i+1:k,i) :=
391: // * - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)'
392: // *
393: this._dgemv.Run("No transpose", K - I, N - K + I, - TAU[I + o_tau], V, I + 1+1 * LDV + o_v, LDV
394: , V, I+1 * LDV + o_v, LDV, ZERO, ref T, I + 1+I * LDT + o_t, 1);
395: V[I+(N - K + I) * LDV + o_v] = VII;
396: }
397: // *
398: // * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
399: // *
400: this._dtrmv.Run("Lower", "No transpose", "Non-unit", K - I, T, I + 1+(I + 1) * LDT + o_t, LDT
401: , ref T, I + 1+I * LDT + o_t, 1);
402: }
403: T[I+I * LDT + o_t] = TAU[I + o_tau];
404: }
405: }
406: }
407: return;
408: // *
409: // * End of DLARFT
410: // *
411:
412: #endregion
413:
414: }
415: }
416: }