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: /// DLARZB applies a real block reflector H or its transpose H**T to
27: /// a real distributed M-by-N C from the left or the right.
28: ///
29: /// Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
30: ///
31: ///</summary>
32: public class DLARZB
33: {
34:
35:
36: #region Dependencies
37:
38: LSAME _lsame; DCOPY _dcopy; DGEMM _dgemm; DTRMM _dtrmm; XERBLA _xerbla;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ONE = 1.0E+0; string TRANST = new string(' ', 1); int I = 0; int INFO = 0; int J = 0;
46:
47: #endregion
48:
49: public DLARZB(LSAME lsame, DCOPY dcopy, DGEMM dgemm, DTRMM dtrmm, XERBLA xerbla)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; this._xerbla = xerbla;
56:
57: #endregion
58:
59: }
60:
61: public DLARZB()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: LSAME lsame = new LSAME();
68: DCOPY dcopy = new DCOPY();
69: XERBLA xerbla = new XERBLA();
70: DGEMM dgemm = new DGEMM(lsame, xerbla);
71: DTRMM dtrmm = new DTRMM(lsame, xerbla);
72:
73: #endregion
74:
75:
76: #region Set Dependencies
77:
78: this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; this._xerbla = xerbla;
79:
80: #endregion
81:
82: }
83: /// <summary>
84: /// Purpose
85: /// =======
86: ///
87: /// DLARZB applies a real block reflector H or its transpose H**T to
88: /// a real distributed M-by-N C from the left or the right.
89: ///
90: /// Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
91: ///
92: ///</summary>
93: /// <param name="SIDE">
94: /// (input) CHARACTER*1
95: /// = 'L': apply H or H' from the Left
96: /// = 'R': apply H or H' from the Right
97: ///</param>
98: /// <param name="TRANS">
99: /// (input) CHARACTER*1
100: /// = 'N': apply H (No transpose)
101: /// = 'C': apply H' (Transpose)
102: ///</param>
103: /// <param name="DIRECT">
104: /// (input) CHARACTER*1
105: /// Indicates how H is formed from a product of elementary
106: /// reflectors
107: /// = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
108: /// = 'B': H = H(k) . . . H(2) H(1) (Backward)
109: ///</param>
110: /// <param name="STOREV">
111: /// (input) CHARACTER*1
112: /// Indicates how the vectors which define the elementary
113: /// reflectors are stored:
114: /// = 'C': Columnwise (not supported yet)
115: /// = 'R': Rowwise
116: ///</param>
117: /// <param name="M">
118: /// (input) INTEGER
119: /// The number of rows of the matrix C.
120: ///</param>
121: /// <param name="N">
122: /// (input) INTEGER
123: /// The number of columns of the matrix C.
124: ///</param>
125: /// <param name="K">
126: /// (input) INTEGER
127: /// The order of the matrix T (= the number of elementary
128: /// reflectors whose product defines the block reflector).
129: ///</param>
130: /// <param name="L">
131: /// (input) INTEGER
132: /// The number of columns of the matrix V containing the
133: /// meaningful part of the Householder reflectors.
134: /// If SIDE = 'L', M .GE. L .GE. 0, if SIDE = 'R', N .GE. L .GE. 0.
135: ///</param>
136: /// <param name="V">
137: /// (input) DOUBLE PRECISION array, dimension (LDV,NV).
138: /// If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
139: ///</param>
140: /// <param name="LDV">
141: /// (input) INTEGER
142: /// The leading dimension of the array V.
143: /// If STOREV = 'C', LDV .GE. L; if STOREV = 'R', LDV .GE. K.
144: ///</param>
145: /// <param name="T">
146: /// (input) DOUBLE PRECISION array, dimension (LDT,K)
147: /// The triangular K-by-K matrix T in the representation of the
148: /// block reflector.
149: ///</param>
150: /// <param name="LDT">
151: /// (input) INTEGER
152: /// The leading dimension of the array T. LDT .GE. K.
153: ///</param>
154: /// <param name="C">
155: /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
156: /// On entry, the M-by-N matrix C.
157: /// On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
158: ///</param>
159: /// <param name="LDC">
160: /// (input) INTEGER
161: /// The leading dimension of the array C. LDC .GE. max(1,M).
162: ///</param>
163: /// <param name="WORK">
164: /// (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
165: ///</param>
166: /// <param name="LDWORK">
167: /// (input) INTEGER
168: /// The leading dimension of the array WORK.
169: /// If SIDE = 'L', LDWORK .GE. max(1,N);
170: /// if SIDE = 'R', LDWORK .GE. max(1,M).
171: ///</param>
172: public void Run(string SIDE, string TRANS, string DIRECT, string STOREV, int M, int N
173: , int K, int L, double[] V, int offset_v, int LDV, double[] T, int offset_t, int LDT
174: , ref double[] C, int offset_c, int LDC, ref double[] WORK, int offset_work, int LDWORK)
175: {
176:
177: #region Array Index Correction
178:
179: int o_v = -1 - LDV + offset_v; int o_t = -1 - LDT + offset_t; int o_c = -1 - LDC + offset_c;
180: int o_work = -1 - LDWORK + offset_work;
181:
182: #endregion
183:
184:
185: #region Strings
186:
187: SIDE = SIDE.Substring(0, 1); TRANS = TRANS.Substring(0, 1); DIRECT = DIRECT.Substring(0, 1);
188: STOREV = STOREV.Substring(0, 1);
189:
190: #endregion
191:
192:
193: #region Prolog
194:
195: // *
196: // * -- LAPACK routine (version 3.1) --
197: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
198: // * November 2006
199: // *
200: // * .. Scalar Arguments ..
201: // * ..
202: // * .. Array Arguments ..
203: // * ..
204: // *
205: // * Purpose
206: // * =======
207: // *
208: // * DLARZB applies a real block reflector H or its transpose H**T to
209: // * a real distributed M-by-N C from the left or the right.
210: // *
211: // * Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
212: // *
213: // * Arguments
214: // * =========
215: // *
216: // * SIDE (input) CHARACTER*1
217: // * = 'L': apply H or H' from the Left
218: // * = 'R': apply H or H' from the Right
219: // *
220: // * TRANS (input) CHARACTER*1
221: // * = 'N': apply H (No transpose)
222: // * = 'C': apply H' (Transpose)
223: // *
224: // * DIRECT (input) CHARACTER*1
225: // * Indicates how H is formed from a product of elementary
226: // * reflectors
227: // * = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
228: // * = 'B': H = H(k) . . . H(2) H(1) (Backward)
229: // *
230: // * STOREV (input) CHARACTER*1
231: // * Indicates how the vectors which define the elementary
232: // * reflectors are stored:
233: // * = 'C': Columnwise (not supported yet)
234: // * = 'R': Rowwise
235: // *
236: // * M (input) INTEGER
237: // * The number of rows of the matrix C.
238: // *
239: // * N (input) INTEGER
240: // * The number of columns of the matrix C.
241: // *
242: // * K (input) INTEGER
243: // * The order of the matrix T (= the number of elementary
244: // * reflectors whose product defines the block reflector).
245: // *
246: // * L (input) INTEGER
247: // * The number of columns of the matrix V containing the
248: // * meaningful part of the Householder reflectors.
249: // * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
250: // *
251: // * V (input) DOUBLE PRECISION array, dimension (LDV,NV).
252: // * If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
253: // *
254: // * LDV (input) INTEGER
255: // * The leading dimension of the array V.
256: // * If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
257: // *
258: // * T (input) DOUBLE PRECISION array, dimension (LDT,K)
259: // * The triangular K-by-K matrix T in the representation of the
260: // * block reflector.
261: // *
262: // * LDT (input) INTEGER
263: // * The leading dimension of the array T. LDT >= K.
264: // *
265: // * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
266: // * On entry, the M-by-N matrix C.
267: // * On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
268: // *
269: // * LDC (input) INTEGER
270: // * The leading dimension of the array C. LDC >= max(1,M).
271: // *
272: // * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
273: // *
274: // * LDWORK (input) INTEGER
275: // * The leading dimension of the array WORK.
276: // * If SIDE = 'L', LDWORK >= max(1,N);
277: // * if SIDE = 'R', LDWORK >= max(1,M).
278: // *
279: // * Further Details
280: // * ===============
281: // *
282: // * Based on contributions by
283: // * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
284: // *
285: // * =====================================================================
286: // *
287: // * .. Parameters ..
288: // * ..
289: // * .. Local Scalars ..
290: // * ..
291: // * .. External Functions ..
292: // * ..
293: // * .. External Subroutines ..
294: // * ..
295: // * .. Executable Statements ..
296: // *
297: // * Quick return if possible
298: // *
299:
300: #endregion
301:
302:
303: #region Body
304:
305: if (M <= 0 || N <= 0) return;
306: // *
307: // * Check for currently supported options
308: // *
309: INFO = 0;
310: if (!this._lsame.Run(DIRECT, "B"))
311: {
312: INFO = - 3;
313: }
314: else
315: {
316: if (!this._lsame.Run(STOREV, "R"))
317: {
318: INFO = - 4;
319: }
320: }
321: if (INFO != 0)
322: {
323: this._xerbla.Run("DLARZB", - INFO);
324: return;
325: }
326: // *
327: if (this._lsame.Run(TRANS, "N"))
328: {
329: FortranLib.Copy(ref TRANST , "T");
330: }
331: else
332: {
333: FortranLib.Copy(ref TRANST , "N");
334: }
335: // *
336: if (this._lsame.Run(SIDE, "L"))
337: {
338: // *
339: // * Form H * C or H' * C
340: // *
341: // * W( 1:n, 1:k ) = C( 1:k, 1:n )'
342: // *
343: for (J = 1; J <= K; J++)
344: {
345: this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
346: }
347: // *
348: // * W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
349: // * C( m-l+1:m, 1:n )' * V( 1:k, 1:l )'
350: // *
351: if (L > 0)
352: {
353: this._dgemm.Run("Transpose", "Transpose", N, K, L, ONE
354: , C, M - L + 1+1 * LDC + o_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
355: , LDWORK);
356: }
357: // *
358: // * W( 1:n, 1:k ) = W( 1:n, 1:k ) * T' or W( 1:m, 1:k ) * T
359: // *
360: this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
361: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
362: // *
363: // * C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )'
364: // *
365: for (J = 1; J <= N; J++)
366: {
367: for (I = 1; I <= K; I++)
368: {
369: C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[J+I * LDWORK + o_work];
370: }
371: }
372: // *
373: // * C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
374: // * V( 1:k, 1:l )' * W( 1:n, 1:k )'
375: // *
376: if (L > 0)
377: {
378: this._dgemm.Run("Transpose", "Transpose", L, N, K, - ONE
379: , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, M - L + 1+1 * LDC + o_c
380: , LDC);
381: }
382: // *
383: }
384: else
385: {
386: if (this._lsame.Run(SIDE, "R"))
387: {
388: // *
389: // * Form C * H or C * H'
390: // *
391: // * W( 1:m, 1:k ) = C( 1:m, 1:k )
392: // *
393: for (J = 1; J <= K; J++)
394: {
395: this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
396: }
397: // *
398: // * W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
399: // * C( 1:m, n-l+1:n ) * V( 1:k, 1:l )'
400: // *
401: if (L > 0)
402: {
403: this._dgemm.Run("No transpose", "Transpose", M, K, L, ONE
404: , C, 1+(N - L + 1) * LDC + o_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
405: , LDWORK);
406: }
407: // *
408: // * W( 1:m, 1:k ) = W( 1:m, 1:k ) * T or W( 1:m, 1:k ) * T'
409: // *
410: this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
411: , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
412: // *
413: // * C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
414: // *
415: for (J = 1; J <= K; J++)
416: {
417: for (I = 1; I <= M; I++)
418: {
419: C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
420: }
421: }
422: // *
423: // * C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
424: // * W( 1:m, 1:k ) * V( 1:k, 1:l )
425: // *
426: if (L > 0)
427: {
428: this._dgemm.Run("No transpose", "No transpose", M, L, K, - ONE
429: , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, 1+(N - L + 1) * LDC + o_c
430: , LDC);
431: }
432: // *
433: }
434: }
435: // *
436: return;
437: // *
438: // * End of DLARZB
439: // *
440:
441: #endregion
442:
443: }
444: }
445: }