3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
13 * .. Scalar Arguments ..
14 * DOUBLE PRECISION ALPHA,BETA
15 * INTEGER K,LDA,LDB,LDC,M,N
16 * CHARACTER TRANSA,TRANSB
18 * .. Array Arguments ..
19 * DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
28 *> DGEMM 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**T,
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.
45 *> TRANSA is CHARACTER*1
46 *> On entry, TRANSA specifies the form of op( A ) to be used in
47 *> the matrix multiplication as follows:
49 *> TRANSA = 'N' or 'n', op( A ) = A.
51 *> TRANSA = 'T' or 't', op( A ) = A**T.
53 *> TRANSA = 'C' or 'c', op( A ) = A**T.
58 *> TRANSB is CHARACTER*1
59 *> On entry, TRANSB specifies the form of op( B ) to be used in
60 *> the matrix multiplication as follows:
62 *> TRANSB = 'N' or 'n', op( B ) = B.
64 *> TRANSB = 'T' or 't', op( B ) = B**T.
66 *> TRANSB = 'C' or 'c', op( B ) = B**T.
72 *> On entry, M specifies the number of rows of the matrix
73 *> op( A ) and of the matrix C. M must be at least zero.
79 *> On entry, N specifies the number of columns of the matrix
80 *> op( B ) and the number of columns of the matrix C. N must be
87 *> On entry, K specifies the number of columns of the matrix
88 *> op( A ) and the number of rows of the matrix op( B ). K must
94 *> ALPHA is DOUBLE PRECISION.
95 *> On entry, ALPHA specifies the scalar alpha.
100 *> A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
101 *> k when TRANSA = 'N' or 'n', and is m otherwise.
102 *> Before entry with TRANSA = 'N' or 'n', the leading m by k
103 *> part of the array A must contain the matrix A, otherwise
104 *> the leading k by m part of the array A must contain the
111 *> On entry, LDA specifies the first dimension of A as declared
112 *> in the calling (sub) program. When TRANSA = 'N' or 'n' then
113 *> LDA must be at least max( 1, m ), otherwise LDA must be at
114 *> least max( 1, k ).
119 *> B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
120 *> n when TRANSB = 'N' or 'n', and is k otherwise.
121 *> Before entry with TRANSB = 'N' or 'n', the leading k by n
122 *> part of the array B must contain the matrix B, otherwise
123 *> the leading n by k part of the array B must contain the
130 *> On entry, LDB specifies the first dimension of B as declared
131 *> in the calling (sub) program. When TRANSB = 'N' or 'n' then
132 *> LDB must be at least max( 1, k ), otherwise LDB must be at
133 *> least max( 1, n ).
138 *> BETA is DOUBLE PRECISION.
139 *> On entry, BETA specifies the scalar beta. When BETA is
140 *> supplied as zero then C need not be set on input.
145 *> C is DOUBLE PRECISION array of DIMENSION ( LDC, n ).
146 *> Before entry, the leading m by n part of the array C must
147 *> contain the matrix C, except when beta is zero, in which
148 *> case C need not be set on entry.
149 *> On exit, the array C is overwritten by the m by n matrix
150 *> ( alpha*op( A )*op( B ) + beta*C ).
156 *> On entry, LDC specifies the first dimension of C as declared
157 *> in the calling (sub) program. LDC must be at least
164 *> \author Univ. of Tennessee
165 *> \author Univ. of California Berkeley
166 *> \author Univ. of Colorado Denver
169 *> \date November 2015
171 *> \ingroup double_blas_level3
173 *> \par Further Details:
174 * =====================
178 *> Level 3 Blas routine.
180 *> -- Written on 8-February-1989.
181 *> Jack Dongarra, Argonne National Laboratory.
182 *> Iain Duff, AERE Harwell.
183 *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
184 *> Sven Hammarling, Numerical Algorithms Group Ltd.
187 * =====================================================================
188 SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
190 * -- Reference BLAS level3 routine (version 3.6.0) --
191 * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
192 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195 * .. Scalar Arguments ..
196 DOUBLE PRECISION ALPHA,BETA
197 INTEGER K,LDA,LDB,LDC,M,N
198 CHARACTER TRANSA,TRANSB
200 * .. Array Arguments ..
201 DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
204 * =====================================================================
206 * .. External Functions ..
210 * .. External Subroutines ..
213 * .. Intrinsic Functions ..
216 * .. Local Scalars ..
217 DOUBLE PRECISION TEMP
218 INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
222 DOUBLE PRECISION ONE,ZERO
223 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
226 * Set NOTA and NOTB as true if A and B respectively are not
227 * transposed and set NROWA, NCOLA and NROWB as the number of rows
228 * and columns of A and the number of rows of B respectively.
230 NOTA = LSAME(TRANSA,'N')
231 NOTB = LSAME(TRANSB,'N')
245 * Test the input parameters.
248 IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
249 + (.NOT.LSAME(TRANSA,'T'))) THEN
251 ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
252 + (.NOT.LSAME(TRANSB,'T'))) THEN
254 ELSE IF (M.LT.0) THEN
256 ELSE IF (N.LT.0) THEN
258 ELSE IF (K.LT.0) THEN
260 ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
262 ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
264 ELSE IF (LDC.LT.MAX(1,M)) THEN
268 CALL XERBLA('DGEMM ',INFO)
272 * Quick return if possible.
274 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
275 + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
277 * And if alpha.eq.zero.
279 IF (ALPHA.EQ.ZERO) THEN
280 IF (BETA.EQ.ZERO) THEN
296 * Start the operations.
301 * Form C := alpha*A*B + beta*C.
304 IF (BETA.EQ.ZERO) THEN
308 ELSE IF (BETA.NE.ONE) THEN
316 C(I,J) = C(I,J) + TEMP*A(I,L)
322 * Form C := alpha*A**T*B + beta*C
328 TEMP = TEMP + A(L,I)*B(L,J)
330 IF (BETA.EQ.ZERO) THEN
333 C(I,J) = ALPHA*TEMP + BETA*C(I,J)
341 * Form C := alpha*A*B**T + beta*C
344 IF (BETA.EQ.ZERO) THEN
348 ELSE IF (BETA.NE.ONE) THEN
356 C(I,J) = C(I,J) + TEMP*A(I,L)
362 * Form C := alpha*A**T*B**T + beta*C
368 TEMP = TEMP + A(L,I)*B(J,L)
370 IF (BETA.EQ.ZERO) THEN
373 C(I,J) = ALPHA*TEMP + BETA*C(I,J)