Question about the use of GEMM vs MATMUL

[isyegatech]

Hi everyone, I am trying to finding the root of a nonlinear system using the fortran 90 subroutine "newt" provided in the book, Numerical Recipes.

I try to substitute the statement [highlight #D3D7CF]GEMM[/highlight], provided in Math Kernel Library, for the fortran built-in statement [highlight #D3D7CF]MATMUL[/highlight]: that is,


replacing [highlight #FCE94F]g)=matmul(fvec),fjac,)[/highlight]
with [highlight #FCE94F]call gemm(fvec),fjac,,g))[/highlight].

However, this change leads to a compiling error that says "Could not resolve generic procedure gemm." I am sure the linking setup in my visual studio's properties page is correct, because in the meantime the use of other MKL functions works normally.

The [highlight #D3D7CF]fjac[/highlight] is a generic real type array used for the Jacobian matrix. The [highlight #D3D7CF]fvec[/highlight] is a target to which a pointer [highlight #D3D7CF]fmin_fvecp[/highlight] points. The pointer is declared in a module's internal function [highlight #D3D7CF]fmin[/highlight]. The module and the associated external function [highlight #D3D7CF]funcv[/highlight] that defines the nonlinear system are shown as follows:

MODULE model
(...)
REAL(WP), DIMENSION), POINTER :: fmin_fvecp
CONTAINS
FUNCTION fmin(x)
IMPLICIT NONE
REAL(WP), DIMENSION), INTENT(IN) :: x
REAL(WP) :: fmin
INTERFACE
FUNCTION funcv(x)
USE mkl95_precision, ONLY: WP => DP
IMPLICIT NONE
REAL(WP), DIMENSION), INTENT(IN) :: x
REAL(WP), DIMENSION(size(x)) :: funcv
END FUNCTION funcv
END INTERFACE
if (.not. associated(fmin_fvecp)) call &
nrerror('fmin: problem with pointer for returned values')
fmin_fvecp=funcv(x)
fmin=0.5_wp*dot(fmin_fvecp,fmin_fvecp)
END FUNCTION fmin
(...)
END MODULE model

FUNCTION funcv(x)
USE mkl95_precision, ONLY: WP => DP
IMPLICIT NONE
REAL(WP), DIMENSION), INTENT(IN) :: x
REAL(WP), DIMENSION(size(x)) :: funcv
funcv(1)=x(1)**2+x(2)**2-2.0_sp
funcv(2)=exp(x(1)-1.0_wp)+x(2)**3-2.0_wp
END FUNCTION funcv

Please let me know what's wrong with the statement call gemm(fvec),fjac,,g)). Thank you for the help.

Lee

 

[gummibaer]

I do not know about this gemm-routine.

But I would set up a small program that just uses this routine to do the operation you want. Say multiply two 3 by three matrices to see how it works.
As Germán would say: take baby-steps, one at a time.

Norbert

The optimist believes we live in the best of all possible worlds - the pessimist fears this might be true.

 

[isyegatech]

After several experiments on the inputs of GEMM, I find the compiler accepts either the statement GEMM(a,,bb,,c,) or GEMM(a,bb,c) as long as the congruity of dimension of input variables holds. So, the compiling error comes my mistaking a) with a,, which hence violates the congruity of matrix multiplication.

What surprises me is that the build-in math function MATMUL accept the loose expression where, for example, in the matrix multiplication of n X 1 A array with n X n array B, the n X 1 array A can reduce to array of length n C, but not in GEMM: that is, suppose array ANS1 is of size n and ANS2 is of size n X 1. It is legal using ANS1=MATMUL(C,B) instead of ANS2=MATMUL(A,B). GEMM only accepts the second form of statement.

 

[salgerman]

Norbert: you beat me to it.

Lee: Yes, with c = matmul(a,b), three cases are possible:

1.- c(n,k) = matmul( a(n,m) b(m,k) ) ; where c(i,j) = sum( a(i, * b,j) )

2.- c(k) = matmul( a(m) b(m,k) ) ; where c(j) = sum( a*b,j) )

3.- c = matmul( a(n,m) b(m) ) ; where c(i) = sum( a(i,*b )

where the dimensions in the expression on the left have been added for illustration purposes of the expected and resulting matrix dimensions.

Germán

 

[mikrom]

The subbroutine GEMM of Math Kernel Library seems to be the GEMM of the BLAS package -
see page 2-84 of the document

http://www.csar.cfs.ac.uk/user_information/software/maths/mklman61.pdf

So if it's GEMM of BLAS then you provide wrong arguments.

For example for double precision subroutine the call should be:

call DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)

where the arguments are described in the source of subroutine:

http://www.netlib.org/blas/dgemm.f

Here are other infos about GEMM

http://en.wikipedia.org/wiki/General_Matrix_Multiply

http://software.intel.com/sites/pro...GUID-36BFBCE9-EB0A-43B0-ADAF-2B65275726EA.htm

 

[isyegatech]

Hi everyone, thank you for all the replies. I've figured out what's wrong with my statement now.

@mikrom: the simplified GEMM statement is legal as long as the compiler knows where the pre-compiled subroutine package is located (that is, BLAS95). Otherwise, as you suggested, for real-type variables of double precision, we have to use DGEMM to perform matrix multiplication.

 

[mikrom]

Don't understand what you mean with "the simplified GEMM statement".
AFAIK, xGEMM are BLAS-subroutines not statement.

 

[isyegatech]

@mikrom: Sorry that I didn't make it clear in my last post. I mean, if people use BLAS95.mod, then they can use the simplified statement, GEMM, instead of the original, a tad lengthy statements like DGEMM and SGEMM for which more input arguments are required. I totally agree with the information you provided. Thanks.

 

[mikrom]

Hi isyegatech,
Could you please show me, how the usage of the simplified GEMM looks like?
I'm curios how many arguments it has - instead of the 13 arguments, which are necessary for conventional subroutine.

 

[isyegatech]

@mikrom: please take a look at the following link．You would see the difference in the use of GEMM (via Math Kernel Library wrapped up by Intel) and xGEMM. Put briefly, the shortest form could be "GEMM(a,b,c)" as opposed to "call dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc).

 

[mikrom]

Very nice, the usage of the simplified GEMM looks more human than the original GEMM with 13 arguments.
Thank you very much.