SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)
USE NUMERICAL_LIBRARIES
IMPLICIT NONE
INTEGER PLHS(*), PRHS(*)
INTEGER NLHS, NRHS
INTEGER MXCREATEFULL, MXGETPR
INTEGER MXGETM, MXGETN

integer M(2), N(2)

real(8), allocatable :: A(:,:), B(:,:), C(:,:)

IF (NRHS .NE. 2) THEN
	CALL MEXERRMSGTXT('abcg requires two input arguments')
ELSEIF (NLHS .GT. 1) THEN
	CALL MEXERRMSGTXT('abcg requires one output argument')
ENDIF

M(1) = MXGETM(PRHS(1))
N(1) = MXGETN(PRHS(1))
M(2) = MXGETM(PRHS(2))
N(2) = MXGETN(PRHS(2))

ALLOCATE(A(M(1),N(1)))
ALLOCATE(B(M(2),N(2)))
ALLOCATE(C(M(1),N(2)))

call mxcopyptrtoreal8(mxgetpr(prhs(1)),A,M(1)*N(1))
call mxcopyptrtoreal8(mxgetpr(prhs(2)),B,M(2)*N(2))

plhs(1) = mxcreatefull(M(1),N(2),0)


call dgemm ( 'n', 'n', m(1), n(1), m(2), 1.0d0, A, m(1), &
B, m(2), 0.0d0, c, m(1) )

!C = matmul(A,B)

CALL MXCOPYreal8TOPTR(C,mxgetpr(plhs(1)),M(1)*N(2))
RETURN
END