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: /// Purpose
21: /// =======
22: ///
23: /// DTBSV solves one of the systems of equations
24: ///
25: /// A*x = b, or A'*x = b,
26: ///
27: /// where b and x are n element vectors and A is an n by n unit, or
28: /// non-unit, upper or lower triangular band matrix, with ( k + 1 )
29: /// diagonals.
30: ///
31: /// No test for singularity or near-singularity is included in this
32: /// routine. Such tests must be performed before calling this routine.
33: ///
34: ///</summary>
35: public class DTBSV
36: {
37:
38:
39: #region Dependencies
40:
41: LSAME _lsame; XERBLA _xerbla;
42:
43: #endregion
44:
45:
46: #region Fields
47:
48: const double ZERO = 0.0E+0; double TEMP = 0; int I = 0; int INFO = 0; int IX = 0; int J = 0; int JX = 0; int KPLUS1 = 0;
49: int KX = 0;int L = 0; bool NOUNIT = false;
50:
51: #endregion
52:
53: public DTBSV(LSAME lsame, XERBLA xerbla)
54: {
55:
56:
57: #region Set Dependencies
58:
59: this._lsame = lsame; this._xerbla = xerbla;
60:
61: #endregion
62:
63: }
64:
65: public DTBSV()
66: {
67:
68:
69: #region Dependencies (Initialization)
70:
71: LSAME lsame = new LSAME();
72: XERBLA xerbla = new XERBLA();
73:
74: #endregion
75:
76:
77: #region Set Dependencies
78:
79: this._lsame = lsame; this._xerbla = xerbla;
80:
81: #endregion
82:
83: }
84: /// <summary>
85: /// Purpose
86: /// =======
87: ///
88: /// DTBSV solves one of the systems of equations
89: ///
90: /// A*x = b, or A'*x = b,
91: ///
92: /// where b and x are n element vectors and A is an n by n unit, or
93: /// non-unit, upper or lower triangular band matrix, with ( k + 1 )
94: /// diagonals.
95: ///
96: /// No test for singularity or near-singularity is included in this
97: /// routine. Such tests must be performed before calling this routine.
98: ///
99: ///</summary>
100: /// <param name="UPLO">
101: /// - CHARACTER*1.
102: /// On entry, UPLO specifies whether the matrix is an upper or
103: /// lower triangular matrix as follows:
104: ///
105: /// UPLO = 'U' or 'u' A is an upper triangular matrix.
106: ///
107: /// UPLO = 'L' or 'l' A is a lower triangular matrix.
108: ///
109: /// Unchanged on exit.
110: ///</param>
111: /// <param name="TRANS">
112: /// - CHARACTER*1.
113: /// On entry, TRANS specifies the equations to be solved as
114: /// follows:
115: ///
116: /// TRANS = 'N' or 'n' A*x = b.
117: ///
118: /// TRANS = 'T' or 't' A'*x = b.
119: ///
120: /// TRANS = 'C' or 'c' A'*x = b.
121: ///
122: /// Unchanged on exit.
123: ///</param>
124: /// <param name="DIAG">
125: /// - CHARACTER*1.
126: /// On entry, DIAG specifies whether or not A is unit
127: /// triangular as follows:
128: ///
129: /// DIAG = 'U' or 'u' A is assumed to be unit triangular.
130: ///
131: /// DIAG = 'N' or 'n' A is not assumed to be unit
132: /// triangular.
133: ///
134: /// Unchanged on exit.
135: ///</param>
136: /// <param name="N">
137: /// - INTEGER.
138: /// On entry, N specifies the order of the matrix A.
139: /// N must be at least zero.
140: /// Unchanged on exit.
141: ///</param>
142: /// <param name="K">
143: /// - INTEGER.
144: /// On entry with UPLO = 'U' or 'u', K specifies the number of
145: /// super-diagonals of the matrix A.
146: /// On entry with UPLO = 'L' or 'l', K specifies the number of
147: /// sub-diagonals of the matrix A.
148: /// K must satisfy 0 .le. K.
149: /// Unchanged on exit.
150: ///</param>
151: /// <param name="A">
152: /// - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
153: /// Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
154: /// by n part of the array A must contain the upper triangular
155: /// band part of the matrix of coefficients, supplied column by
156: /// column, with the leading diagonal of the matrix in row
157: /// ( k + 1 ) of the array, the first super-diagonal starting at
158: /// position 2 in row k, and so on. The top left k by k triangle
159: /// of the array A is not referenced.
160: /// The following program segment will transfer an upper
161: /// triangular band matrix from conventional full matrix storage
162: /// to band storage:
163: ///
164: /// DO 20, J = 1, N
165: /// M = K + 1 - J
166: /// DO 10, I = MAX( 1, J - K ), J
167: /// A( M + I, J ) = matrix( I, J )
168: /// 10 CONTINUE
169: /// 20 CONTINUE
170: ///
171: /// Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
172: /// by n part of the array A must contain the lower triangular
173: /// band part of the matrix of coefficients, supplied column by
174: /// column, with the leading diagonal of the matrix in row 1 of
175: /// the array, the first sub-diagonal starting at position 1 in
176: /// row 2, and so on. The bottom right k by k triangle of the
177: /// array A is not referenced.
178: /// The following program segment will transfer a lower
179: /// triangular band matrix from conventional full matrix storage
180: /// to band storage:
181: ///
182: /// DO 20, J = 1, N
183: /// M = 1 - J
184: /// DO 10, I = J, MIN( N, J + K )
185: /// A( M + I, J ) = matrix( I, J )
186: /// 10 CONTINUE
187: /// 20 CONTINUE
188: ///
189: /// Note that when DIAG = 'U' or 'u' the elements of the array A
190: /// corresponding to the diagonal elements of the matrix are not
191: /// referenced, but are assumed to be unity.
192: /// Unchanged on exit.
193: ///</param>
194: /// <param name="LDA">
195: /// - INTEGER.
196: /// On entry, LDA specifies the first dimension of A as declared
197: /// in the calling (sub) program. LDA must be at least
198: /// ( k + 1 ).
199: /// Unchanged on exit.
200: ///</param>
201: /// <param name="X">
202: /// - DOUBLE PRECISION array of dimension at least
203: /// ( 1 + ( n - 1 )*abs( INCX ) ).
204: /// Before entry, the incremented array X must contain the n
205: /// element right-hand side vector b. On exit, X is overwritten
206: /// with the solution vector x.
207: ///</param>
208: /// <param name="INCX">
209: /// - INTEGER.
210: /// On entry, INCX specifies the increment for the elements of
211: /// X. INCX must not be zero.
212: /// Unchanged on exit.
213: ///
214: ///</param>
215: public void Run(string UPLO, string TRANS, string DIAG, int N, int K, double[] A, int offset_a
216: , int LDA, ref double[] X, int offset_x, int INCX)
217: {
218:
219: #region Array Index Correction
220:
221: int o_a = -1 - LDA + offset_a; int o_x = -1 + offset_x;
222:
223: #endregion
224:
225:
226: #region Strings
227:
228: UPLO = UPLO.Substring(0, 1); TRANS = TRANS.Substring(0, 1); DIAG = DIAG.Substring(0, 1);
229:
230: #endregion
231:
232:
233: #region Prolog
234:
235: // * .. Scalar Arguments ..
236: // * ..
237: // * .. Array Arguments ..
238: // * ..
239: // *
240: // * Purpose
241: // * =======
242: // *
243: // * DTBSV solves one of the systems of equations
244: // *
245: // * A*x = b, or A'*x = b,
246: // *
247: // * where b and x are n element vectors and A is an n by n unit, or
248: // * non-unit, upper or lower triangular band matrix, with ( k + 1 )
249: // * diagonals.
250: // *
251: // * No test for singularity or near-singularity is included in this
252: // * routine. Such tests must be performed before calling this routine.
253: // *
254: // * Arguments
255: // * ==========
256: // *
257: // * UPLO - CHARACTER*1.
258: // * On entry, UPLO specifies whether the matrix is an upper or
259: // * lower triangular matrix as follows:
260: // *
261: // * UPLO = 'U' or 'u' A is an upper triangular matrix.
262: // *
263: // * UPLO = 'L' or 'l' A is a lower triangular matrix.
264: // *
265: // * Unchanged on exit.
266: // *
267: // * TRANS - CHARACTER*1.
268: // * On entry, TRANS specifies the equations to be solved as
269: // * follows:
270: // *
271: // * TRANS = 'N' or 'n' A*x = b.
272: // *
273: // * TRANS = 'T' or 't' A'*x = b.
274: // *
275: // * TRANS = 'C' or 'c' A'*x = b.
276: // *
277: // * Unchanged on exit.
278: // *
279: // * DIAG - CHARACTER*1.
280: // * On entry, DIAG specifies whether or not A is unit
281: // * triangular as follows:
282: // *
283: // * DIAG = 'U' or 'u' A is assumed to be unit triangular.
284: // *
285: // * DIAG = 'N' or 'n' A is not assumed to be unit
286: // * triangular.
287: // *
288: // * Unchanged on exit.
289: // *
290: // * N - INTEGER.
291: // * On entry, N specifies the order of the matrix A.
292: // * N must be at least zero.
293: // * Unchanged on exit.
294: // *
295: // * K - INTEGER.
296: // * On entry with UPLO = 'U' or 'u', K specifies the number of
297: // * super-diagonals of the matrix A.
298: // * On entry with UPLO = 'L' or 'l', K specifies the number of
299: // * sub-diagonals of the matrix A.
300: // * K must satisfy 0 .le. K.
301: // * Unchanged on exit.
302: // *
303: // * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
304: // * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
305: // * by n part of the array A must contain the upper triangular
306: // * band part of the matrix of coefficients, supplied column by
307: // * column, with the leading diagonal of the matrix in row
308: // * ( k + 1 ) of the array, the first super-diagonal starting at
309: // * position 2 in row k, and so on. The top left k by k triangle
310: // * of the array A is not referenced.
311: // * The following program segment will transfer an upper
312: // * triangular band matrix from conventional full matrix storage
313: // * to band storage:
314: // *
315: // * DO 20, J = 1, N
316: // * M = K + 1 - J
317: // * DO 10, I = MAX( 1, J - K ), J
318: // * A( M + I, J ) = matrix( I, J )
319: // * 10 CONTINUE
320: // * 20 CONTINUE
321: // *
322: // * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
323: // * by n part of the array A must contain the lower triangular
324: // * band part of the matrix of coefficients, supplied column by
325: // * column, with the leading diagonal of the matrix in row 1 of
326: // * the array, the first sub-diagonal starting at position 1 in
327: // * row 2, and so on. The bottom right k by k triangle of the
328: // * array A is not referenced.
329: // * The following program segment will transfer a lower
330: // * triangular band matrix from conventional full matrix storage
331: // * to band storage:
332: // *
333: // * DO 20, J = 1, N
334: // * M = 1 - J
335: // * DO 10, I = J, MIN( N, J + K )
336: // * A( M + I, J ) = matrix( I, J )
337: // * 10 CONTINUE
338: // * 20 CONTINUE
339: // *
340: // * Note that when DIAG = 'U' or 'u' the elements of the array A
341: // * corresponding to the diagonal elements of the matrix are not
342: // * referenced, but are assumed to be unity.
343: // * Unchanged on exit.
344: // *
345: // * LDA - INTEGER.
346: // * On entry, LDA specifies the first dimension of A as declared
347: // * in the calling (sub) program. LDA must be at least
348: // * ( k + 1 ).
349: // * Unchanged on exit.
350: // *
351: // * X - DOUBLE PRECISION array of dimension at least
352: // * ( 1 + ( n - 1 )*abs( INCX ) ).
353: // * Before entry, the incremented array X must contain the n
354: // * element right-hand side vector b. On exit, X is overwritten
355: // * with the solution vector x.
356: // *
357: // * INCX - INTEGER.
358: // * On entry, INCX specifies the increment for the elements of
359: // * X. INCX must not be zero.
360: // * Unchanged on exit.
361: // *
362: // *
363: // * Level 2 Blas routine.
364: // *
365: // * -- Written on 22-October-1986.
366: // * Jack Dongarra, Argonne National Lab.
367: // * Jeremy Du Croz, Nag Central Office.
368: // * Sven Hammarling, Nag Central Office.
369: // * Richard Hanson, Sandia National Labs.
370: // *
371: // *
372: // * .. Parameters ..
373: // * ..
374: // * .. Local Scalars ..
375: // * ..
376: // * .. External Functions ..
377: // * ..
378: // * .. External Subroutines ..
379: // * ..
380: // * .. Intrinsic Functions ..
381: // INTRINSIC MAX,MIN;
382: // * ..
383: // *
384: // * Test the input parameters.
385: // *
386:
387: #endregion
388:
389:
390: #region Body
391:
392: INFO = 0;
393: if (!this._lsame.Run(UPLO, "U") && !this._lsame.Run(UPLO, "L"))
394: {
395: INFO = 1;
396: }
397: else
398: {
399: if (!this._lsame.Run(TRANS, "N") && !this._lsame.Run(TRANS, "T") && !this._lsame.Run(TRANS, "C"))
400: {
401: INFO = 2;
402: }
403: else
404: {
405: if (!this._lsame.Run(DIAG, "U") && !this._lsame.Run(DIAG, "N"))
406: {
407: INFO = 3;
408: }
409: else
410: {
411: if (N < 0)
412: {
413: INFO = 4;
414: }
415: else
416: {
417: if (K < 0)
418: {
419: INFO = 5;
420: }
421: else
422: {
423: if (LDA < (K + 1))
424: {
425: INFO = 7;
426: }
427: else
428: {
429: if (INCX == 0)
430: {
431: INFO = 9;
432: }
433: }
434: }
435: }
436: }
437: }
438: }
439: if (INFO != 0)
440: {
441: this._xerbla.Run("DTBSV ", INFO);
442: return;
443: }
444: // *
445: // * Quick return if possible.
446: // *
447: if (N == 0) return;
448: // *
449: NOUNIT = this._lsame.Run(DIAG, "N");
450: // *
451: // * Set up the start point in X if the increment is not unity. This
452: // * will be ( N - 1 )*INCX too small for descending loops.
453: // *
454: if (INCX <= 0)
455: {
456: KX = 1 - (N - 1) * INCX;
457: }
458: else
459: {
460: if (INCX != 1)
461: {
462: KX = 1;
463: }
464: }
465: // *
466: // * Start the operations. In this version the elements of A are
467: // * accessed by sequentially with one pass through A.
468: // *
469: if (this._lsame.Run(TRANS, "N"))
470: {
471: // *
472: // * Form x := inv( A )*x.
473: // *
474: if (this._lsame.Run(UPLO, "U"))
475: {
476: KPLUS1 = K + 1;
477: if (INCX == 1)
478: {
479: for (J = N; J >= 1; J += - 1)
480: {
481: if (X[J + o_x] != ZERO)
482: {
483: L = KPLUS1 - J;
484: if (NOUNIT) X[J + o_x] = X[J + o_x] / A[KPLUS1+J * LDA + o_a];
485: TEMP = X[J + o_x];
486: for (I = J - 1; I >= Math.Max(1, J - K); I += - 1)
487: {
488: X[I + o_x] = X[I + o_x] - TEMP * A[L + I+J * LDA + o_a];
489: }
490: }
491: }
492: }
493: else
494: {
495: KX = KX + (N - 1) * INCX;
496: JX = KX;
497: for (J = N; J >= 1; J += - 1)
498: {
499: KX = KX - INCX;
500: if (X[JX + o_x] != ZERO)
501: {
502: IX = KX;
503: L = KPLUS1 - J;
504: if (NOUNIT) X[JX + o_x] = X[JX + o_x] / A[KPLUS1+J * LDA + o_a];
505: TEMP = X[JX + o_x];
506: for (I = J - 1; I >= Math.Max(1, J - K); I += - 1)
507: {
508: X[IX + o_x] = X[IX + o_x] - TEMP * A[L + I+J * LDA + o_a];
509: IX = IX - INCX;
510: }
511: }
512: JX = JX - INCX;
513: }
514: }
515: }
516: else
517: {
518: if (INCX == 1)
519: {
520: for (J = 1; J <= N; J++)
521: {
522: if (X[J + o_x] != ZERO)
523: {
524: L = 1 - J;
525: if (NOUNIT) X[J + o_x] = X[J + o_x] / A[1+J * LDA + o_a];
526: TEMP = X[J + o_x];
527: for (I = J + 1; I <= Math.Min(N, J + K); I++)
528: {
529: X[I + o_x] = X[I + o_x] - TEMP * A[L + I+J * LDA + o_a];
530: }
531: }
532: }
533: }
534: else
535: {
536: JX = KX;
537: for (J = 1; J <= N; J++)
538: {
539: KX = KX + INCX;
540: if (X[JX + o_x] != ZERO)
541: {
542: IX = KX;
543: L = 1 - J;
544: if (NOUNIT) X[JX + o_x] = X[JX + o_x] / A[1+J * LDA + o_a];
545: TEMP = X[JX + o_x];
546: for (I = J + 1; I <= Math.Min(N, J + K); I++)
547: {
548: X[IX + o_x] = X[IX + o_x] - TEMP * A[L + I+J * LDA + o_a];
549: IX = IX + INCX;
550: }
551: }
552: JX = JX + INCX;
553: }
554: }
555: }
556: }
557: else
558: {
559: // *
560: // * Form x := inv( A')*x.
561: // *
562: if (this._lsame.Run(UPLO, "U"))
563: {
564: KPLUS1 = K + 1;
565: if (INCX == 1)
566: {
567: for (J = 1; J <= N; J++)
568: {
569: TEMP = X[J + o_x];
570: L = KPLUS1 - J;
571: for (I = Math.Max(1, J - K); I <= J - 1; I++)
572: {
573: TEMP = TEMP - A[L + I+J * LDA + o_a] * X[I + o_x];
574: }
575: if (NOUNIT) TEMP = TEMP / A[KPLUS1+J * LDA + o_a];
576: X[J + o_x] = TEMP;
577: }
578: }
579: else
580: {
581: JX = KX;
582: for (J = 1; J <= N; J++)
583: {
584: TEMP = X[JX + o_x];
585: IX = KX;
586: L = KPLUS1 - J;
587: for (I = Math.Max(1, J - K); I <= J - 1; I++)
588: {
589: TEMP = TEMP - A[L + I+J * LDA + o_a] * X[IX + o_x];
590: IX = IX + INCX;
591: }
592: if (NOUNIT) TEMP = TEMP / A[KPLUS1+J * LDA + o_a];
593: X[JX + o_x] = TEMP;
594: JX = JX + INCX;
595: if (J > K) KX = KX + INCX;
596: }
597: }
598: }
599: else
600: {
601: if (INCX == 1)
602: {
603: for (J = N; J >= 1; J += - 1)
604: {
605: TEMP = X[J + o_x];
606: L = 1 - J;
607: for (I = Math.Min(N, J + K); I >= J + 1; I += - 1)
608: {
609: TEMP = TEMP - A[L + I+J * LDA + o_a] * X[I + o_x];
610: }
611: if (NOUNIT) TEMP = TEMP / A[1+J * LDA + o_a];
612: X[J + o_x] = TEMP;
613: }
614: }
615: else
616: {
617: KX = KX + (N - 1) * INCX;
618: JX = KX;
619: for (J = N; J >= 1; J += - 1)
620: {
621: TEMP = X[JX + o_x];
622: IX = KX;
623: L = 1 - J;
624: for (I = Math.Min(N, J + K); I >= J + 1; I += - 1)
625: {
626: TEMP = TEMP - A[L + I+J * LDA + o_a] * X[IX + o_x];
627: IX = IX - INCX;
628: }
629: if (NOUNIT) TEMP = TEMP / A[1+J * LDA + o_a];
630: X[JX + o_x] = TEMP;
631: JX = JX - INCX;
632: if ((N - J) >= K) KX = KX - INCX;
633: }
634: }
635: }
636: }
637: // *
638: return;
639: // *
640: // * End of DTBSV .
641: // *
642:
643: #endregion
644:
645: }
646: }
647: }