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

**URL:** https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999
**Category:** General Usage
**Created:** [October 8, 2020, 8:02pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999 "2020-10-08T20:02:41Z")
**Posts on this page:** 9
**Page:** 2

<div class="post-metadata">

### Author: ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [October 9, 2020, 8:07pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/21 "2020-10-09T20:07:38Z")

</div>

> [@apo383](#):
>
> It’s a pain point for Matlab people, but worth it to survive.

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.

---

<div class="post-metadata">

### Author: ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [October 9, 2020, 8:15pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/22 "2020-10-09T20:15:22Z")

</div>

> [@lmiq](#):
>
> I think it is somewhat hard to document the fact that an array of sizes `(n,1)` is not a 1D array.

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.

---

<div class="post-metadata">

### Author: ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [October 9, 2020, 8:20pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/23 "2020-10-09T20:20:17Z")

</div>

Actually, the Fortran example works in Julia as well:

```julia
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](https://gcc.gnu.org/onlinedocs/gfortran/MATMUL.html) will also choke on two vector inputs: “MATRIX\_A and MATRIX\_B shall not both be rank one arrays”.

---

<div class="post-metadata">

### Author: ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [October 9, 2020, 10:11pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/24 "2020-10-09T22:11:17Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [October 10, 2020, 1:12am UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/25 "2020-10-10T01:12:14Z")

</div>

> [@apo383](#):
>
> I suspect [gfortran’s MATMUL](https://gcc.gnu.org/onlinedocs/gfortran/MATMUL.html) will also choke on two vector input

Actually it works:

```fortran
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)`.

---

<div class="post-metadata">

### Author: ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [October 10, 2020, 1:37am UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/26 "2020-10-10T01:37:22Z")

</div>

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…)

---

<div class="post-metadata">

### Author: ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [October 10, 2020, 1:49am UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/27 "2020-10-10T01:49:30Z")

</div>

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

```julia
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:

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

```

---

<div class="post-metadata">

### Author: ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)
#### Post date: [October 10, 2020, 4:10am UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/28 "2020-10-10T04:10:18Z")

</div>

> …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](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/11).

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.

---

<div class="post-metadata">

### Author: ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [October 10, 2020, 2:50pm UTC](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999/29 "2020-10-10T14:50:09Z")

</div>

> [@lmiq](#):
>
> it is somewhat hard to document the fact that an array of sizes `(n,1)` is not a 1D array.

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.

[Previous page](https://discourse.julialang.org/t/error-in-mul-c-a-b-if-size-b-n-1/47999.md?page=1)
