Fortran array rank for matmul intrinsic

Question

The following link https://gcc.gnu.org/onlinedocs/gcc-5.4.0/gfortran/MATMUL.html clearly states that gfortran expects matrices input to matmul to be of rank 1 OR 2. However the following snippet wont compile:

Program scratch
  real(kind=8) :: A(10)=(/0,1,2,3,4,5,6,7,8,9/)
  real(kind=8) :: B(10)=(/0,1,2,3,4,5,6,7,8,9/)
  real(kind=8) :: C(10,10)
  print *,rank(A),rank(B)
  C=matmul(A,B)  
End Program scratch

gfortran gives the error:

$gfortran scratch.f90 
scratch.f90:6:13:

   C=matmul(A,B)
         1
Error: ‘matrix_b’ argument of ‘matmul’ intrinsic at (1) must be of rank 2

My gfortran is 5.4.0 (compatible with the link above). Am I doing something really stupid?

@VladimirF Thanks for that, unfortunately I changed the declarations to `real(kind=4)` and gfortran reports the same error. — Clinton Winant, Jun 02 '17 at 15:05
That was not the point, see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter and https://stackoverflow.com/documentation/fortran/939/data-types/4390/precision-of-floating-point-numbers#t=20170602150748974393 and my answer. It has nothing to do with the reported error. — Vladimir F Героям слава, Jun 02 '17 at 15:07
Very closely related, but not an *exact* duplicate https://stackoverflow.com/questions/22305554/vector-multiplication-using-matmul-in-fortran — Vladimir F Героям слава, Jun 02 '17 at 15:18

score 2 · Accepted Answer · answered Jun 02 '17 at 16:07

2

You can use RESHAPE to get them into a form MATMUL will like:

Program scratch
  real(kind=8) :: A(10)=(/0,1,2,3,4,5,6,7,8,9/)
  real(kind=8) :: B(10)=(/0,1,2,3,4,5,6,7,8,9/)
  real(kind=8) :: C(10,10)
  print *,rank(A),rank(B)
  C = matmul( RESHAPE(A,(/10,1/)), RESHAPE(B,(/1,10/)) )
  WRITE(*,"(10F7.2)") C
End Program scratch

answered Jun 02 '17 at 16:07

Jack

5,801
1
15
20

Thanks, that is helpful, I was really wondering about the disconnect between the gfortran documentation and the compiler. – Clinton Winant Jun 02 '17 at 16:12
@ClintonWinant The reshape is a way how to view a (10) array as (1,10) or (10,1) array. The explanation remains the same. – Vladimir F Героям слава Jun 02 '17 at 16:37
@ClintonWinant You are correct -- the gfortran documentation is incorrect. Here is the relevant line from the Fortran Standard: "MATRIX A and MATRIX B shall not both have rank one." http://www.j3-fortran.org/doc/year/10/10-007.pdf – Jack Jun 04 '17 at 17:10
@Jack Thanks, just to close the loop, gfortran ONLY accepts both arrays having rank 2, otherwise the compiler returns: `Error: Incompatible ranks 2 and 1 in assignment at (1)` even though the standard says otherwise. I believe we have exhausted the subject...... – Clinton Winant Jun 05 '17 at 08:57
Actually, no. With `A(10)` and `B(10)`, gfortran will take BOTH `C = matmul( reshape(A,(/1,10/)), B )` and `C = matmul( A, reshape(B,(/10,1/)) )`. However, in both cases, C has to be a single-element array of rank 1: `C(1)` – Jack Jun 05 '17 at 11:41
@ClintonWinant The new error message is about the assignment, not about `matmul`. If you multiply matrix by a vector, you will get a vector. You cannot assign a vector to a matrix. – Vladimir F Героям слава Jun 05 '17 at 15:37

Vladimir F Героям слава · Answer 2 · 2017-06-02T20:19:07.463

1

You must do this for a tensor product of two vectors

Program scratch
  integer, parameter :: dp = kind(1.d0)
  real(dp) :: A(10,1)=reshape((/0,1,2,3,4,5,6,7,8,9/), (/ 10, 1 /))
  real(dp) :: B(1,10)=reshape((/0,1,2,3,4,5,6,7,8,9/), (/ 1, 10 /))
  real(dp) :: C(10,10)
  print *,rank(A),rank(B)
  C=matmul(A,B)
  print *, C
End Program scratch

If you do

   A(1,10)
   B(10,1)

you will get a scalar product. With just two 1D arrays it is not clear which of the two products you want (although for a dot product there is a special function available).

A or B can be a 1D array when you are multiplying a matrix by a vector.

edited Jun 02 '17 at 20:19

answered Jun 02 '17 at 15:07

Vladimir F Героям слава

57,977
4
76
119

so the point is that in spite of what the manual says about matrix inputs of rank 1 being admissible, the compiler is smarter, knows that this is ambiguous, and insists on rank 2. Because the dot_product intrinsic exists, I thought matmul would default to the dyadic product, wrong... Thank you – Clinton Winant Jun 02 '17 at 16:04
1

Rank 1 is aloowed, but not both arguments. – Vladimir F Героям слава Jun 02 '17 at 16:38
1

This line `real(dp) :: A(10,1)=(/0,1,2,3,4,5,6,7,8,9/)` gives me the error **Error: Incompatible ranks 2 and 1 in assignment at (1)** in gfortran, and **PGF90-S-0155-Shape of initializer does not match shape of a (scratch.f90: 4)** in pgfortran. Am I doing something wrong here? – Jack Jun 02 '17 at 17:14
No, I think I have an error there, but my compiler allowed that for some reason. I did test it. – Vladimir F Героям слава Jun 02 '17 at 20:17
What compiler is that? – Jack Jun 02 '17 at 21:25
Intel Fortran . – Vladimir F Героям слава Jun 02 '17 at 21:35

Fortran array rank for matmul intrinsic

2 Answers2

Linked

Related