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: /// DGBTF2 computes an LU factorization of a real m-by-n band matrix A
27: /// using partial pivoting with row interchanges.
28: ///
29: /// This is the unblocked version of the algorithm, calling Level 2 BLAS.
30: ///
31: ///</summary>
32: public class DGBTF2
33: {
34:
35:
36: #region Dependencies
37:
38: IDAMAX _idamax; DGER _dger; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int I = 0; int J = 0; int JP = 0; int JU = 0; int KM = 0;
46: int KV = 0;
47:
48: #endregion
49:
50: public DGBTF2(IDAMAX idamax, DGER dger, DSCAL dscal, DSWAP dswap, XERBLA xerbla)
51: {
52:
53:
54: #region Set Dependencies
55:
56: this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; this._xerbla = xerbla;
57:
58: #endregion
59:
60: }
61:
62: public DGBTF2()
63: {
64:
65:
66: #region Dependencies (Initialization)
67:
68: IDAMAX idamax = new IDAMAX();
69: XERBLA xerbla = new XERBLA();
70: DSCAL dscal = new DSCAL();
71: DSWAP dswap = new DSWAP();
72: DGER dger = new DGER(xerbla);
73:
74: #endregion
75:
76:
77: #region Set Dependencies
78:
79: this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; this._xerbla = xerbla;
80:
81: #endregion
82:
83: }
84: /// <summary>
85: /// Purpose
86: /// =======
87: ///
88: /// DGBTF2 computes an LU factorization of a real m-by-n band matrix A
89: /// using partial pivoting with row interchanges.
90: ///
91: /// This is the unblocked version of the algorithm, calling Level 2 BLAS.
92: ///
93: ///</summary>
94: /// <param name="M">
95: /// (input) INTEGER
96: /// The number of rows of the matrix A. M .GE. 0.
97: ///</param>
98: /// <param name="N">
99: /// (input) INTEGER
100: /// The number of columns of the matrix A. N .GE. 0.
101: ///</param>
102: /// <param name="KL">
103: /// (input) INTEGER
104: /// The number of subdiagonals within the band of A. KL .GE. 0.
105: ///</param>
106: /// <param name="KU">
107: /// (input) INTEGER
108: /// The number of superdiagonals within the band of A. KU .GE. 0.
109: ///</param>
110: /// <param name="AB">
111: /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
112: /// On entry, the matrix A in band storage, in rows KL+1 to
113: /// 2*KL+KU+1; rows 1 to KL of the array need not be set.
114: /// The j-th column of A is stored in the j-th column of the
115: /// array AB as follows:
116: /// AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku).LE.i.LE.min(m,j+kl)
117: ///
118: /// On exit, details of the factorization: U is stored as an
119: /// upper triangular band matrix with KL+KU superdiagonals in
120: /// rows 1 to KL+KU+1, and the multipliers used during the
121: /// factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
122: /// See below for further details.
123: ///</param>
124: /// <param name="LDAB">
125: /// (input) INTEGER
126: /// The leading dimension of the array AB. LDAB .GE. 2*KL+KU+1.
127: ///</param>
128: /// <param name="IPIV">
129: /// (output) INTEGER array, dimension (min(M,N))
130: /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
131: /// matrix was interchanged with row IPIV(i).
132: ///</param>
133: /// <param name="INFO">
134: /// (output) INTEGER
135: /// = 0: successful exit
136: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
137: /// .GT. 0: if INFO = +i, U(i,i) is exactly zero. The factorization
138: /// has been completed, but the factor U is exactly
139: /// singular, and division by zero will occur if it is used
140: /// to solve a system of equations.
141: ///</param>
142: public void Run(int M, int N, int KL, int KU, ref double[] AB, int offset_ab, int LDAB
143: , ref int[] IPIV, int offset_ipiv, ref int INFO)
144: {
145:
146: #region Array Index Correction
147:
148: int o_ab = -1 - LDAB + offset_ab; int o_ipiv = -1 + offset_ipiv;
149:
150: #endregion
151:
152:
153: #region Prolog
154:
155: // *
156: // * -- LAPACK routine (version 3.1) --
157: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
158: // * November 2006
159: // *
160: // * .. Scalar Arguments ..
161: // * ..
162: // * .. Array Arguments ..
163: // * ..
164: // *
165: // * Purpose
166: // * =======
167: // *
168: // * DGBTF2 computes an LU factorization of a real m-by-n band matrix A
169: // * using partial pivoting with row interchanges.
170: // *
171: // * This is the unblocked version of the algorithm, calling Level 2 BLAS.
172: // *
173: // * Arguments
174: // * =========
175: // *
176: // * M (input) INTEGER
177: // * The number of rows of the matrix A. M >= 0.
178: // *
179: // * N (input) INTEGER
180: // * The number of columns of the matrix A. N >= 0.
181: // *
182: // * KL (input) INTEGER
183: // * The number of subdiagonals within the band of A. KL >= 0.
184: // *
185: // * KU (input) INTEGER
186: // * The number of superdiagonals within the band of A. KU >= 0.
187: // *
188: // * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
189: // * On entry, the matrix A in band storage, in rows KL+1 to
190: // * 2*KL+KU+1; rows 1 to KL of the array need not be set.
191: // * The j-th column of A is stored in the j-th column of the
192: // * array AB as follows:
193: // * AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
194: // *
195: // * On exit, details of the factorization: U is stored as an
196: // * upper triangular band matrix with KL+KU superdiagonals in
197: // * rows 1 to KL+KU+1, and the multipliers used during the
198: // * factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
199: // * See below for further details.
200: // *
201: // * LDAB (input) INTEGER
202: // * The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
203: // *
204: // * IPIV (output) INTEGER array, dimension (min(M,N))
205: // * The pivot indices; for 1 <= i <= min(M,N), row i of the
206: // * matrix was interchanged with row IPIV(i).
207: // *
208: // * INFO (output) INTEGER
209: // * = 0: successful exit
210: // * < 0: if INFO = -i, the i-th argument had an illegal value
211: // * > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
212: // * has been completed, but the factor U is exactly
213: // * singular, and division by zero will occur if it is used
214: // * to solve a system of equations.
215: // *
216: // * Further Details
217: // * ===============
218: // *
219: // * The band storage scheme is illustrated by the following example, when
220: // * M = N = 6, KL = 2, KU = 1:
221: // *
222: // * On entry: On exit:
223: // *
224: // * * * * + + + * * * u14 u25 u36
225: // * * * + + + + * * u13 u24 u35 u46
226: // * * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
227: // * a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
228: // * a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
229: // * a31 a42 a53 a64 * * m31 m42 m53 m64 * *
230: // *
231: // * Array elements marked * are not used by the routine; elements marked
232: // * + need not be set on entry, but are required by the routine to store
233: // * elements of U, because of fill-in resulting from the row
234: // * interchanges.
235: // *
236: // * =====================================================================
237: // *
238: // * .. Parameters ..
239: // * ..
240: // * .. Local Scalars ..
241: // * ..
242: // * .. External Functions ..
243: // * ..
244: // * .. External Subroutines ..
245: // * ..
246: // * .. Intrinsic Functions ..
247: // INTRINSIC MAX, MIN;
248: // * ..
249: // * .. Executable Statements ..
250: // *
251: // * KV is the number of superdiagonals in the factor U, allowing for
252: // * fill-in.
253: // *
254:
255: #endregion
256:
257:
258: #region Body
259:
260: KV = KU + KL;
261: // *
262: // * Test the input parameters.
263: // *
264: INFO = 0;
265: if (M < 0)
266: {
267: INFO = - 1;
268: }
269: else
270: {
271: if (N < 0)
272: {
273: INFO = - 2;
274: }
275: else
276: {
277: if (KL < 0)
278: {
279: INFO = - 3;
280: }
281: else
282: {
283: if (KU < 0)
284: {
285: INFO = - 4;
286: }
287: else
288: {
289: if (LDAB < KL + KV + 1)
290: {
291: INFO = - 6;
292: }
293: }
294: }
295: }
296: }
297: if (INFO != 0)
298: {
299: this._xerbla.Run("DGBTF2", - INFO);
300: return;
301: }
302: // *
303: // * Quick return if possible
304: // *
305: if (M == 0 || N == 0) return;
306: // *
307: // * Gaussian elimination with partial pivoting
308: // *
309: // * Set fill-in elements in columns KU+2 to KV to zero.
310: // *
311: for (J = KU + 2; J <= Math.Min(KV, N); J++)
312: {
313: for (I = KV - J + 2; I <= KL; I++)
314: {
315: AB[I+J * LDAB + o_ab] = ZERO;
316: }
317: }
318: // *
319: // * JU is the index of the last column affected by the current stage
320: // * of the factorization.
321: // *
322: JU = 1;
323: // *
324: for (J = 1; J <= Math.Min(M, N); J++)
325: {
326: // *
327: // * Set fill-in elements in column J+KV to zero.
328: // *
329: if (J + KV <= N)
330: {
331: for (I = 1; I <= KL; I++)
332: {
333: AB[I+(J + KV) * LDAB + o_ab] = ZERO;
334: }
335: }
336: // *
337: // * Find pivot and test for singularity. KM is the number of
338: // * subdiagonal elements in the current column.
339: // *
340: KM = Math.Min(KL, M - J);
341: JP = this._idamax.Run(KM + 1, AB, KV + 1+J * LDAB + o_ab, 1);
342: IPIV[J + o_ipiv] = JP + J - 1;
343: if (AB[KV + JP+J * LDAB + o_ab] != ZERO)
344: {
345: JU = Math.Max(JU, Math.Min(J + KU + JP - 1, N));
346: // *
347: // * Apply interchange to columns J to JU.
348: // *
349: if (JP != 1) this._dswap.Run(JU - J + 1, ref AB, KV + JP+J * LDAB + o_ab, LDAB - 1, ref AB, KV + 1+J * LDAB + o_ab, LDAB - 1);
350: // *
351: if (KM > 0)
352: {
353: // *
354: // * Compute multipliers.
355: // *
356: this._dscal.Run(KM, ONE / AB[KV + 1+J * LDAB + o_ab], ref AB, KV + 2+J * LDAB + o_ab, 1);
357: // *
358: // * Update trailing submatrix within the band.
359: // *
360: if (JU > J)
361: {
362: this._dger.Run(KM, JU - J, - ONE, AB, KV + 2+J * LDAB + o_ab, 1, AB, KV+(J + 1) * LDAB + o_ab
363: , LDAB - 1, ref AB, KV + 1+(J + 1) * LDAB + o_ab, LDAB - 1);
364: }
365: }
366: }
367: else
368: {
369: // *
370: // * If pivot is zero, set INFO to the index of the pivot
371: // * unless a zero pivot has already been found.
372: // *
373: if (INFO == 0) INFO = J;
374: }
375: }
376: return;
377: // *
378: // * End of DGBTF2
379: // *
380:
381: #endregion
382:
383: }
384: }
385: }