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