1

My question is similar to this one Multiply a 3D matrix with a 2D matrix. However, I'm coding in Fortran.

Say, if I have a RxSxT matrix A and an SxU matrix B, where R,S,T,U are integers, and I want to multiply A(:,:,0) with B. How can I do this with matmul? When I do something like

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

The compiler (gfortran) gives:

    Warning:Array reference at (1) is out of bounds (0 < 1) in dimension 3 
    f951: internal compiler error: Segmentation fault

Is there a way around this? Thanks.

EDIT: I should add that I'm actually transposing the second matrix. Say, A a RxSxT matrix and B a UxS matrix. Then

 C(:,:,0) = matmul(B,transpose(A(:,:,0))

That transpose might be part of the problem. Does it convert A(i,j,k) to A(k,i,j)?

Community
  • 1
  • 1
Samuel Tan
  • 1,700
  • 5
  • 22
  • 35

2 Answers2

8

Remember that in Fortran your array indices start from 1 by default. So unless you have specified your array A to have a non-default lower bound on the 3rd dimension, gfortran is entirely correct in pointing out your error.

Of course, an internal compiler error is always a compiler bug; unless you have some ancient version of gfortran please file a bug at http://gcc.gnu.org/bugzilla

janneb
  • 36,249
  • 2
  • 81
  • 97
  • So if I have A(:,:,1) instead it should be fine? It seems like that problem is resolved, but I still get the seg fault. My version is 4.6, which I think is the latest stable (non-development) version. – Samuel Tan Dec 15 '11 at 01:16
2

transpose (A(:,:,0)) should interchange the indices A(i,j,0) to A(j,i,0). A(:,:,0) is a rank two matrix.

The compiler should never crash, whether or not the input source code is correct. Are you using the latest version of gfortran? You could report this "internal compiler error: Segmentation fault" to the gfortran development team: http://gcc.gnu.org/wiki/GFortran#bugs

M. S. B.
  • 28,968
  • 2
  • 46
  • 73
  • That's interesting. So transpose(A(i,j,k)) gives A(j,i,k)? What if I want to go from A(i,j,k) to A(k,j,i)? Or to A(i,k,j)? What command should I use? – Samuel Tan Dec 15 '11 at 01:17