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