Error in mul!(c,a,b) if size(b) = (n,1)

General Usage
21

As a Matlab person, I can assure you that it is a blessed relief that Julia distinguishes scalars, vectors and matrices.

The way Matlab works in this regard is a terrible headache, with zero advantages.

3 Likes
22

The easiest way to understand it is to simply count the dimensions. (n, 1) are 1, 2 dimensions. (n,) is 1 dimension. The fact that the length of the last dimension is 1 is incidental. It could even be (n, 0) and it’s still 2D.

2 Likes
23

Actually, the Fortran example works in Julia as well:

julia> x = collect([1 1 1;]')
3×1 Array{Int64,2}:
 1
 1
 1

julia> mysum(x) = sum(x[i] for i in 1:length(x));

julia> mysum(x)
3

You can access a Julia multi-dimensional array by a single (column-major) index, but that’s not the same as their types being the same. I suspect gfortran’s MATMUL will also choke on two vector inputs: “MATRIX_A and MATRIX_B shall not both be rank one arrays”.

2 Likes
24

Well, that is why I think that the error thrown by those functions is an imposition of the type verification only. If a method accepted the input unidimensional matrix and feed the same calculations with the correct pointers, the result must be the same.

25

Actually it works:

program main
  double precision :: x(3), y(3,1)
  do i = 1, 3
    x(i) = 1.
    y(i,1) = 1.
  end do
  write(*,*) matmul(x,y)
end program main

$ gfortran matmul.f90 -o matmul
$ ./matmul
   3.0000000000000000

Edit: it does choke if x is declared as x(3,1).

26

I sort of agree with you, because for this case the memory layout is what you are thinking. But as you found, Fortran’s MATMUL can also choke from its type verification. (BTW your example doesn’t show two vectors x(3) and y(3) which the documentation says should fail.)

The admonition against thinking in terms of pointers and memory layout deserves some explanation. The type system allows for more abstraction than a Fortran real array, and allows things like sparse arrays, or LU decompositions, or adjoint arrays to be used with the same common syntax like *, and things will usually just work. The abstraction is incredibly powerful, and enables automatic differentiation to propagate through packages that were never written with it in mind. Enabled by the type system. Similarly, curve_fit should just work if fed abstractions of 2d arrays with entirely different memory layout. (But alas, not 1d arrays…)

1 Like
27

I think my expectation on this also comes from the fact that a naive algorithm would work in these cases:

julia> function matmul(a,b)
         na = size(a,1); ma = size(a,2)
         nb = size(b,1); mb = size(b,2)
         if ( ma != nb ); error(" ma != nb "); end
         c = zeros(na,mb)
         for i in 1:na
           for j in 1:mb
             for k in 1:ma
               c[i,j] = c[i,j] + a[i,k]*b[k,j]
             end
           end
         end
         c
       end

julia> a = ones(3); b = ones(3,1) ;

julia> matmul(transpose(a),b)
1×1 Array{Float64,2}:
 3.0

julia> a = ones(3,1); b = ones(3) ;

julia> matmul(transpose(a),b)
1×1 Array{Float64,2}:
 3.0

julia> a = ones(3,1); b = ones(3,1);

julia> matmul(transpose(a),b)
1×1 Array{Float64,2}:
 3.0

Edit: Yes, fortran complains at compiling time if both are vectors:

Error: ‘matrix_b’ argument of ‘matmul’ intrinsic at (1) must be of rank 2
28

…a naive algorithm would work in these cases

Yes, but the issue is not mul! per se. As @mcabbott explains, curve_fit doesn’t anticipate matrix input, with a potentially very easy fix.

In fact, mul! does handle vectors or matrices, it just doesn’t like mismatched types (i.e. Matlab ambiguity). Also, the documentation calls mul! a low-level function, which signals that higher level functions should take responsibility for massaging its various arguments before calling it. mul! is specific to get its high performance, whereas your example or even a plain * would work more generally at the cost of new allocation. So I would look to curve_fit (or perhaps lmfit) to handle edge cases if it wants to, or have clearer documentation.

1 Like
29

But in Julia, it isn’t. Also, it is not something that needs to be explicitly documented, as the manual gives you no reason to assume that it is.

1 Like