3 /* Subroutine */ int dgemm_(char *transa, char *transb, integer *m, integer *
4 n, integer *k, doublereal *alpha, doublereal *a, integer *lda,
5 doublereal *b, integer *ldb, doublereal *beta, doublereal *c__,
8 /* System generated locals */
9 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
13 integer i__, j, l, info;
17 extern logical lsame_(char *, char *);
19 extern /* Subroutine */ int xerbla_(char *, integer *);
21 /* .. Scalar Arguments .. */
23 /* .. Array Arguments .. */
29 /* DGEMM performs one of the matrix-matrix operations */
31 /* C := alpha*op( A )*op( B ) + beta*C, */
33 /* where op( X ) is one of */
35 /* op( X ) = X or op( X ) = X', */
37 /* alpha and beta are scalars, and A, B and C are matrices, with op( A ) */
38 /* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. */
43 /* TRANSA - CHARACTER*1. */
44 /* On entry, TRANSA specifies the form of op( A ) to be used in */
45 /* the matrix multiplication as follows: */
47 /* TRANSA = 'N' or 'n', op( A ) = A. */
49 /* TRANSA = 'T' or 't', op( A ) = A'. */
51 /* TRANSA = 'C' or 'c', op( A ) = A'. */
53 /* Unchanged on exit. */
55 /* TRANSB - CHARACTER*1. */
56 /* On entry, TRANSB specifies the form of op( B ) to be used in */
57 /* the matrix multiplication as follows: */
59 /* TRANSB = 'N' or 'n', op( B ) = B. */
61 /* TRANSB = 'T' or 't', op( B ) = B'. */
63 /* TRANSB = 'C' or 'c', op( B ) = B'. */
65 /* Unchanged on exit. */
68 /* On entry, M specifies the number of rows of the matrix */
69 /* op( A ) and of the matrix C. M must be at least zero. */
70 /* Unchanged on exit. */
73 /* On entry, N specifies the number of columns of the matrix */
74 /* op( B ) and the number of columns of the matrix C. N must be */
76 /* Unchanged on exit. */
79 /* On entry, K specifies the number of columns of the matrix */
80 /* op( A ) and the number of rows of the matrix op( B ). K must */
81 /* be at least zero. */
82 /* Unchanged on exit. */
84 /* ALPHA - DOUBLE PRECISION. */
85 /* On entry, ALPHA specifies the scalar alpha. */
86 /* Unchanged on exit. */
88 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is */
89 /* k when TRANSA = 'N' or 'n', and is m otherwise. */
90 /* Before entry with TRANSA = 'N' or 'n', the leading m by k */
91 /* part of the array A must contain the matrix A, otherwise */
92 /* the leading k by m part of the array A must contain the */
94 /* Unchanged on exit. */
97 /* On entry, LDA specifies the first dimension of A as declared */
98 /* in the calling (sub) program. When TRANSA = 'N' or 'n' then */
99 /* LDA must be at least max( 1, m ), otherwise LDA must be at */
100 /* least max( 1, k ). */
101 /* Unchanged on exit. */
103 /* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is */
104 /* n when TRANSB = 'N' or 'n', and is k otherwise. */
105 /* Before entry with TRANSB = 'N' or 'n', the leading k by n */
106 /* part of the array B must contain the matrix B, otherwise */
107 /* the leading n by k part of the array B must contain the */
109 /* Unchanged on exit. */
112 /* On entry, LDB specifies the first dimension of B as declared */
113 /* in the calling (sub) program. When TRANSB = 'N' or 'n' then */
114 /* LDB must be at least max( 1, k ), otherwise LDB must be at */
115 /* least max( 1, n ). */
116 /* Unchanged on exit. */
118 /* BETA - DOUBLE PRECISION. */
119 /* On entry, BETA specifies the scalar beta. When BETA is */
120 /* supplied as zero then C need not be set on input. */
121 /* Unchanged on exit. */
123 /* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). */
124 /* Before entry, the leading m by n part of the array C must */
125 /* contain the matrix C, except when beta is zero, in which */
126 /* case C need not be set on entry. */
127 /* On exit, the array C is overwritten by the m by n matrix */
128 /* ( alpha*op( A )*op( B ) + beta*C ). */
131 /* On entry, LDC specifies the first dimension of C as declared */
132 /* in the calling (sub) program. LDC must be at least */
134 /* Unchanged on exit. */
137 /* Level 3 Blas routine. */
139 /* -- Written on 8-February-1989. */
140 /* Jack Dongarra, Argonne National Laboratory. */
141 /* Iain Duff, AERE Harwell. */
142 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
143 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
146 /* .. External Functions .. */
148 /* .. External Subroutines .. */
150 /* .. Intrinsic Functions .. */
152 /* .. Local Scalars .. */
154 /* .. Parameters .. */
157 /* Set NOTA and NOTB as true if A and B respectively are not */
158 /* transposed and set NROWA, NCOLA and NROWB as the number of rows */
159 /* and columns of A and the number of rows of B respectively. */
161 /* Parameter adjustments */
163 a_offset = 1 + a_dim1;
166 b_offset = 1 + b_dim1;
169 c_offset = 1 + c_dim1;
173 nota = lsame_(transa, "N");
174 notb = lsame_(transb, "N");
188 /* Test the input parameters. */
191 if (! nota && ! lsame_(transa, "C") && ! lsame_(
194 } else if (! notb && ! lsame_(transb, "C") && !
195 lsame_(transb, "T")) {
203 } else if (*lda < max(1,nrowa)) {
205 } else if (*ldb < max(1,nrowb)) {
207 } else if (*ldc < max(1,*m)) {
211 xerbla_("DGEMM ", &info);
215 /* Quick return if possible. */
217 if (*m == 0 || *n == 0 || (*alpha == 0. || *k == 0) && *beta == 1.) {
221 /* And if alpha.eq.zero. */
226 for (j = 1; j <= i__1; ++j) {
228 for (i__ = 1; i__ <= i__2; ++i__) {
229 c__[i__ + j * c_dim1] = 0.;
236 for (j = 1; j <= i__1; ++j) {
238 for (i__ = 1; i__ <= i__2; ++i__) {
239 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
248 /* Start the operations. */
253 /* Form C := alpha*A*B + beta*C. */
256 for (j = 1; j <= i__1; ++j) {
259 for (i__ = 1; i__ <= i__2; ++i__) {
260 c__[i__ + j * c_dim1] = 0.;
263 } else if (*beta != 1.) {
265 for (i__ = 1; i__ <= i__2; ++i__) {
266 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
271 for (l = 1; l <= i__2; ++l) {
272 if (b[l + j * b_dim1] != 0.) {
273 temp = *alpha * b[l + j * b_dim1];
275 for (i__ = 1; i__ <= i__3; ++i__) {
276 c__[i__ + j * c_dim1] += temp * a[i__ + l *
287 /* Form C := alpha*A'*B + beta*C */
290 for (j = 1; j <= i__1; ++j) {
292 for (i__ = 1; i__ <= i__2; ++i__) {
295 for (l = 1; l <= i__3; ++l) {
296 temp += a[l + i__ * a_dim1] * b[l + j * b_dim1];
300 c__[i__ + j * c_dim1] = *alpha * temp;
302 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[
313 /* Form C := alpha*A*B' + beta*C */
316 for (j = 1; j <= i__1; ++j) {
319 for (i__ = 1; i__ <= i__2; ++i__) {
320 c__[i__ + j * c_dim1] = 0.;
323 } else if (*beta != 1.) {
325 for (i__ = 1; i__ <= i__2; ++i__) {
326 c__[i__ + j * c_dim1] = *beta * c__[i__ + j * c_dim1];
331 for (l = 1; l <= i__2; ++l) {
332 if (b[j + l * b_dim1] != 0.) {
333 temp = *alpha * b[j + l * b_dim1];
335 for (i__ = 1; i__ <= i__3; ++i__) {
336 c__[i__ + j * c_dim1] += temp * a[i__ + l *
347 /* Form C := alpha*A'*B' + beta*C */
350 for (j = 1; j <= i__1; ++j) {
352 for (i__ = 1; i__ <= i__2; ++i__) {
355 for (l = 1; l <= i__3; ++l) {
356 temp += a[l + i__ * a_dim1] * b[j + l * b_dim1];
360 c__[i__ + j * c_dim1] = *alpha * temp;
362 c__[i__ + j * c_dim1] = *alpha * temp + *beta * c__[