Non intuitive MATMUL behaviour

[GerritGroot]

Hi,

I'v got some strange unexpected behaviour of MATMUL in F90, which seems to be compiler independent, so I think I didn't get the point of how MATMUL works

I got:

PROGRAM TestMatMul
IMPLICIT NONE
REAL, DIMENSION(4,4) :: A
REAL, DIMENSION(4) :: x,y

A=RESHAPE(SOURCE=(/1.0000000,0.0000000,0.0000000,-30.000000, &
                   0.0000000,1.0000000,0.0000000,  0.0000000, &
                   0.0000000,0.0000000,1.0000000,  0.0000000, &
                   0.0000000,0.0000000,0.0000000,  1.0000000/),SHAPE=(/4,4/))
x =RESHAPE(SOURCE=(/30,0,0,1/),SHAPE=(/4/))

y=MATMUL(A,x)
WRITE(*,*)y

END PROGRAM TestMatMul

Gives as output

      30.000000    0.000000E+00    0.000000E+00     -899.000000

While I expected it to give the same result as Octave. The following code in Octave or Matlab

x=[30
   0
   0
   1];
M=[1.0000000       0.0000000       0.0000000      -30.000000
   0.0000000       1.0000000       0.0000000       0.0000000
   0.0000000       0.0000000       1.0000000       0.0000000
   0.0000000       0.0000000       0.0000000       1.0000000];
y=M*x;

Gives the result:

ans = 
       0
       0
       0
       1

How canI get it right in MATMUL?

 

[FJacq]

The problem is not MATMUL but how terms in a matrix are organized.

FORTRAN orders them by column whereas C and apparently OCTAVE organize them by row.

Try the following pogram :

PROGRAM TestMatMul
IMPLICIT NONE
REAL, DIMENSION(4,4) :: A
REAL, DIMENSION(4) :: x,y

a(1,:)=(/1.0000000,0.0000000,0.0000000,-30.000000/)
a(2,:)=(/0.0000000,1.0000000,0.0000000, 0.0000000/)
a(3,:)=(/0.0000000,0.0000000,1.0000000, 0.0000000/)
a(4,:)=(/0.0000000,0.0000000,0.0000000, 1.0000000/)

x=(/30.,0.,0.,1./)

y=MATMUL(A,x)
WRITE(*,*)y

END PROGRAM TestMatMul

François Jacq

 

[GerritGroot]

Thanks, this one does the trick. Ye have to admit that it is not intuitive indeed.

Now that you say so, I vaguely remember having read this for fortran 77, but I certainly didn't expect to find this legacy in a standard Fortran 90 routine.

This choice had something to do with memory efficiency, right?

Anyway, I don't like to adapt myself to the computer. Computers are there to help us and adapt to us So I'll make my own MATMULT where I can fill in the matrix the way they are being written on paper.

Gerrit

 

[FJacq]

FORTRAN has never changed : the present standard (FORTRAN-2008) is still compatible with oldest FORTRAN codes (with few exceptions) : backward compatibility over more than 40 years.

Why doing that choice ? Really I don't know. And we cannot ask John Bakus anymore (the IBM's guy who inverted FORTRAN in 1957 and who recently died).

I think that MATLAB follows the same rule as FORTRAN.

In practice today, this rule has no particular consequence except :
- when one wants to fill a matrix from a vector (with RESHAPE or DATA).
- when one takes care about performance (the loop on the row index must be the most internal one).

François Jacq

 

[GerritGroot]

Yes, you're right, but I prefer to have my data in human readable format, so I prefer to program them also in human readable order, to prevent mistakes (if not one may switch to assembly )

Anyway, it's not such a big deal. I didn't know about the existense of the TRANSPOSE function for example, so I didn't even have to make a routine to get MATMUL working in a human way, I just took y=MATMUL(TRANSPOSE(M),x)

 

[FJacq]

Yes, you're right, but I prefer to have my data in human readable format, so I prefer to program them also in human readable order, to prevent mistakes (if not one may switch to assembly smile )

Look into the modified program I proposed : it is in human readable format.

I repeat that a potential trouble only occurs when one wants to construct a matrix from a list of values (vector). My programming avoids that difficulty.

François Jacq

 

[xwb]

It is a purely mathematical concept. If you have an xy graph and drew it out, the x would normally be on the horizontal and y on the vertical. So reading horizontally, you'd get

x=1,y=1  x=2,y=1 x=3,y=1
x=1,y=2  x=2,y=2 x=3,y=2
x=1,y=3  x=2,y=3 x=3,y=3

Matrices were thought of as rows and columns: not columns and rows so naturally the rows come first. The result is a yx way of thinking instead of an xy way of thinking.

Also with the advent of other languages, there is the ownership concept where the index on the left owns the ones on the right. In Fortran, the ownership concept was only brought in in F77: about 24 years too late. Since it had to be backward compatible, the storage mechanism never changed.