Multidimensional arrays

C++ has no multidimensional array type (MDA) in its standard library. However, MDAs are crucial for the developmemt of scientific applications. One of the reasons for the continuing success of languages like Fortran or Python is their excellent support for MDAsfootnote{For Python arrays are introduced by the texttt{numpy} package.}. The lack of an MDA type in C++ was indeed the spark that initiated the developement of libpninexus. Before discussing libpninexus s array facilities some terminology should be defined:

element type (ET)

referes to the data type of the individual elements stored in an MDA. For MDAs this will typically be a numeric type like an integer or a floating point number.

rank \(r\)

denotes the number of dimensions of an MDA

shape \(\mathbf{s}\)

is a vector of dimension \(r\) whose elements are the number of elements along each dimension. The elements of \(\mathbf{s}\) are denoted as \(s_i\) with \(i=0,\dots,r-1\)

The hart of libpninexus s MDA support is the mdarray template. mdarray is extremely powerfull. Thus, using mdarray directly to define array types is not for the faint harted. To simplify the usage of multidimensional arrays the library provides three templates derived form mdarray which are easy to use.

template alias

description

static_array<T,size_t …N>

a static arrays whose shape, rank, and element type are fixed at compile time.

fixed_dim_array<T,R>

element type and rank are fixed at compile time but the shape can be changed at runtime.

dynamic_array<T>

a fully dynamic array type where only the element type must be known at compile time. The rank as well as the shape can be altered at runtime

These types are all defined in pni/arrays.hpp. In addition to this two basic templates there are several utility classes and templates like

template

functionality

arrayview

this template provides a view on a particular part of an array

arrayerasure

type erasure for an array

All array types derived from mdarray provide the following features

  • unary and binary arithmetics if the element type is a numeric type.

  • slicing to extract only a part of a large array

  • simple access to data elements using variadic operators.

  • all array types are full STL compliant containers and thus can be used along with STL algorithms.

Array construction and inquery

Constructing arrays is rather simple by means of the create() function provided by the array templates. As shown in array_create.cpp

see examples/core/array_create.cpp
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <vector>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

//some usefull type definitions
typedef dynamic_array<float64> darray_type;

int main(int ,char **)
{
    //construction from shape
    auto a1 = darray_type::create(shape_t{1024,2048});

    //construction from shape and buffer
    auto a2 = darray_type::create(shape_t{1024,2048},
                                  darray_type::storage_type(1024*2048));

    //construction from initializer lists
    auto a3 = darray_type::create({5},{1,2,3,4,5});

    return 0;
}

In line 35 a concrete array type is defined from the dynamic_array utility template. The create() function comes in three flavors as shown in the previous example

line

description

40

Only the shape is passed as an argument. The storage container is allocated automatically by the constructor.

43

The storage container with initial data is passed as a second argument to create() along with the shape of the array.

47

shape as well as initial data are provided via initializer lists.

After an array has been created we may want to retrieve some of its basic properties. In the next example we do exactly this

see examples/core/array_inquery.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <vector>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

//some usefull type definitions
typedef dynamic_array<float64> darray_type;

template<typename ArrayT>
void show_info(const ArrayT &a)
{
    std::cout<<"Data type: "<<type_id(a)<<std::endl;
    std::cout<<"Rank     : "<<a.rank()<<std::endl;
    std::cout<<"Shape    : (";
    auto s = a.template shape<shape_t>();
    for(auto n: s) std::cout<<" "<<n<<" ";
    std::cout<<")"<<std::endl;
    std::cout<<"Size     : "<<a.size()<<std::endl;
}

int main(int ,char **)
{
    auto a1 = darray_type::create(shape_t{1024,2048});

    show_info(a1);

    return 0;
}

The important part is the implementation of the show_info() function template starting at line 37. The function template type_id() is used in line 40 to retrieve the type ID of the arrays element type. rank() in line 41 returns the number of dimension and size() in line 46 the total number of elements stored in the array. The shape() template function in line 43 returns the number of elements along each dimension stored in a user provided container type.

