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