about matmul()

[IDZ_A_Intel]

Hi,

My question is about the function matmul(). As we know, matmul() is used to multiply two matrixs, A, B. However, I have the question that how to do when I want to multiply one row of A and one column of B? For example, I define three matrixs, A, B and C

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

allocate(A(5,3),B(3,4),C(1,1))

Now, if I run

C=matmul(A(2,:),B(:,2))

I get the error message,

test1.f90(40): error #6366: The shapes of the array expressions do not conform.
C=matmul(A(2,:),B(:,2))

However, if I use

C=matmul(reshape(A(2,:),(/1,3/)),reshape(B(:,2),(/3,1/)))

There is no problem. Moreover, I find if I have two matrix, A(a,b) and B(c,d), I must reshape the sub-matrix when either a or b or c or d=1. This is very boring (I must use reshape() very time). Why cannot the compiler automatically detect the shape of sub-matrix in matmul()? In matlab, it is very easy, and I just use C=A(2,:)*B(:,2).

Thanks very much.

[IDZ_A_Intel]

This is simply the way the Fortran standard defines the intrinsic MATMUL, which is not really a general purpose matrix multiply. A(2,:) ad B(:,2) are rank 1 arrays, and the standard says that if either of the arguments are rank 1, the other must be rank 2. This is why you have to use reshape, and you have to reshape both to rank 2 because otherwise the result will be rank 1 and won't conform to C.

You might want to look at the matrix multiply routines in Intel Math Kernel Library to see if they better meet your needs (though for small arrays they will be slower than using MATMUL.)

[IDZ_A_Intel]

Hi Steve,

Could I ask what matrix multiply routines the MKL have provided? Do you mean BLAS and LAPACK routines? Thanks very much.

Ying

[IDZ_A_Intel]

I was thinking of GEMV in MKL, but I'll admit that this is a topic I'm not intimately familiar with. I'm sure there are others here, such as tim18, who can probably offer more cogent advice.

[TimP]

MKL supports BLAS functions ?GEMM for rank 2 by rank 2 multiplication (including OpenMP parallelization and much additional optimization for Xeon), and ?GEMV for rank 2 by rank 1 multiplication.

For rank 1 by rank 1 multiplication, DOT_PRODUCT usually is optimized better by ifort than what is possible with BLAS ?DOT. Also, for problems too small to benefit from threading, MATMUL, compiled at -O3, is often better than BLAS ?GEMV or ?GEMM.

It may be possible to write a Fortran generic interface replacement for MATMUL which could recognize various BLAS cases and select among BLAS95 GEMM, GEMV, and DOT. Not having seen it done, I suppose there are more difficulties than meet the eye.

[IDZ_A_Intel]

Hi tim18,

Yes, that is what I want if there could be a generic interface replacement for MATMUL which could recognize various BLAS cases. That will be very useful, because when I see the user manual for BLAS functions, I find the functions need many many parameters. How can I remember all the functions with all the parameters? Besides, there has not been a fortran IDE to remind me. It should be easier for C++, because VS2008+visual Assistant is very comfortable.

Anyway, thanks for your very helpful information. I am using -O3, which is very fast to run the code. And I will try openMP to make use of the rest of cores in my PC.

Ying

[Jugoslav_Dujic · Accepted solution]

Quoting Steve Lionel (Intel)

This is simply the way the Fortran standard defines the intrinsic MATMUL, which is not really a general purpose matrix multiply. A(2,:) ad B(:,2) are rank 1 arrays, and the standard says that if either of the arguments are rank 1, the other must be rank 2. This is why you have to use reshape, and you have to reshape both to rank 2 because otherwise the result will be rank 1 and won't conform to C.

You might want to look at the matrix multiply routines in Intel Math Kernel Library to see if they better meet your needs (though for small arrays they will be slower than using MATMUL.)

But you missed a much simpler way to make a rank-2 column- (or row-) matrix than RESHAPE; namely, if X(:) is an array, then X(1) is a scalar, but X(1:1) is still an array (and, for that matter, X(1:0) is also a zero-sized array). Thus:

[fortran]real, allocatable:: A(:,:), B(:,:),C(:,:)

allocate(A(5,3),B(3,4),C(1,1))

C=matmul(A(2:2,:),B(:,2:2))
[/fortran]

works as expected.

[IDZ_A_Intel]

Yes, it is a very good idea. Very simple but workable. Thanks very much.