The content of arrays can be copied to and from containers using the standard std::copy() template function from the STL. In addition a version of the assignment operator is provided which allows assignment of values from an initializer list. This is particularly useful for static arrays which basically do not require construction.

using static_array_type = ....;

static_array_type a;
a = {1,2,3,4,5};

Linear access to data

As already mentioned in the first section of this chapter, the array types provided by libpninexus are fully STL compliant containers. They provided all the iterators required by the STL. Before we have a look on STL lets first investigate how to simply access data elements in an array

see examples/core/array_linear_access.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef uint16                          channel_type;
typedef fixed_dim_array<channel_type,1> mca_type;

int main(int ,char **)
{
    auto mca = mca_type::create(shape_t{128}); 

    //initialize
    std::random_device rdev;
    std::mt19937 generator(rdev());
    std::uniform_int_distribution<channel_type> dist(0,65000);

    //generate data
    for(auto &channel: mca)  channel = dist(generator);

    //subtract some number 
    for(size_t i=0;i<mca.size();++i) if(mca[i]>=10) mca[i] -= 10;

    //set the first and last element to 0 
    mca.front() = 0;
    mca.back()  = 0;
    
    //output data
    for(auto channel: mca) std::cout<<channel<<std::endl;

    return 0;
}

For all array types the new C++ for-each construction can be used as shown in lines 47 and 57. Unchecked access (no index bounds are checked) is provided via the [] operator as demonstrated in line 50 (use at() for checked access). Some of the operations in this example can be done much more efficient with STL algorithms as demonstrated in the next example

see examples/core/array_stl.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef uint16                        pixel_type;
typedef fixed_dim_array<pixel_type,2> image_type;

int main(int ,char **)
{
    auto image = image_type::create(shape_t{1024,512});

    std::fill(image.begin(),image.end(),0); //use STL std::fill to initialize

    std::random_device rdev;
    std::mt19937 generator(rdev());
    std::uniform_int_distribution<pixel_type> dist(0,65000);

    std::generate(image.begin(),image.end(),
                  [&generator,&dist](){ return dist(generator); });

    //get min and max
    std::cout<<"Min: "<<*std::min_element(image.begin(),image.end())<<std::endl;
    std::cout<<"Max: "<<*std::max_element(image.begin(),image.end())<<std::endl;
    std::cout<<"Sum: ";
    std::cout<<std::accumulate(image.begin(),image.end(),pixel_type(0));
    std::cout<<std::endl;
                               
    return 0;
}

In line 42 std::fill() is used to initialize the array to 0 and std::generate() in line 48 fills it with random numbers using a lambda expression . The rest of the example should be trivial (if not, please lookup a good C++ STL reference).

Multidimensional access

Though being an important feature, linear access to multidimensional arrays is not always useful. In particular the last example where we pretended to work on image data implementing algorithms would be rather tedious if we would have had only linear access. It is natural for such objects to think in pixel coordinates (i,j) rather than the linear offset in memory. libpninexus provides easy multidimensional access to the data stored in an array. The next example shows how to use this feature to work only on a small region of the image data as defined in the last example

see examples/core/array_multiindex.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef uint16                        pixel_type;
typedef std::array<size_t,2>          index_type;
typedef fixed_dim_array<pixel_type,2> image_type;

int main(int ,char **)
{
    auto image = image_type::create(shape_t{1024,512});

    std::fill(image.begin(),image.end(),0); //use STL std::fill to initialize

    std::random_device rdev;
    std::mt19937 generator(rdev());
    std::uniform_int_distribution<pixel_type> dist(0,65000);

    std::generate(image.begin(),image.end(),
                  [&generator,&dist](){ return dist(generator); });

    
    size_t zero_count = 0;
    size_t max_count  = 0;
    for(size_t i=512;i<934;++i)
    {
        for(size_t j=128;j<414;++j)
        {
            if(image(i,j) == 0) zero_count++;
            if(image(index_type{{i,j}}) >= 10000) max_count++;
                
        }
    }

    std::cout<<"Found 0 in "<<zero_count<<" pixels!"<<std::endl;
    std::cout<<"Found max in "<<max_count<<" pixels!"<<std::endl;
                               
    return 0;
}

