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