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