The interesting part here are lines 59 and 60. You can pass the multidimensional indexes either as a variadic argument list to the () operator of the array type (as in line 59) or you can use a container like in line 60. The former approach might look a bit more familiar, however, in some cases when decisions have to made at runtime the container approach might fits better. However, passing containers reduces access performance approximately by a factor of 2. Thus, as a rule of thumb you should always use the variadic form when you know the number of dimensions the array has and containers only in those cases where this information is only available at runtime.

Array views and slicing

In the previous example multiindex access was used to do work on only a small part of the image data. libpninexus provides view types for arrays which would make these operations easier. Views are created by passing instances of slice to the () operator of an array type. Slices in libpninexus work pretty much the same as in python. Lets have a look on the following example

see examples/core/array_view.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef uint16                        pixel_type;
typedef std::array<size_t,2>          index_type;
typedef fixed_dim_array<pixel_type,2> image_type;

int main(int ,char **)
{
    auto image = image_type::create(shape_t{1024,512});

    std::fill(image.begin(),image.end(),0); //use STL std::fill to initialize

    std::random_device rdev;
    std::mt19937 generator(rdev());
    std::uniform_int_distribution<pixel_type> dist(0,65000);

    std::generate(image.begin(),image.end(),
                  [&generator,&dist](){ return dist(generator); });

    
    auto roi = image(slice(512,934),slice(128,414));
    auto zero_count = std::count_if(roi.begin(),roi.end(),
                                    [](pixel_type &p){return p==0;});
    auto max_count  = std::count_if(roi.begin(),roi.end(),
                                    [](pixel_type &p){return p>= 10000; });

    std::cout<<"Found 0 in "<<zero_count<<" pixels!"<<std::endl;
    std::cout<<"Found max in "<<max_count<<" pixels!"<<std::endl;
                               
    return 0;
}

The view is created in line 53 where the slices are passed instead of integer indices to the () operator. A slice selects an entire index range along a dimension. The first argument to the slice constructor is the starting index and the last the stop index of the range. The stop index is not included (just as it is the case with Python slices). If the () operator of an array is called with any of its arguments being a slice a view object is returned instead of a single value or reference to a single value. View objects are pretty much like arrays themselves. However, they do not hold data by themselves but only a reference to the original array. Like arrays they are fully STL compliant containers and thus can be used with STL algorithms as shown in lines 54 and 56.

View types can be copied and moved and thus can be stored in STL containers as shown in the next example

see examples/core/array_view_container.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef uint16                        pixel_type;
typedef std::array<size_t,2>          index_type;
typedef fixed_dim_array<pixel_type,2> image_type;
typedef image_type::view_type         roi_type;
typedef std::vector<roi_type>         roi_vector;

int main(int ,char **)
{
    auto image = image_type::create(shape_t{1024,512});

    std::fill(image.begin(),image.end(),0); //use STL std::fill to initialize

    std::random_device rdev;
    std::mt19937 generator(rdev());
    std::uniform_int_distribution<pixel_type> dist(0,65000);

    std::generate(image.begin(),image.end(),
                  [&generator,&dist](){ return dist(generator); });

    roi_vector rois; 
    rois.push_back(image(slice(512,934),slice(128,414)));
    rois.push_back(image(slice(0,128),slice(4,100)));
    rois.push_back(image(200,slice(450,512)));

    for(auto roi: rois)
    {
        auto zero_count = std::count_if(roi.begin(),roi.end(),
                                        [](pixel_type &p){return p==0;});
        auto max_count  = std::count_if(roi.begin(),roi.end(),
                                        [](pixel_type &p){return p>= 10000; });

        std::cout<<std::endl;
        std::cout<<"Found 0 in "<<zero_count<<" pixels!"<<std::endl;
        std::cout<<"Found max in "<<max_count<<" pixels!"<<std::endl;
        std::cout<<std::endl;
    }
                               
    return 0;
}

