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