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