Here we apply the algorithms from the previous example not to a single but to several selections in the image. As shown in lines 55 to 57 we can safely store views in a container and later iterate over it.

In general views make algorithm development much easier as we have to develop algorithms only for entire arrays. If it should be applied to only a part of an array we can use a view and pass it to the algorithm. As views expose the same interface as an array the algorithm should work on views too.

Arithmetic expressions

Array and view types fully support the common arithmetic operators +, *, /, and - in their binary and unary forms. The binary versions are implemented as expression templates avoiding the allocation of unnecessary temporary and giving the compiler more possibilities to optimize the code. Views, arrays and scalars can be mixed within all arithmetic expressions. There is nothing magical with expression templates as they work entirely transparent to the user. Just use the arithmetic expressions as you are used to

see examples/core/array_arithmetic1.cpp
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef float64                        number_type;
typedef fixed_dim_array<number_type,2> image_type;

int main(int ,char **)
{
    shape_t s{1024,512};
    auto image      = image_type::create(s);
    auto background = image_type::create(s);
    number_type exp_time = 1.234;
    number_type current  = 98.3445;

    //compute the corrected image
    image = (image-background)/exp_time/current;
                               
    return 0;
}

The important line here is 47 where arrays and scalars are mixed in an arithmetic expression. One can also mix arrays, selections, and scalars as the next examples shows

see examples/core/array_arithmetic2.cpp for full code
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

typedef float64                        number_type;
typedef fixed_dim_array<number_type,3> stack_type;
typedef fixed_dim_array<number_type,2> image_type;
typedef fixed_dim_array<number_type,1> data_type;

int main(int ,char **)
{
    shape_t frame_shape{1024,512};
    shape_t data_shape{100};
    shape_t stack_shape{100,1024,512};
    auto image_stack = stack_type::create(stack_shape);
    auto background  = image_type::create(frame_shape);
    auto exp_time    = data_type::create(data_shape);
    auto current     = data_type::create(data_shape);

    for(size_t i = 0;i<data_shape[0];++i)
    {
        std::cout<<"Processing frame "<<i<<" ..."<<std::endl;
        auto curr_frame = image_stack(i,slice(0,1024),slice(0,512));
        curr_frame = (curr_frame - background)/exp_time[i]/current[i];
    }
                               
    return 0;
}

In line 53 a single image frame is selected from a stack of images and used in line 54 in an arithmetic expression. In fact, what we are doing here is, we are writing the corrected data back on the stack since curr_frame is just a view on the particular image in the stack.

Example: matrix-vector and matrix-matrix multiplication

In the last example matrix vector multiplications are treated. The full code can be viewed in array_arithmetic3.cpp in the source distribution. But lets first start with the header

see examples/core/array_arithmetic3.cpp for full code
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include <iostream>
#include <algorithm>
#include <random>
#include <pni/types.hpp>
#include <pni/arrays.hpp>

using namespace pni;

// define the matrix and  vector types
template<typename ElementT,size_t TDimN> using matrix_temp = static_array<ElementT,TDimN,TDimN>;
template<typename ElementT,size_t TDimN> using vector_temp = static_array<ElementT,TDimN>;

//-----------------------------------------------------------------------------
// print a matrix
template<typename ElementT,size_t TDimN>
std::ostream &operator<<(std::ostream &o,const matrix_temp<ElementT,TDimN> &m)
{
    for(size_t i = 0;i<3;++i)
    {
        o<<"| ";
        for(size_t j=0;j<3;++j) o<<m(i,j)<<" ";
        o<<"|"<<std::endl;
    }

    return o;
}

