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: /// DGTSV solves the equation
27: ///
28: /// A*X = B,
29: ///
30: /// where A is an n by n tridiagonal matrix, by Gaussian elimination with
31: /// partial pivoting.
32: ///
33: /// Note that the equation A'*X = B may be solved by interchanging the
34: /// order of the arguments DU and DL.
35: ///
36: ///</summary>
37: public class DGTSV
38: {
39:
40:
41: #region Dependencies
42:
43: XERBLA _xerbla;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const double ZERO = 0.0E+0; int I = 0; int J = 0; double FACT = 0; double TEMP = 0;
51:
52: #endregion
53:
54: public DGTSV(XERBLA xerbla)
55: {
56:
57:
58: #region Set Dependencies
59:
60: this._xerbla = xerbla;
61:
62: #endregion
63:
64: }
65:
66: public DGTSV()
67: {
68:
69:
70: #region Dependencies (Initialization)
71:
72: XERBLA xerbla = new XERBLA();
73:
74: #endregion
75:
76:
77: #region Set Dependencies
78:
79: this._xerbla = xerbla;
80:
81: #endregion
82:
83: }
84: /// <summary>
85: /// Purpose
86: /// =======
87: ///
88: /// DGTSV solves the equation
89: ///
90: /// A*X = B,
91: ///
92: /// where A is an n by n tridiagonal matrix, by Gaussian elimination with
93: /// partial pivoting.
94: ///
95: /// Note that the equation A'*X = B may be solved by interchanging the
96: /// order of the arguments DU and DL.
97: ///
98: ///</summary>
99: /// <param name="N">
100: /// (input) INTEGER
101: /// The order of the matrix A. N .GE. 0.
102: ///</param>
103: /// <param name="NRHS">
104: /// (input) INTEGER
105: /// The number of right hand sides, i.e., the number of columns
106: /// of the matrix B. NRHS .GE. 0.
107: ///</param>
108: /// <param name="DL">
109: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
110: /// On entry, DL must contain the (n-1) sub-diagonal elements of
111: /// A.
112: ///
113: /// On exit, DL is overwritten by the (n-2) elements of the
114: /// second super-diagonal of the upper triangular matrix U from
115: /// the LU factorization of A, in DL(1), ..., DL(n-2).
116: ///</param>
117: /// <param name="D">
118: /// (input/output) DOUBLE PRECISION array, dimension (N)
119: /// On entry, D must contain the diagonal elements of A.
120: ///
121: /// On exit, D is overwritten by the n diagonal elements of U.
122: ///</param>
123: /// <param name="DU">
124: /// (input/output) DOUBLE PRECISION array, dimension (N-1)
125: /// On entry, DU must contain the (n-1) super-diagonal elements
126: /// of A.
127: ///
128: /// On exit, DU is overwritten by the (n-1) elements of the first
129: /// super-diagonal of U.
130: ///</param>
131: /// <param name="B">
132: /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
133: /// On entry, the N by NRHS matrix of right hand side matrix B.
134: /// On exit, if INFO = 0, the N by NRHS solution matrix X.
135: ///</param>
136: /// <param name="LDB">
137: /// (input) INTEGER
138: /// The leading dimension of the array B. LDB .GE. max(1,N).
139: ///</param>
140: /// <param name="INFO">
141: /// (output) INTEGER
142: /// = 0: successful exit
143: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
144: /// .GT. 0: if INFO = i, U(i,i) is exactly zero, and the solution
145: /// has not been computed. The factorization has not been
146: /// completed unless i = N.
147: ///</param>
148: public void Run(int N, int NRHS, ref double[] DL, int offset_dl, ref double[] D, int offset_d, ref double[] DU, int offset_du, ref double[] B, int offset_b
149: , int LDB, ref int INFO)
150: {
151:
152: #region Array Index Correction
153:
154: int o_dl = -1 + offset_dl; int o_d = -1 + offset_d; int o_du = -1 + offset_du; int o_b = -1 - LDB + offset_b;
155:
156: #endregion
157:
158:
159: #region Prolog
160:
161: // *
162: // * -- LAPACK routine (version 3.1) --
163: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
164: // * November 2006
165: // *
166: // * .. Scalar Arguments ..
167: // * ..
168: // * .. Array Arguments ..
169: // * ..
170: // *
171: // * Purpose
172: // * =======
173: // *
174: // * DGTSV solves the equation
175: // *
176: // * A*X = B,
177: // *
178: // * where A is an n by n tridiagonal matrix, by Gaussian elimination with
179: // * partial pivoting.
180: // *
181: // * Note that the equation A'*X = B may be solved by interchanging the
182: // * order of the arguments DU and DL.
183: // *
184: // * Arguments
185: // * =========
186: // *
187: // * N (input) INTEGER
188: // * The order of the matrix A. N >= 0.
189: // *
190: // * NRHS (input) INTEGER
191: // * The number of right hand sides, i.e., the number of columns
192: // * of the matrix B. NRHS >= 0.
193: // *
194: // * DL (input/output) DOUBLE PRECISION array, dimension (N-1)
195: // * On entry, DL must contain the (n-1) sub-diagonal elements of
196: // * A.
197: // *
198: // * On exit, DL is overwritten by the (n-2) elements of the
199: // * second super-diagonal of the upper triangular matrix U from
200: // * the LU factorization of A, in DL(1), ..., DL(n-2).
201: // *
202: // * D (input/output) DOUBLE PRECISION array, dimension (N)
203: // * On entry, D must contain the diagonal elements of A.
204: // *
205: // * On exit, D is overwritten by the n diagonal elements of U.
206: // *
207: // * DU (input/output) DOUBLE PRECISION array, dimension (N-1)
208: // * On entry, DU must contain the (n-1) super-diagonal elements
209: // * of A.
210: // *
211: // * On exit, DU is overwritten by the (n-1) elements of the first
212: // * super-diagonal of U.
213: // *
214: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
215: // * On entry, the N by NRHS matrix of right hand side matrix B.
216: // * On exit, if INFO = 0, the N by NRHS solution matrix X.
217: // *
218: // * LDB (input) INTEGER
219: // * The leading dimension of the array B. LDB >= max(1,N).
220: // *
221: // * INFO (output) INTEGER
222: // * = 0: successful exit
223: // * < 0: if INFO = -i, the i-th argument had an illegal value
224: // * > 0: if INFO = i, U(i,i) is exactly zero, and the solution
225: // * has not been computed. The factorization has not been
226: // * completed unless i = N.
227: // *
228: // * =====================================================================
229: // *
230: // * .. Parameters ..
231: // * ..
232: // * .. Local Scalars ..
233: // * ..
234: // * .. Intrinsic Functions ..
235: // INTRINSIC ABS, MAX;
236: // * ..
237: // * .. External Subroutines ..
238: // * ..
239: // * .. Executable Statements ..
240: // *
241:
242: #endregion
243:
244:
245: #region Body
246:
247: INFO = 0;
248: if (N < 0)
249: {
250: INFO = - 1;
251: }
252: else
253: {
254: if (NRHS < 0)
255: {
256: INFO = - 2;
257: }
258: else
259: {
260: if (LDB < Math.Max(1, N))
261: {
262: INFO = - 7;
263: }
264: }
265: }
266: if (INFO != 0)
267: {
268: this._xerbla.Run("DGTSV ", - INFO);
269: return;
270: }
271: // *
272: if (N == 0) return;
273: // *
274: if (NRHS == 1)
275: {
276: for (I = 1; I <= N - 2; I++)
277: {
278: if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))
279: {
280: // *
281: // * No row interchange required
282: // *
283: if (D[I + o_d] != ZERO)
284: {
285: FACT = DL[I + o_dl] / D[I + o_d];
286: D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];
287: B[I + 1+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b] - FACT * B[I+1 * LDB + o_b];
288: }
289: else
290: {
291: INFO = I;
292: return;
293: }
294: DL[I + o_dl] = ZERO;
295: }
296: else
297: {
298: // *
299: // * Interchange rows I and I+1
300: // *
301: FACT = D[I + o_d] / DL[I + o_dl];
302: D[I + o_d] = DL[I + o_dl];
303: TEMP = D[I + 1 + o_d];
304: D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;
305: DL[I + o_dl] = DU[I + 1 + o_du];
306: DU[I + 1 + o_du] = - FACT * DL[I + o_dl];
307: DU[I + o_du] = TEMP;
308: TEMP = B[I+1 * LDB + o_b];
309: B[I+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b];
310: B[I + 1+1 * LDB + o_b] = TEMP - FACT * B[I + 1+1 * LDB + o_b];
311: }
312: }
313: if (N > 1)
314: {
315: I = N - 1;
316: if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))
317: {
318: if (D[I + o_d] != ZERO)
319: {
320: FACT = DL[I + o_dl] / D[I + o_d];
321: D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];
322: B[I + 1+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b] - FACT * B[I+1 * LDB + o_b];
323: }
324: else
325: {
326: INFO = I;
327: return;
328: }
329: }
330: else
331: {
332: FACT = D[I + o_d] / DL[I + o_dl];
333: D[I + o_d] = DL[I + o_dl];
334: TEMP = D[I + 1 + o_d];
335: D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;
336: DU[I + o_du] = TEMP;
337: TEMP = B[I+1 * LDB + o_b];
338: B[I+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b];
339: B[I + 1+1 * LDB + o_b] = TEMP - FACT * B[I + 1+1 * LDB + o_b];
340: }
341: }
342: if (D[N + o_d] == ZERO)
343: {
344: INFO = N;
345: return;
346: }
347: }
348: else
349: {
350: for (I = 1; I <= N - 2; I++)
351: {
352: if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))
353: {
354: // *
355: // * No row interchange required
356: // *
357: if (D[I + o_d] != ZERO)
358: {
359: FACT = DL[I + o_dl] / D[I + o_d];
360: D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];
361: for (J = 1; J <= NRHS; J++)
362: {
363: B[I + 1+J * LDB + o_b] = B[I + 1+J * LDB + o_b] - FACT * B[I+J * LDB + o_b];
364: }
365: }
366: else
367: {
368: INFO = I;
369: return;
370: }
371: DL[I + o_dl] = ZERO;
372: }
373: else
374: {
375: // *
376: // * Interchange rows I and I+1
377: // *
378: FACT = D[I + o_d] / DL[I + o_dl];
379: D[I + o_d] = DL[I + o_dl];
380: TEMP = D[I + 1 + o_d];
381: D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;
382: DL[I + o_dl] = DU[I + 1 + o_du];
383: DU[I + 1 + o_du] = - FACT * DL[I + o_dl];
384: DU[I + o_du] = TEMP;
385: for (J = 1; J <= NRHS; J++)
386: {
387: TEMP = B[I+J * LDB + o_b];
388: B[I+J * LDB + o_b] = B[I + 1+J * LDB + o_b];
389: B[I + 1+J * LDB + o_b] = TEMP - FACT * B[I + 1+J * LDB + o_b];
390: }
391: }
392: }
393: if (N > 1)
394: {
395: I = N - 1;
396: if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))
397: {
398: if (D[I + o_d] != ZERO)
399: {
400: FACT = DL[I + o_dl] / D[I + o_d];
401: D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];
402: for (J = 1; J <= NRHS; J++)
403: {
404: B[I + 1+J * LDB + o_b] = B[I + 1+J * LDB + o_b] - FACT * B[I+J * LDB + o_b];
405: }
406: }
407: else
408: {
409: INFO = I;
410: return;
411: }
412: }
413: else
414: {
415: FACT = D[I + o_d] / DL[I + o_dl];
416: D[I + o_d] = DL[I + o_dl];
417: TEMP = D[I + 1 + o_d];
418: D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;
419: DU[I + o_du] = TEMP;
420: for (J = 1; J <= NRHS; J++)
421: {
422: TEMP = B[I+J * LDB + o_b];
423: B[I+J * LDB + o_b] = B[I + 1+J * LDB + o_b];
424: B[I + 1+J * LDB + o_b] = TEMP - FACT * B[I + 1+J * LDB + o_b];
425: }
426: }
427: }
428: if (D[N + o_d] == ZERO)
429: {
430: INFO = N;
431: return;
432: }
433: }
434: // *
435: // * Back solve with the matrix U from the factorization.
436: // *
437: if (NRHS <= 2)
438: {
439: J = 1;
440: LABEL70:;
441: B[N+J * LDB + o_b] = B[N+J * LDB + o_b] / D[N + o_d];
442: if (N > 1) B[N - 1+J * LDB + o_b] = (B[N - 1+J * LDB + o_b] - DU[N - 1 + o_du] * B[N+J * LDB + o_b]) / D[N - 1 + o_d];
443: for (I = N - 2; I >= 1; I += - 1)
444: {
445: B[I+J * LDB + o_b] = (B[I+J * LDB + o_b] - DU[I + o_du] * B[I + 1+J * LDB + o_b] - DL[I + o_dl] * B[I + 2+J * LDB + o_b]) / D[I + o_d];
446: }
447: if (J < NRHS)
448: {
449: J = J + 1;
450: goto LABEL70;
451: }
452: }
453: else
454: {
455: for (J = 1; J <= NRHS; J++)
456: {
457: B[N+J * LDB + o_b] = B[N+J * LDB + o_b] / D[N + o_d];
458: if (N > 1) B[N - 1+J * LDB + o_b] = (B[N - 1+J * LDB + o_b] - DU[N - 1 + o_du] * B[N+J * LDB + o_b]) / D[N - 1 + o_d];
459: for (I = N - 2; I >= 1; I += - 1)
460: {
461: B[I+J * LDB + o_b] = (B[I+J * LDB + o_b] - DU[I + o_du] * B[I + 1+J * LDB + o_b] - DL[I + o_dl] * B[I + 2+J * LDB + o_b]) / D[I + o_d];
462: }
463: }
464: }
465: // *
466: return;
467: // *
468: // * End of DGTSV
469: // *
470:
471: #endregion
472:
473: }
474: }
475: }