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: /// DLATRZ factors the M-by-(M+L) real upper trapezoidal matrix
27: /// [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R 0 ) * Z, by means
28: /// of orthogonal transformations. Z is an (M+L)-by-(M+L) orthogonal
29: /// matrix and, R and A1 are M-by-M upper triangular matrices.
30: ///
31: ///</summary>
32: public class DLATRZ
33: {
34:
35:
36: #region Dependencies
37:
38: DLARFG _dlarfg; DLARZ _dlarz;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ZERO = 0.0E+0; int I = 0;
46:
47: #endregion
48:
49: public DLATRZ(DLARFG dlarfg, DLARZ dlarz)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dlarfg = dlarfg; this._dlarz = dlarz;
56:
57: #endregion
58:
59: }
60:
61: public DLATRZ()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: LSAME lsame = new LSAME();
68: DLAMC3 dlamc3 = new DLAMC3();
69: DLAPY2 dlapy2 = new DLAPY2();
70: DNRM2 dnrm2 = new DNRM2();
71: DSCAL dscal = new DSCAL();
72: DAXPY daxpy = new DAXPY();
73: DCOPY dcopy = new DCOPY();
74: XERBLA xerbla = new XERBLA();
75: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
76: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
77: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
78: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
79: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
80: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
81: DGEMV dgemv = new DGEMV(lsame, xerbla);
82: DGER dger = new DGER(xerbla);
83: DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);
84:
85: #endregion
86:
87:
88: #region Set Dependencies
89:
90: this._dlarfg = dlarfg; this._dlarz = dlarz;
91:
92: #endregion
93:
94: }
95: /// <summary>
96: /// Purpose
97: /// =======
98: ///
99: /// DLATRZ factors the M-by-(M+L) real upper trapezoidal matrix
100: /// [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R 0 ) * Z, by means
101: /// of orthogonal transformations. Z is an (M+L)-by-(M+L) orthogonal
102: /// matrix and, R and A1 are M-by-M upper triangular matrices.
103: ///
104: ///</summary>
105: /// <param name="M">
106: /// (input) INTEGER
107: /// The number of rows of the matrix A. M .GE. 0.
108: ///</param>
109: /// <param name="N">
110: /// (input) INTEGER
111: /// The number of columns of the matrix A. N .GE. 0.
112: ///</param>
113: /// <param name="L">
114: /// (input) INTEGER
115: /// The number of columns of the matrix A containing the
116: /// meaningful part of the Householder vectors. N-M .GE. L .GE. 0.
117: ///</param>
118: /// <param name="A">
119: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
120: /// On entry, the leading M-by-N upper trapezoidal part of the
121: /// array A must contain the matrix to be factorized.
122: /// On exit, the leading M-by-M upper triangular part of A
123: /// contains the upper triangular matrix R, and elements N-L+1 to
124: /// N of the first M rows of A, with the array TAU, represent the
125: /// orthogonal matrix Z as a product of M elementary reflectors.
126: ///</param>
127: /// <param name="LDA">
128: /// (input) INTEGER
129: /// The leading dimension of the array A. LDA .GE. max(1,M).
130: ///</param>
131: /// <param name="TAU">
132: /// (output) DOUBLE PRECISION array, dimension (M)
133: /// The scalar factors of the elementary reflectors.
134: ///</param>
135: /// <param name="WORK">
136: /// (workspace) DOUBLE PRECISION array, dimension (M)
137: ///</param>
138: public void Run(int M, int N, int L, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
139: , ref double[] WORK, int offset_work)
140: {
141:
142: #region Array Index Correction
143:
144: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
145:
146: #endregion
147:
148:
149: #region Prolog
150:
151: // *
152: // * -- LAPACK routine (version 3.1) --
153: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
154: // * November 2006
155: // *
156: // * .. Scalar Arguments ..
157: // * ..
158: // * .. Array Arguments ..
159: // * ..
160: // *
161: // * Purpose
162: // * =======
163: // *
164: // * DLATRZ factors the M-by-(M+L) real upper trapezoidal matrix
165: // * [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R 0 ) * Z, by means
166: // * of orthogonal transformations. Z is an (M+L)-by-(M+L) orthogonal
167: // * matrix and, R and A1 are M-by-M upper triangular matrices.
168: // *
169: // * Arguments
170: // * =========
171: // *
172: // * M (input) INTEGER
173: // * The number of rows of the matrix A. M >= 0.
174: // *
175: // * N (input) INTEGER
176: // * The number of columns of the matrix A. N >= 0.
177: // *
178: // * L (input) INTEGER
179: // * The number of columns of the matrix A containing the
180: // * meaningful part of the Householder vectors. N-M >= L >= 0.
181: // *
182: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
183: // * On entry, the leading M-by-N upper trapezoidal part of the
184: // * array A must contain the matrix to be factorized.
185: // * On exit, the leading M-by-M upper triangular part of A
186: // * contains the upper triangular matrix R, and elements N-L+1 to
187: // * N of the first M rows of A, with the array TAU, represent the
188: // * orthogonal matrix Z as a product of M elementary reflectors.
189: // *
190: // * LDA (input) INTEGER
191: // * The leading dimension of the array A. LDA >= max(1,M).
192: // *
193: // * TAU (output) DOUBLE PRECISION array, dimension (M)
194: // * The scalar factors of the elementary reflectors.
195: // *
196: // * WORK (workspace) DOUBLE PRECISION array, dimension (M)
197: // *
198: // * Further Details
199: // * ===============
200: // *
201: // * Based on contributions by
202: // * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
203: // *
204: // * The factorization is obtained by Householder's method. The kth
205: // * transformation matrix, Z( k ), which is used to introduce zeros into
206: // * the ( m - k + 1 )th row of A, is given in the form
207: // *
208: // * Z( k ) = ( I 0 ),
209: // * ( 0 T( k ) )
210: // *
211: // * where
212: // *
213: // * T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ),
214: // * ( 0 )
215: // * ( z( k ) )
216: // *
217: // * tau is a scalar and z( k ) is an l element vector. tau and z( k )
218: // * are chosen to annihilate the elements of the kth row of A2.
219: // *
220: // * The scalar tau is returned in the kth element of TAU and the vector
221: // * u( k ) in the kth row of A2, such that the elements of z( k ) are
222: // * in a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in
223: // * the upper triangular part of A1.
224: // *
225: // * Z is given by
226: // *
227: // * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ).
228: // *
229: // * =====================================================================
230: // *
231: // * .. Parameters ..
232: // * ..
233: // * .. Local Scalars ..
234: // * ..
235: // * .. External Subroutines ..
236: // * ..
237: // * .. Executable Statements ..
238: // *
239: // * Test the input arguments
240: // *
241: // * Quick return if possible
242: // *
243:
244: #endregion
245:
246:
247: #region Body
248:
249: if (M == 0)
250: {
251: return;
252: }
253: else
254: {
255: if (M == N)
256: {
257: for (I = 1; I <= N; I++)
258: {
259: TAU[I + o_tau] = ZERO;
260: }
261: return;
262: }
263: }
264: // *
265: for (I = M; I >= 1; I += - 1)
266: {
267: // *
268: // * Generate elementary reflector H(i) to annihilate
269: // * [ A(i,i) A(i,n-l+1:n) ]
270: // *
271: this._dlarfg.Run(L + 1, ref A[I+I * LDA + o_a], ref A, I+(N - L + 1) * LDA + o_a, LDA, ref TAU[I + o_tau]);
272: // *
273: // * Apply H(i) to A(1:i-1,i:n) from the right
274: // *
275: this._dlarz.Run("Right", I - 1, N - I + 1, L, A, I+(N - L + 1) * LDA + o_a, LDA
276: , TAU[I + o_tau], ref A, 1+I * LDA + o_a, LDA, ref WORK, offset_work);
277: // *
278: }
279: // *
280: return;
281: // *
282: // * End of DLATRZ
283: // *
284:
285: #endregion
286:
287: }
288: }
289: }