//-----------------------------------------------------------------------------
// print a vector
template<typename ElementT,size_t TDimN>
std::ostream &operator<<(std::ostream &o,const vector_temp<ElementT,TDimN> &v)
{
    for(auto x: v) o<<"| "<<x<<" |"<<std::endl;
    return o;
}

//-----------------------------------------------------------------------------
// matrix-vector multiplication
template<typename ElementT,size_t TDimN>
vector_temp<ElementT,TDimN> mv_mult(const matrix_temp<ElementT,TDimN> &m,const vector_temp<ElementT,TDimN> &v)
{
    vector_temp<ElementT,TDimN> result;

    size_t i = 0;
    for(auto &r: result)
    {
        const auto row = m(i++,slice(0,TDimN));
        r = std::inner_product(v.begin(),v.end(),row.begin(),ElementT(0));
    }
    return result;
}

//-----------------------------------------------------------------------------
// vector-matrix multiplication
template<typename ElementT,size_t TDimN>
vector_temp<ElementT,TDimN> mv_mult(const vector_temp<ElementT,TDimN> &v,const matrix_temp<ElementT,TDimN> &m)
{
    vector_temp<ElementT,TDimN> result;

    size_t i = 0;
    for(auto &r: result)
    {
        const auto col = m(slice(0,TDimN),i++);
        r = std::inner_product(col.begin(),col.end(),v.begin(),ElementT(0));
    }
    return result;
}

//-----------------------------------------------------------------------------
// matrix-matrix multiplication
template<typename ElementT,size_t TDimN>
matrix_temp<ElementT,TDimN> mv_mult(const matrix_temp<ElementT,TDimN> &m1,const matrix_temp<ElementT,TDimN> &m2)
{
    matrix_temp<ElementT,TDimN> result;

    for(size_t i=0;i<TDimN;++i)
    {
        for(size_t j=0;j<TDimN;++j)
        {
            const auto row = m1(i,slice(0,TDimN));
            const auto col = m2(slice(0,TDimN),j);
            result(i,j) = std::inner_product(row.begin(),row.end(),
                                             col.begin(),ElementT(0));
        }
    }
    return result;
}

//-----------------------------------------------------------------------------
// define some local types
typedef float64                    number_type;
typedef vector_temp<number_type,3> vector_type;
typedef matrix_temp<number_type,3> matrix_type;


int main(int ,char **)
{
    vector_type v;
    matrix_type m1,m2; 
    m1 = {1,2,3,4,5,6,7,8,9};
    m2 = {9,8,7,6,5,4,3,2,1};
    v  = {1,2,3};

    std::cout<<"m1 = "<<std::endl<<m1<<std::endl;
    std::cout<<"m2 = "<<std::endl<<m2<<std::endl;
    std::cout<<"v  = "<<std::endl<<v<<std::endl;
    std::cout<<"m1.v = "<<std::endl<<mv_mult(m1,v)<<std::endl;
    std::cout<<"v.m1 = "<<std::endl<<mv_mult(v,m1)<<std::endl;
    std::cout<<"m1.m2 = "<<std::endl<<mv_mult(m1,m2)<<std::endl;
                               
    return 0;
}

Besides including all required header files matrix and vector templates are defined in lines 37 and 38 using the new C++11 template aliasing.

Matrix vector multiplication

The implementation of the matrix vector multiplication is shown in the next block.

\[\mathbf{r} = A\mathbf{v} \mbox{ or } r_j = A_{j,i}v_i\]

with \(\mathbf{A}\) denoting a \(N\times N\) matrix and \(\mathbf{r}\) and \(\mathbf{v}\) are \(N\)-dimensional vectors. In all formulas Einsteins sum convention is used.

