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