PROGRAM CALL_SGEMM IMPLICIT NONE INTEGER, PARAMETER :: M=4,N=3,K=2 INTEGER, PARAMETER :: LDA=M+1, LDB=K+1, LDC=M+1 REAL ALPHA, BETA REAL A(LDA,K), B(LDB,N), C(LDC,N), D(LDC,N) ! A is M by K, B is K by N and C is M by N. INTEGER I, J INTERFACE SUBROUTINE SGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) IMPLICIT NONE ! .. Scalar Arguments .. REAL ALPHA,BETA INTEGER K,LDA,LDB,LDC,M,N CHARACTER(LEN=*) TRANSA,TRANSB ! .. ! .. Array Arguments .. REAL A(LDA,*),B(LDB,*),C(LDC,*) END END INTERFACE DO J=1,K DO I=1,M A(I,J)=real(I)+real(J)/10. END DO END DO DO J=1,N DO I=1,K B(I,J)=real(I)+real(J)/10. END DO END DO ALPHA=0.5 BETA=0.1 C=0.01 D(1:M,1:N)=alpha*matmul(A(1:M,1:K),B(1:K,1:N))+beta*C(1:M,1:N) CALL SGEMM('No','no',M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) D(1:M,1:N)=D(1:M,1:N)-C(1:M,1:N) ! Should have D(1:M,1:N) = C(1:M,1:N) except for rounding errors. WRITE(*,*) 'Largest difference: ', MAXVAL(ABS(D(1:M,1:N))) WRITE(*,*) 'Relative error: ', MAXVAL(ABS(D(1:M,1:N)))/MAXVAL(ABS(C(1:M,1:N))) END PROGRAM ! OUTPUT: ! Values of adjoint flags TRANSA, TRANSB = N,n ! Values of sizes m,n,k,lda,ldb and ldc = 4,3,2,5,3,4 ! Values of alpha, beta = 5.000000e-001, 1.000000e-001 !Value of i,j and A(i,j) = 1, 1, 1.100000e+000 !Value of i,j and A(i,j) = 2, 1, 2.100000e+000 !Value of i,j and A(i,j) = 3, 1, 3.100000e+000 !Value of i,j and A(i,j) = 4, 1, 4.100000e+000 !Value of i,j and A(i,j) = 1, 2, 1.200000e+000 !Value of i,j and A(i,j) = 2, 2, 2.200000e+000 !Value of i,j and A(i,j) = 3, 2, 3.200000e+000 !Value of i,j and A(i,j) = 4, 2, 4.200000e+000 ! !Value of i,j and B(i,j) = 1, 1, 1.100000e+000 !Value of i,j and B(i,j) = 2, 1, 2.100000e+000 !Value of i,j and B(i,j) = 1, 2, 1.200000e+000 !Value of i,j and B(i,j) = 2, 2, 2.200000e+000 !Value of i,j and B(i,j) = 1, 3, 1.300000e+000 !Value of i,j and B(i,j) = 2, 3, 2.300000e+000 !Value of i,j and C(i,j) = 1, 1, 1.866000e+00 !Value of i,j and C(i,j) = 2, 1, 3.466000e+00 !Value of i,j and C(i,j) = 3, 1, 5.066000e+00 !Value of i,j and C(i,j) = 4, 1, 6.665999e+00 !Value of i,j and C(i,j) = 1, 2, 1.981000e+00 !Value of i,j and C(i,j) = 2, 2, 3.681000e+00 !Value of i,j and C(i,j) = 3, 2, 5.381001e+00 !Value of i,j and C(i,j) = 4, 2, 7.081000e+00 !Value of i,j and C(i,j) = 1, 3, 2.096000e+00 !Value of i,j and C(i,j) = 2, 3, 3.896000e+00 !Value of i,j and C(i,j) = 3, 3, 5.696000e+00 !Value of i,j and C(i,j) = 4, 3, 7.495999e+00 ! Largest difference: 4.7683715E-7 ! Relative error: 6.361222E-8 ! These values are implementation dependent because ! of rounding effects. ! ! Purpose of SGEMM ! ================ ! SGEMM performs one of the matrix-matrix operations ! C := alpha*op( A )*op( B ) + beta*C, ! where op( X ) is one of ! op( X ) = X or op( X ) = X', ! alpha and beta are scalars, and A, B and C are matrices, with op( A ) ! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. ! Arguments of SGEMM ! ================== ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! TRANSA = 'N' or 'n', op( A ) = A. ! TRANSA = 'T' or 't', op( A ) = A'. ! TRANSA = 'C' or 'c', op( A ) = A'. ! Unchanged on exit. ! TRANSB - CHARACTER*1. ! On entry, TRANSB specifies the form of op( B ) to be used in ! the matrix multiplication as follows: ! TRANSB = 'N' or 'n', op( B ) = B. ! TRANSB = 'T' or 't', op( B ) = B'. ! TRANSB = 'C' or 'c', op( B ) = B'. ! Unchanged on exit. ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix ! op( A ) and of the matrix C. M must be at least zero. ! Unchanged on exit. ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix ! op( B ) and the number of columns of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! K - INTEGER. ! On entry, K specifies the number of columns of the matrix ! op( A ) and the number of rows of the matrix op( B ). K must ! be at least zero. ! Unchanged on exit. ! ALPHA - REAL . ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! A - REAL array of DIMENSION ( LDA, ka ), where ka is ! k when TRANSA = 'N' or 'n', and is m otherwise. ! Before entry with TRANSA = 'N' or 'n', the leading m by k ! part of the array A must contain the matrix A, otherwise ! the leading k by m part of the array A must contain the ! matrix A. ! Unchanged on exit. ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANSA = 'N' or 'n' then ! LDA must be at least max( 1, m ), otherwise LDA must be at ! least max( 1, k ). ! Unchanged on exit. ! B - REAL array of DIMENSION ( LDB, kb ), where kb is ! n when TRANSB = 'N' or 'n', and is k otherwise. ! Before entry with TRANSB = 'N' or 'n', the leading k by n ! part of the array B must contain the matrix B, otherwise ! the leading n by k part of the array B must contain the ! matrix B. ! Unchanged on exit. ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. When TRANSB = 'N' or 'n' then ! LDB must be at least max( 1, k ), otherwise LDB must be at ! least max( 1, n ). ! Unchanged on exit. ! BETA - REAL . ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then C need not be set on input. ! Unchanged on exit. ! C - REAL array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n matrix ! ( alpha*op( A )*op( B ) + beta*C ). ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! Level 3 Blas routine.