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: /// DGETF2 computes an LU factorization of a general m-by-n matrix A
27: /// using partial pivoting with row interchanges.
28: ///
29: /// The factorization has the form
30: /// A = P * L * U
31: /// where P is a permutation matrix, L is lower triangular with unit
32: /// diagonal elements (lower trapezoidal if m .GT. n), and U is upper
33: /// triangular (upper trapezoidal if m .LT. n).
34: ///
35: /// This is the right-looking Level 2 BLAS version of the algorithm.
36: ///
37: ///</summary>
38: public class DGETF2
39: {
40:
41:
42: #region Dependencies
43:
44: DLAMCH _dlamch; IDAMAX _idamax; DGER _dger; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla;
45:
46: #endregion
47:
48:
49: #region Fields
50:
51: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; double SFMIN = 0; int I = 0; int J = 0; int JP = 0;
52:
53: #endregion
54:
55: public DGETF2(DLAMCH dlamch, IDAMAX idamax, DGER dger, DSCAL dscal, DSWAP dswap, XERBLA xerbla)
56: {
57:
58:
59: #region Set Dependencies
60:
61: this._dlamch = dlamch; this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap;
62: this._xerbla = xerbla;
63:
64: #endregion
65:
66: }
67:
68: public DGETF2()
69: {
70:
71:
72: #region Dependencies (Initialization)
73:
74: LSAME lsame = new LSAME();
75: DLAMC3 dlamc3 = new DLAMC3();
76: IDAMAX idamax = new IDAMAX();
77: XERBLA xerbla = new XERBLA();
78: DSCAL dscal = new DSCAL();
79: DSWAP dswap = new DSWAP();
80: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
81: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
82: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
83: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
84: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
85: DGER dger = new DGER(xerbla);
86:
87: #endregion
88:
89:
90: #region Set Dependencies
91:
92: this._dlamch = dlamch; this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap;
93: this._xerbla = xerbla;
94:
95: #endregion
96:
97: }
98: /// <summary>
99: /// Purpose
100: /// =======
101: ///
102: /// DGETF2 computes an LU factorization of a general m-by-n matrix A
103: /// using partial pivoting with row interchanges.
104: ///
105: /// The factorization has the form
106: /// A = P * L * U
107: /// where P is a permutation matrix, L is lower triangular with unit
108: /// diagonal elements (lower trapezoidal if m .GT. n), and U is upper
109: /// triangular (upper trapezoidal if m .LT. n).
110: ///
111: /// This is the right-looking Level 2 BLAS version of the algorithm.
112: ///
113: ///</summary>
114: /// <param name="M">
115: /// (input) INTEGER
116: /// The number of rows of the matrix A. M .GE. 0.
117: ///</param>
118: /// <param name="N">
119: /// (input) INTEGER
120: /// The number of columns of the matrix A. N .GE. 0.
121: ///</param>
122: /// <param name="A">
123: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
124: /// On entry, the m by n matrix to be factored.
125: /// On exit, the factors L and U from the factorization
126: /// A = P*L*U; the unit diagonal elements of L are not stored.
127: ///</param>
128: /// <param name="LDA">
129: /// (input) INTEGER
130: /// The leading dimension of the array A. LDA .GE. max(1,M).
131: ///</param>
132: /// <param name="IPIV">
133: /// (output) INTEGER array, dimension (min(M,N))
134: /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
135: /// matrix was interchanged with row IPIV(i).
136: ///</param>
137: /// <param name="INFO">
138: /// (output) INTEGER
139: /// = 0: successful exit
140: /// .LT. 0: if INFO = -k, the k-th argument had an illegal value
141: /// .GT. 0: if INFO = k, U(k,k) is exactly zero. The factorization
142: /// has been completed, but the factor U is exactly
143: /// singular, and division by zero will occur if it is used
144: /// to solve a system of equations.
145: ///</param>
146: public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref int[] IPIV, int offset_ipiv, ref int INFO)
147: {
148:
149: #region Array Index Correction
150:
151: int o_a = -1 - LDA + offset_a; int o_ipiv = -1 + offset_ipiv;
152:
153: #endregion
154:
155:
156: #region Prolog
157:
158: // *
159: // * -- LAPACK routine (version 3.1) --
160: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
161: // * November 2006
162: // *
163: // * .. Scalar Arguments ..
164: // * ..
165: // * .. Array Arguments ..
166: // * ..
167: // *
168: // * Purpose
169: // * =======
170: // *
171: // * DGETF2 computes an LU factorization of a general m-by-n matrix A
172: // * using partial pivoting with row interchanges.
173: // *
174: // * The factorization has the form
175: // * A = P * L * U
176: // * where P is a permutation matrix, L is lower triangular with unit
177: // * diagonal elements (lower trapezoidal if m > n), and U is upper
178: // * triangular (upper trapezoidal if m < n).
179: // *
180: // * This is the right-looking Level 2 BLAS version of the algorithm.
181: // *
182: // * Arguments
183: // * =========
184: // *
185: // * M (input) INTEGER
186: // * The number of rows of the matrix A. M >= 0.
187: // *
188: // * N (input) INTEGER
189: // * The number of columns of the matrix A. N >= 0.
190: // *
191: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
192: // * On entry, the m by n matrix to be factored.
193: // * On exit, the factors L and U from the factorization
194: // * A = P*L*U; the unit diagonal elements of L are not stored.
195: // *
196: // * LDA (input) INTEGER
197: // * The leading dimension of the array A. LDA >= max(1,M).
198: // *
199: // * IPIV (output) INTEGER array, dimension (min(M,N))
200: // * The pivot indices; for 1 <= i <= min(M,N), row i of the
201: // * matrix was interchanged with row IPIV(i).
202: // *
203: // * INFO (output) INTEGER
204: // * = 0: successful exit
205: // * < 0: if INFO = -k, the k-th argument had an illegal value
206: // * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
207: // * has been completed, but the factor U is exactly
208: // * singular, and division by zero will occur if it is used
209: // * to solve a system of equations.
210: // *
211: // * =====================================================================
212: // *
213: // * .. Parameters ..
214: // * ..
215: // * .. Local Scalars ..
216: // * ..
217: // * .. External Functions ..
218: // * ..
219: // * .. External Subroutines ..
220: // * ..
221: // * .. Intrinsic Functions ..
222: // INTRINSIC MAX, MIN;
223: // * ..
224: // * .. Executable Statements ..
225: // *
226: // * Test the input parameters.
227: // *
228:
229: #endregion
230:
231:
232: #region Body
233:
234: INFO = 0;
235: if (M < 0)
236: {
237: INFO = - 1;
238: }
239: else
240: {
241: if (N < 0)
242: {
243: INFO = - 2;
244: }
245: else
246: {
247: if (LDA < Math.Max(1, M))
248: {
249: INFO = - 4;
250: }
251: }
252: }
253: if (INFO != 0)
254: {
255: this._xerbla.Run("DGETF2", - INFO);
256: return;
257: }
258: // *
259: // * Quick return if possible
260: // *
261: if (M == 0 || N == 0) return;
262: // *
263: // * Compute machine safe minimum
264: // *
265: SFMIN = this._dlamch.Run("S");
266: // *
267: for (J = 1; J <= Math.Min(M, N); J++)
268: {
269: // *
270: // * Find pivot and test for singularity.
271: // *
272: JP = J - 1 + this._idamax.Run(M - J + 1, A, J+J * LDA + o_a, 1);
273: IPIV[J + o_ipiv] = JP;
274: if (A[JP+J * LDA + o_a] != ZERO)
275: {
276: // *
277: // * Apply the interchange to columns 1:N.
278: // *
279: if (JP != J) this._dswap.Run(N, ref A, J+1 * LDA + o_a, LDA, ref A, JP+1 * LDA + o_a, LDA);
280: // *
281: // * Compute elements J+1:M of J-th column.
282: // *
283: if (J < M)
284: {
285: if (Math.Abs(A[J+J * LDA + o_a]) >= SFMIN)
286: {
287: this._dscal.Run(M - J, ONE / A[J+J * LDA + o_a], ref A, J + 1+J * LDA + o_a, 1);
288: }
289: else
290: {
291: for (I = 1; I <= M - J; I++)
292: {
293: A[J + I+J * LDA + o_a] = A[J + I+J * LDA + o_a] / A[J+J * LDA + o_a];
294: }
295: }
296: }
297: // *
298: }
299: else
300: {
301: if (INFO == 0)
302: {
303: // *
304: INFO = J;
305: }
306: }
307: // *
308: if (J < Math.Min(M, N))
309: {
310: // *
311: // * Update trailing submatrix.
312: // *
313: this._dger.Run(M - J, N - J, - ONE, A, J + 1+J * LDA + o_a, 1, A, J+(J + 1) * LDA + o_a
314: , LDA, ref A, J + 1+(J + 1) * LDA + o_a, LDA);
315: }
316: }
317: return;
318: // *
319: // * End of DGETF2
320: // *
321:
322: #endregion
323:
324: }
325: }
326: }