3 /* Subroutine */ int sgemm_(char *transa, char *transb, integer *m, integer *
4 n, integer *k, real *alpha, real *a, integer *lda, real *b, integer *
5 ldb, real *beta, real *c__, integer *ldc)
7 /* System generated locals */
8 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
12 integer i__, j, l, info;
16 extern logical lsame_(char *, char *);
18 extern /* Subroutine */ int xerbla_(char *, integer *);
20 /* .. Scalar Arguments .. */
22 /* .. Array Arguments .. */
28 /* SGEMM performs one of the matrix-matrix operations */
30 /* C := alpha*op( A )*op( B ) + beta*C, */
32 /* where op( X ) is one of */
34 /* op( X ) = X or op( X ) = X', */
36 /* alpha and beta are scalars, and A, B and C are matrices, with op( A ) */
37 /* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. */
42 /* TRANSA - CHARACTER*1. */
43 /* On entry, TRANSA specifies the form of op( A ) to be used in */
44 /* the matrix multiplication as follows: */
46 /* TRANSA = 'N' or 'n', op( A ) = A. */
48 /* TRANSA = 'T' or 't', op( A ) = A'. */
50 /* TRANSA = 'C' or 'c', op( A ) = A'. */
52 /* Unchanged on exit. */
54 /* TRANSB - CHARACTER*1. */
55 /* On entry, TRANSB specifies the form of op( B ) to be used in */
56 /* the matrix multiplication as follows: */
58 /* TRANSB = 'N' or 'n', op( B ) = B. */
60 /* TRANSB = 'T' or 't', op( B ) = B'. */
62 /* TRANSB = 'C' or 'c', op( B ) = B'. */
64 /* Unchanged on exit. */
67 /* On entry, M specifies the number of rows of the matrix */
68 /* op( A ) and of the matrix C. M must be at least zero. */
69 /* Unchanged on exit. */
72 /* On entry, N specifies the number of columns of the matrix */
73 /* op( B ) and the number of columns of the matrix C. N must be */
75 /* Unchanged on exit. */
78 /* On entry, K specifies the number of columns of the matrix */
79 /* op( A ) and the number of rows of the matrix op( B ). K must */
80 /* be at least zero. */
81 /* Unchanged on exit. */
84 /* On entry, ALPHA specifies the scalar alpha. */
85 /* Unchanged on exit. */
87 /* A - REAL array of DIMENSION ( LDA, ka ), where ka is */
88 /* k when TRANSA = 'N' or 'n', and is m otherwise. */
89 /* Before entry with TRANSA = 'N' or 'n', the leading m by k */
90 /* part of the array A must contain the matrix A, otherwise */
91 /* the leading k by m part of the array A must contain the */
93 /* Unchanged on exit. */
96 /* On entry, LDA specifies the first dimension of A as declared */
97 /* in the calling (sub) program. When TRANSA = 'N' or 'n' then */
98 /* LDA must be at least max( 1, m ), otherwise LDA must be at */
99 /* least max( 1, k ). */
100 /* Unchanged on exit. */
102 /* B - REAL array of DIMENSION ( LDB, kb ), where kb is */
103 /* n when TRANSB = 'N' or 'n', and is k otherwise. */
104 /* Before entry with TRANSB = 'N' or 'n', the leading k by n */
105 /* part of the array B must contain the matrix B, otherwise */
106 /* the leading n by k part of the array B must contain the */
108 /* Unchanged on exit. */
111 /* On entry, LDB specifies the first dimension of B as declared */
112 /* in the calling (sub) program. When TRANSB = 'N' or 'n' then */
113 /* LDB must be at least max( 1, k ), otherwise LDB must be at */
114 /* least max( 1, n ). */
115 /* Unchanged on exit. */
118 /* On entry, BETA specifies the scalar beta. When BETA is */
119 /* supplied as zero then C need not be set on input. */
120 /* Unchanged on exit. */
122 /* C - REAL array of DIMENSION ( LDC, n ). */
123 /* Before entry, the leading m by n part of the array C must */
124 /* contain the matrix C, except when beta is zero, in which */
125 /* case C need not be set on entry. */
126 /* On exit, the array C is overwritten by the m by n matrix */
127 /* ( alpha*op( A )*op( B ) + beta*C ). */
130 /* On entry, LDC specifies the first dimension of C as declared */
131 /* in the calling (sub) program. LDC must be at least */
133 /* Unchanged on exit. */
136 /* Level 3 Blas routine. */
138 /* -- Written on 8-February-1989. */
139 /* Jack Dongarra, Argonne National Laboratory. */
140 /* Iain Duff, AERE Harwell. */
141 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
142 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
145 /* .. External Functions .. */
147 /* .. External Subroutines .. */
149 /* .. Intrinsic Functions .. */
151 /* .. Local Scalars .. */
153 /* .. Parameters .. */
156 /* Set NOTA and NOTB as true if A and B respectively are not */
157 /* transposed and set NROWA, NCOLA and NROWB as the number of rows */
158 /* and columns of A and the number of rows of B respectively. */
160 /* Parameter adjustments */
162 a_offset = 1 + a_dim1;
165 b_offset = 1 + b_dim1;
168 c_offset = 1 + c_dim1;
172 nota = lsame_(transa, "N");
173 notb = lsame_(transb, "N");
187 /* Test the input parameters. */
190 if (! nota && ! lsame_(transa, "C") && ! lsame_(
193 } else if (! notb && ! lsame_(transb, "C") && !
194 lsame_(transb, "T")) {
202 } else if (*lda < max(1,nrowa)) {
204 } else if (*ldb < max(1,nrowb)) {
206 } else if (*ldc < max(1,*m)) {
210 xerbla_("SGEMM ", &info);
214 /* Quick return if possible. */
216 if (*m == 0 || *n == 0 || (*alpha == 0.f || *k == 0) && *beta == 1.f) {
220 /* And if alpha.eq.zero. */
225 for (j = 1; j <= i__1; ++j) {
227 for (i__ = 1; i__ <= i__2; ++i__) {
228 c__[i__ + j * c_dim1] = 0.f;
235 for (j = 1; j <= i__1; ++j) {
237 for (i__ = 1; i__ <= i__2; ++i__) {
238 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
247 /* Start the operations. */
252 /* Form C := alpha*A*B + beta*C. */
255 for (j = 1; j <= i__1; ++j) {
258 for (i__ = 1; i__ <= i__2; ++i__) {
259 c__[i__ + j * c_dim1] = 0.f;
262 } else if (*beta != 1.f) {
264 for (i__ = 1; i__ <= i__2; ++i__) {
265 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
270 for (l = 1; l <= i__2; ++l) {
271 if (b[l + j * b_dim1] != 0.f) {
272 temp = *alpha * b[l + j * b_dim1];
274 for (i__ = 1; i__ <= i__3; ++i__) {
275 c__[i__ + j * c_dim1] += temp * a[i__ + l *
286 /* Form C := alpha*A'*B + beta*C */
289 for (j = 1; j <= i__1; ++j) {
291 for (i__ = 1; i__ <= i__2; ++i__) {
294 for (l = 1; l <= i__3; ++l) {
295 temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
299 c__[i__ + j * c_dim1] = *alpha * temp;
301 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
312 /* Form C := alpha*A*B' + beta*C */
315 for (j = 1; j <= i__1; ++j) {
318 for (i__ = 1; i__ <= i__2; ++i__) {
319 c__[i__ + j * c_dim1] = 0.f;
322 } else if (*beta != 1.f) {
324 for (i__ = 1; i__ <= i__2; ++i__) {
325 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
330 for (l = 1; l <= i__2; ++l) {
331 if (b[j + l * b_dim1] != 0.f) {
332 temp = *alpha * b[j + l * b_dim1];
334 for (i__ = 1; i__ <= i__3; ++i__) {
335 c__[i__ + j * c_dim1] += temp * a[i__ + l *
346 /* Form C := alpha*A'*B' + beta*C */
349 for (j = 1; j <= i__1; ++j) {
351 for (i__ = 1; i__ <= i__2; ++i__) {
354 for (l = 1; l <= i__3; ++l) {
355 temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
359 c__[i__ + j * c_dim1] = *alpha * temp;
361 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[