66
67
68
69
70
71
72
73
74
75
76
77
78
template<typename ElementT,size_t TDimN>
vector_temp<ElementT,TDimN> mv_mult(const matrix_temp<ElementT,TDimN> &m,const vector_temp<ElementT,TDimN> &v)
{
    vector_temp<ElementT,TDimN> result;

    size_t i = 0;
    for(auto &r: result)
    {
        const auto row = m(i++,slice(0,TDimN));
        r = std::inner_product(v.begin(),v.end(),row.begin(),ElementT(0));
    }
    return result;
}

In line 74 we select the \(i\)-th row of the matrix and compute the inner product of the row vector and the input vector in line 75.

Vector matrix multiplication

The vector matrix multiplication

\[\mathbf{r} = \mathbf{v}A\mbox{ or } r_i = v_j A_{j,i}\]

is computed analogously

82
83
84
85
86
87
88
89
90
91
92
93
94
template<typename ElementT,size_t TDimN>
vector_temp<ElementT,TDimN> mv_mult(const vector_temp<ElementT,TDimN> &v,const matrix_temp<ElementT,TDimN> &m)
{
    vector_temp<ElementT,TDimN> result;

    size_t i = 0;
    for(auto &r: result)
    {
        const auto col = m(slice(0,TDimN),i++);
        r = std::inner_product(col.begin(),col.end(),v.begin(),ElementT(0));
    }
    return result;
}

despite the fact that we are choosing the appropriate column instead of a row in line 90.

Matrix matrix multiplication

Finally we need an implementation for the matrix - matrix multiplication

\[C = AB \mbox{ or } C_{i,j} = A_{i,k}B_{k,j}\]
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
template<typename ElementT,size_t TDimN>
matrix_temp<ElementT,TDimN> mv_mult(const matrix_temp<ElementT,TDimN> &m1,const matrix_temp<ElementT,TDimN> &m2)
{
    matrix_temp<ElementT,TDimN> result;

    for(size_t i=0;i<TDimN;++i)
    {
        for(size_t j=0;j<TDimN;++j)
        {
            const auto row = m1(i,slice(0,TDimN));
            const auto col = m2(slice(0,TDimN),j);
            result(i,j) = std::inner_product(row.begin(),row.end(),
                                             col.begin(),ElementT(0));
        }
    }
    return result;
}

The rows and columns are selected in lines 107 and 108 respectively. Line 109 finally computes the inner product of the row and column vector.

Putting it all together: the main function

Finally the main program shows a simple application of these template functions.

123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
int main(int ,char **)
{
    vector_type v;
    matrix_type m1,m2; 
    m1 = {1,2,3,4,5,6,7,8,9};
    m2 = {9,8,7,6,5,4,3,2,1};
    v  = {1,2,3};

    std::cout<<"m1 = "<<std::endl<<m1<<std::endl;
    std::cout<<"m2 = "<<std::endl<<m2<<std::endl;
    std::cout<<"v  = "<<std::endl<<v<<std::endl;
    std::cout<<"m1.v = "<<std::endl<<mv_mult(m1,v)<<std::endl;
    std::cout<<"v.m1 = "<<std::endl<<mv_mult(v,m1)<<std::endl;
    std::cout<<"m1.m2 = "<<std::endl<<mv_mult(m1,m2)<<std::endl;
                               
    return 0;
}

It is important to understand that the appropriate function is determined by the types of the arguments (vector or matrix). This is a rather nice example of how to use the typing system of C++ to add meaning to objects.

The output of the program is

>./array_arithmetic3
m1 =
| 1 2 3 |
| 4 5 6 |
| 7 8 9 |

m2 =
| 9 8 7 |
| 6 5 4 |
| 3 2 1 |

v  =
| 1 |
| 2 |
| 3 |

m1.v =
| 14 |
| 32 |
| 50 |

v.m1 =
| 30 |
| 36 |
| 42 |

m1.m2 =
| 30 24 18 |
| 84 69 54 |
| 138 114 90 |