# Implement matrix operations

## In this learning path

- [Introduction](https://learn.arm.com/learning-paths/cross-platform/matrix/)
- [Laying the foundations](https://learn.arm.com/learning-paths/cross-platform/matrix/1-foundations/)
- [Test the library](https://learn.arm.com/learning-paths/cross-platform/matrix/2-testing/)
- [Start coding](https://learn.arm.com/learning-paths/cross-platform/matrix/3-code-1/)
- [Implement matrix operations](https://learn.arm.com/learning-paths/cross-platform/matrix/4-code-2/)
- [Next Steps](https://learn.arm.com/learning-paths/cross-platform/matrix/_next-steps/)

In the previous section, you created the core, or the boilerplate, needed for `Matrix` objects. They can now be constructed, assigned to/from, and dumped to the screen, but the library is still missing actual matrix operations: add, subtract, and multiply.

In this section, you add the missing operations.

## Operations architecture

When designing matrix processing, it’s desirable to separate two concerns:

- The traversal of the matrices.
- The actual data processing.

This allows you to develop the independent parts and then compose them, allowing for easy testing and extension of the library.

## Unary scalar operations

Unary operations are operations that take a single `Matrix` as argument. The very first unary operations to consider when processing a `Matrix` are:

- `Neg` to negate a value.
- `Abs` to get the absolute value.
- `Sqrt` to get the square root of a value.
- `Log` to get the natural logarithm of a value.

You add these unary scalar operations, in the `MatComp` namespace, but in a separate file.

Add the code below to a new file `include/Matrix/Operators.h`:

```cpp
#pragma once

#include <cmath>
#include <type_traits>

namespace MatComp {

struct unaryOperation { };

/// Negate \p v.
template <typename Ty>
class Neg : public unaryOperation {
  public:
    constexpr Ty operator()(const Ty &v) const { return -v; }
};

/// Get the absolute value of \p v.
template <typename Ty>
class Abs : public unaryOperation {
  public:
    template <typename T = Ty>
    constexpr std::enable_if_t<std::is_unsigned<T>::value, Ty> operator()(
        const Ty &v) const {
        return v;
    }
    template <typename T = Ty>
    constexpr std::enable_if_t<std::is_signed<T>::value, Ty> operator()(
        const Ty &v) const {
        return std::abs(v);
    }
};

/// Get the square root of \p v.
template <typename Ty>
class Sqrt : public unaryOperation {
  public:
    constexpr Ty operator()(const Ty &v) const { return std::sqrt(v); }
};

/// Get the natural logarithm of \p v.
template <typename Ty>
class Log : public unaryOperation {
  public:
    constexpr Ty operator()(const Ty &v) const { return std::log(v); }
};

} // namespace MatComp
```

All these unary operations derive from the empty `unaryOperation` `struct`, which enable code bases to type/group those operations. This essentially *tags* these operations. They are all implemented as templates, in order for them to be type-generic. Although none of these carry a state, they have been implemented as classes with the function operator (`operator()`) defined. This makes them suitable for using in bigger algorithms and is a common pattern used in the C++ standard library.

One point worth mentioning is related to the `Abs` class: depending on the type used at instantiation, the compiler selects an optimized implementation for unsigned types, as there is no need to compute the absolute value of an always positive value. This optimization is transparent to users.

Those operators are marked as `constexpr` so that the compiler can optimize the computation by performing them at compile time rather than runtime if the actual values are known at compile time. They are marked `const`, because the operator does not modify the object state - as these objects don’t have a state.

The new operations must have tests, so add those into a separate `tests/Operators.cpp` file:

```cpp
#include "Matrix/Operators.h"

#include <limits>

#include "gtest/gtest.h"

using std::numeric_limits;

TEST(unaryOperator, Neg) {
    EXPECT_EQ(MatComp::Neg<uint8_t>()(1), uint8_t(-1));
    EXPECT_EQ(MatComp::Neg<uint8_t>()(-2), uint8_t(2));
    EXPECT_EQ(MatComp::Neg<uint16_t>()(1), uint16_t(-1));
    EXPECT_EQ(MatComp::Neg<uint16_t>()(-2), uint16_t(2));
    EXPECT_EQ(MatComp::Neg<uint32_t>()(1), uint32_t(-1));
    EXPECT_EQ(MatComp::Neg<uint32_t>()(-2), uint32_t(2));
    EXPECT_EQ(MatComp::Neg<uint64_t>()(1), uint64_t(-1));
    EXPECT_EQ(MatComp::Neg<uint64_t>()(-2), uint64_t(2));
    
    EXPECT_EQ(MatComp::Neg<int8_t>()(1), int8_t(-1));
    EXPECT_EQ(MatComp::Neg<int8_t>()(-2), int8_t(2));
    EXPECT_EQ(MatComp::Neg<int16_t>()(1), int16_t(-1));
    EXPECT_EQ(MatComp::Neg<int16_t>()(-2), int16_t(2));
    EXPECT_EQ(MatComp::Neg<int32_t>()(1), int32_t(-1));
    EXPECT_EQ(MatComp::Neg<int32_t>()(-2), int32_t(2));
    EXPECT_EQ(MatComp::Neg<int64_t>()(1), int64_t(-1));
    EXPECT_EQ(MatComp::Neg<int64_t>()(-2), int64_t(2));

    EXPECT_FLOAT_EQ(MatComp::Neg<float>()(1.0), float(-1.0));
    EXPECT_FLOAT_EQ(MatComp::Neg<float>()(-2.0), float(2.0));
    EXPECT_DOUBLE_EQ(MatComp::Neg<double>()(1.0), double(-1.0));
    EXPECT_DOUBLE_EQ(MatComp::Neg<double>()(-2.0), double(2.0));
}

TEST(unaryOperator, Abs) {
    EXPECT_EQ(MatComp::Abs<uint8_t>()(1), uint8_t(1));
    EXPECT_EQ(MatComp::Abs<uint8_t>()(-2), numeric_limits<uint8_t>::max() - 2 + 1);
    EXPECT_EQ(MatComp::Abs<uint16_t>()(1), uint16_t(1));
    EXPECT_EQ(MatComp::Abs<uint16_t>()(-2), numeric_limits<uint16_t>::max() - 2 + 1);
    EXPECT_EQ(MatComp::Abs<uint32_t>()(1), uint32_t(1));
    EXPECT_EQ(MatComp::Abs<uint32_t>()(-2), numeric_limits<uint32_t>::max() - 2 + 1);
    EXPECT_EQ(MatComp::Abs<uint64_t>()(1), uint64_t(1));
    EXPECT_EQ(MatComp::Abs<uint64_t>()(-2), numeric_limits<uint64_t>::max() - 2 + 1);
    
    EXPECT_EQ(MatComp::Abs<int8_t>()(1), int8_t(1));
    EXPECT_EQ(MatComp::Abs<int8_t>()(-2), int8_t(2));
    EXPECT_EQ(MatComp::Abs<int16_t>()(1), int16_t(1));
    EXPECT_EQ(MatComp::Abs<int16_t>()(-2), int16_t(2));
    EXPECT_EQ(MatComp::Abs<int32_t>()(1), int32_t(1));
    EXPECT_EQ(MatComp::Abs<int32_t>()(-2), int32_t(2));
    EXPECT_EQ(MatComp::Abs<int64_t>()(1), int64_t(1));
    EXPECT_EQ(MatComp::Abs<int64_t>()(-2), int64_t(2));

    EXPECT_FLOAT_EQ(MatComp::Abs<float>()(1.0), float(1.0));
    EXPECT_FLOAT_EQ(MatComp::Abs<float>()(-2.0), float(2.0));
    EXPECT_DOUBLE_EQ(MatComp::Abs<double>()(1.0), double(1.0));
    EXPECT_DOUBLE_EQ(MatComp::Abs<double>()(-2.0), double(2.0));
}

TEST(unaryOperator, Sqrt) {
    EXPECT_EQ(MatComp::Sqrt<uint8_t>()(4), uint8_t(2));
    EXPECT_EQ(MatComp::Sqrt<uint16_t>()(16), uint16_t(4));
    EXPECT_EQ(MatComp::Sqrt<uint32_t>()(25), uint32_t(5));
    EXPECT_EQ(MatComp::Sqrt<uint64_t>()(36), uint64_t(6));
    
    EXPECT_FLOAT_EQ(MatComp::Sqrt<float>()(4.0), float(2.0));
    EXPECT_FLOAT_EQ(MatComp::Sqrt<float>()(9.0), float(3.0));
    EXPECT_DOUBLE_EQ(MatComp::Sqrt<double>()(4.0), double(2.0));
    EXPECT_DOUBLE_EQ(MatComp::Sqrt<double>()(9.0), double(3.0));
}

TEST(unaryOperator, Log) {
    EXPECT_EQ(MatComp::Log<uint8_t>()(1), uint8_t(0));
    EXPECT_EQ(MatComp::Log<uint16_t>()(64), uint16_t(6));
    EXPECT_EQ(MatComp::Log<uint32_t>()(100), uint32_t(2));
    EXPECT_EQ(MatComp::Log<uint64_t>()(256), uint64_t(8));

    EXPECT_FLOAT_EQ(MatComp::Log<float>()(1.0), float(0.0));
    EXPECT_FLOAT_EQ(MatComp::Log<float>()(64.0), float(4.158883));
    EXPECT_DOUBLE_EQ(MatComp::Log<double>()(1.0), double(0.0));
    EXPECT_DOUBLE_EQ(MatComp::Log<double>()(64.0), double(4.158883));
}
```

## Whole matrix operations: neg, abs, sqrt, log

With the unary scalar operations in place, you can now add some matrix processing.

All the whole matrix operations have the same pattern of operation: each element in the matrix needs to be transformed with the unary scalar operation, the order in which the elements are modified does not matter, so there is no need to apply the scalar operations taking into account rows and columns. Giving the compiler one simple big loop gives it the most optimization opportunities.

First, create a `applyEltWiseUnaryOp` helper routine in the public section of `Matrix` in `include/Matrix/Matrix.h` which is templated on the desired scalar operation as follows:

```cpp
/// Apply element wise unary scalar operator \p op to each element.
template <template <typename> class uOp>
Matrix &applyEltWiseUnaryOp(const uOp<Ty> &op) {
    static_assert(std::is_base_of<unaryOperation, uOp<Ty>>::value,
                  "op must be a unaryOperation");
    for (size_t i = 0; i < getNumElements(); i++)
        data[i] = op(data[i]);
    return *this;
}
```

Although `applyEltWiseUnaryOp` does not need to be a public member, it allows users to add their own operations should they lack one specific to their usage. In other words, this provides easy extendability for the library. `applyEltWiseUnaryOp` also checks at compile time, with the `static_assert` statement, that the operation that it has to perform is a unary operation.

Now include `Matrix/Operators.h` at the top of `Matrix/Matrix.h` so you can use the `Neg`, `Abs`, `Sqrt` and `Log` with `applyEltWiseUnaryOp`. In the `Matrix` public section in `include/Matrix/Matrix.h`, add the `neg`, `abs`, `sqrt` and `log` member operations:

```cpp
/// Apply negate to each element of this Matrix.
Matrix &neg() { return applyEltWiseUnaryOp(Neg<Ty>()); }

/// Apply absolute value to each element of this Matrix.
Matrix &abs() { return applyEltWiseUnaryOp(Abs<Ty>()); }

/// Apply square root to each element of this Matrix.
Matrix &sqrt() { return applyEltWiseUnaryOp(Sqrt<Ty>()); }

/// Apply natural logarithm to each element of this Matrix.
Matrix &log() { return applyEltWiseUnaryOp(Log<Ty>()); }
```

Those member operations modify the object in place, which might not be desirable, so you also provide a functional version of these operators, which is a version that modifies a copy of the object instead. This is achieved by providing functions outside of the `Matrix` class - but still in the `MatComp` namespace. Add these after the `Matrix` class closing `};` in `include/Matrix/Matrix.h`:

```cpp
/// Functional version of neg.
template <typename Ty>
Matrix<Ty> neg(const Matrix<Ty> &m) {
    return Matrix(m).neg();
}

/// Functional version of abs.
template <typename Ty>
Matrix<Ty> abs(const Matrix<Ty> &m) {
    return Matrix(m).abs();
}

/// Functional version of sqrt.
template <typename Ty>
Matrix<Ty> sqrt(const Matrix<Ty> &m) {
    return Matrix(m).sqrt();
}

/// Functional version of log.
template <typename Ty>
Matrix<Ty> log(const Matrix<Ty> &m) {
    return Matrix(m).log();
}
```

Users can now write either:

- `m.abs()` which modifies matrix `m` in-place so that it contains the absolute value of each of its elements.
- `abs(m)` which returns a copy of `m` with the absolute value of each of `m` elements, leaving `m` untouched.

Of course, each of these functions needs to have tests.

Add the following tests in `tests/Matrix.cpp`:

```cpp
TEST(Matrix, neg) {
    Matrix<int16_t> m0(2, 2, {1, -2, -3, 4});
    Matrix<int16_t> m1 = neg(m0); // Functional version
    EXPECT_EQ(m1.get(0, 0), -1);
    EXPECT_EQ(m1.get(0, 1), 2);
    EXPECT_EQ(m1.get(1, 0), 3);
    EXPECT_EQ(m1.get(1, 1), -4);

    m0.neg(); // In-place version
    EXPECT_EQ(m0.get(0, 0), -1);
    EXPECT_EQ(m0.get(0, 1), 2);
    EXPECT_EQ(m0.get(1, 0), 3);
    EXPECT_EQ(m0.get(1, 1), -4);
}
```

## Binary scalar operations

Users also need to have binary operations; operations that have two inputs: addition, subtraction, and multiplication are the most popular representatives of binary operations. Using the same pattern as for unary operations, add those binary operations to `include/Matrix/Operators.h` in the `MatComp` namespace:

```cpp
struct binaryOperation { };

/// Add \p lhs and \p rhs.
template <typename Ty>
class Add : public binaryOperation {
  public:
    constexpr Ty operator()(const Ty &lhs, const Ty &rhs) const {
        return lhs + rhs;
    }
};

/// Subtract \p rhs from \p lhs.
template <typename Ty>
class Sub : public binaryOperation {
  public:
    constexpr Ty operator()(const Ty &lhs, const Ty &rhs) const {
        return lhs - rhs;
    }
};

/// Multiply \p lhs by \p rhs.
template <typename Ty>
class Mul : public binaryOperation {
  public:
    constexpr Ty operator()(const Ty &lhs, const Ty &rhs) const {
        return lhs * rhs;
    }
};
```

Add the tests to `tests/Operators.cpp`:

```cpp
TEST(binaryOperator, Add) {
    EXPECT_EQ(MatComp::Add<uint8_t>()(1, 2), uint8_t(3));
    EXPECT_EQ(MatComp::Add<uint16_t>()(3, 4), uint16_t(7));
    EXPECT_EQ(MatComp::Add<uint32_t>()(5, 6), uint32_t(11));
    EXPECT_EQ(MatComp::Add<uint64_t>()(7, 8), uint64_t(15));

    EXPECT_EQ(MatComp::Add<int8_t>()(9, 10), int8_t(19));
    EXPECT_EQ(MatComp::Add<int16_t>()(11, 12), int16_t(23));
    EXPECT_EQ(MatComp::Add<int32_t>()(13, 14), int32_t(27));
    EXPECT_EQ(MatComp::Add<int64_t>()(15, 16), int64_t(31));

    EXPECT_FLOAT_EQ(MatComp::Add<float>()(8.0, 12.0), float(20.0));
    EXPECT_DOUBLE_EQ(MatComp::Add<double>()(32.0, 64.0), double(96.0));
}

TEST(binaryOperator, Sub) {
    EXPECT_EQ(MatComp::Sub<uint8_t>()(2, 1), uint8_t(1));
    EXPECT_EQ(MatComp::Sub<uint16_t>()(3, 1), uint16_t(2));
    EXPECT_EQ(MatComp::Sub<uint32_t>()(5, 2), uint32_t(3));
    EXPECT_EQ(MatComp::Sub<uint64_t>()(8, 3), uint64_t(5));

    EXPECT_EQ(MatComp::Sub<int8_t>()(9, 10), int8_t(-1));
    EXPECT_EQ(MatComp::Sub<int16_t>()(11, 11), int16_t(0));
    EXPECT_EQ(MatComp::Sub<int32_t>()(14, 14), int32_t(0));
    EXPECT_EQ(MatComp::Sub<int64_t>()(15, 16), int64_t(-1));

    EXPECT_FLOAT_EQ(MatComp::Sub<float>()(8.0, 12.0), float(-4.0));
    EXPECT_DOUBLE_EQ(MatComp::Sub<double>()(32.0, 64.0), double(-32.0));
}

TEST(binaryOperator, Mul) {
    EXPECT_EQ(MatComp::Mul<uint8_t>()(1, 2), uint8_t(2));
    EXPECT_EQ(MatComp::Mul<uint16_t>()(3, 4), uint16_t(12));
    EXPECT_EQ(MatComp::Mul<uint32_t>()(5, 6), uint32_t(30));
    EXPECT_EQ(MatComp::Mul<uint64_t>()(7, 8), uint64_t(56));

    EXPECT_EQ(MatComp::Mul<int8_t>()(9, 10), int8_t(90));
    EXPECT_EQ(MatComp::Mul<int16_t>()(11, 12), int16_t(132));
    EXPECT_EQ(MatComp::Mul<int32_t>()(13, 14), int32_t(182));
    EXPECT_EQ(MatComp::Mul<int64_t>()(15, 16), int64_t(240));

    EXPECT_FLOAT_EQ(MatComp::Mul<float>()(8.0, 12.0), float(96.0));
    EXPECT_DOUBLE_EQ(MatComp::Mul<double>()(32.0, 64.0), double(2048.0));
}
```

Check the code compiles and passes the checks:

```bash
cd build
ninja
ninja check
```

## Element wise Matrix operations: add, subtract, multiply

Now that scalar binary operators are implemented, you can add `Matrix` level binary operators. The pattern is very similar to the one used for unary operators: you add a helper routine that traverses the matrices and applies the scalar operator, but with an extra step. Although the element-wise processing of matrices does not require a specific order for the rows or columns, the matrices must have the same dimensions; if not, this is an error and the program must be terminated.

Add this `applyEltWiseBinaryOp` helper routine to the public section of `Matrix` in `include/Matrix/Matrix.h`:

```cpp
/// Apply element wise binary scalar operator \p op to each element.
template <template <typename> class bOp>
Matrix &applyEltWiseBinaryOp(const bOp<Ty> &op, const Matrix &rhs) {
    static_assert(std::is_base_of<binaryOperation, bOp<Ty>>::value,
                  "op must be a binaryOperation");
    if (getNumRows() != rhs.getNumRows() || getNumColumns() != rhs.getNumColumns())
        die(__FILE__, __LINE__, "Inconsistent Matrix dimensions for elementwise operation");

    for (size_t i = 0; i < getNumElements(); i++)
        data[i] = op(data[i], rhs.data[i]);
    return *this;
}
```

Add to the public section of `Matrix` in `include/Matrix/Matrix.h` the `add`, `sub` and `mul` operators, composing the matrices traversal `applyEltWiseBinaryOp` with the scalar operators you’ve just added:

```cpp
// Add Matrix \p rhs to this Matrix.
Matrix &operator+=(const Matrix &rhs) {
    return applyEltWiseBinaryOp(Add<Ty>(), rhs);
}
// Subtract Matrix \p rhs from this Matrix.
Matrix &operator-=(const Matrix &rhs) {
    return applyEltWiseBinaryOp(Sub<Ty>(), rhs);
}
// Element-wise multiplication of this Matrix by \p rhs (a.k.a dot product).
Matrix &operator*=(const Matrix &rhs) {
    return applyEltWiseBinaryOp(Mul<Ty>(), rhs);
}
```

These operators do modify the objects in place, which explains the `+=`, `-=`, and `*=`. You now add functional versions of these operations. Outside of the `Matrix` class in `include/Matrix/Operators.h`, but in the `MatComp` namespace, add:

```cpp
/// Add Matrix \p lhs and \p rhs.
template <typename Ty>
Matrix<Ty> operator+(const Matrix<Ty> &lhs, const Matrix<Ty> &rhs) {
    return Matrix(lhs) += rhs;
}

/// Subtract Matrix \p rhs from Matrix \p lhs.
template <typename Ty>
Matrix<Ty> operator-(const Matrix<Ty> &lhs, const Matrix<Ty> &rhs) {
    return Matrix(lhs) -= rhs;
}

/// Elementwise multiplication (a.k.a dot product) of Matrix \p lhs by \p rhs.
template <typename Ty>
Matrix<Ty> operator*(const Matrix<Ty> &lhs, const Matrix<Ty> &rhs) {
    return Matrix(lhs) *= rhs;
}
```

You can see that those operations first make a copy of the `lhs` object with `Matrix(lhs)`, then modify the copy in place with `rhs` and then return the result without modifying `lhs` or `rhs`.

Add tests for all these in-place and functional operators in `tests/Matrix.cpp`:

```cpp
TEST(Matrix, inPlaceAdd) {
    Matrix<int16_t> m0 = Matrix<int16_t>::ones(3, 3);
    const Matrix<int16_t> m1 = Matrix<int16_t>::identity(3);
    m0 += m1;
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 3);
    EXPECT_EQ(m0.getNumColumns(), 3);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            if (row == col)
                EXPECT_EQ(m0.get(row, col), 2);
            else
                EXPECT_EQ(m0.get(row, col), 1);
        }
    }
}

TEST(Matrix, inPlaceSub) {
    Matrix<int32_t> m0 = Matrix<int32_t>::ones(4, 4);
    const Matrix<int32_t> m1 = Matrix<int32_t>::identity(4);
    m0 -= m1;
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 4);
    EXPECT_EQ(m0.getNumColumns(), 4);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            if (row == col)
                EXPECT_EQ(m0.get(row, col), 0);
            else
                EXPECT_EQ(m0.get(row, col), 1);
        }
    }
}

TEST(Matrix, inPlaceDotProduct) {
    Matrix<int64_t> m0 = Matrix<int64_t>(2, 2, {1, 2, 3, 4});
    const Matrix<int64_t> m1 = Matrix<int64_t>::identity(2);
    m0 *= m1;
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 2);
    EXPECT_EQ(m0.getNumColumns(), 2);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            if (row == col)
                EXPECT_EQ(m0.get(row, col), row * m0.getNumColumns() + col + 1);
            else
                EXPECT_EQ(m0.get(row, col), 0);
        }
    }
}

TEST(Matrix, functionalAdd) {
    Matrix<int16_t> m0 = Matrix<int16_t>::ones(3, 3) + Matrix<int16_t>::identity(3);
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 3);
    EXPECT_EQ(m0.getNumColumns(), 3);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            if (row == col)
                EXPECT_EQ(m0.get(row, col), 2);
            else
                EXPECT_EQ(m0.get(row, col), 1);
        }
    }
}

TEST(Matrix, functionalSub) {
    Matrix<int32_t> m0 = Matrix<int32_t>::ones(4, 4) - Matrix<int32_t>::identity(4);
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 4);
    EXPECT_EQ(m0.getNumColumns(), 4);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            if (row == col)
                EXPECT_EQ(m0.get(row, col), 0);
            else
                EXPECT_EQ(m0.get(row, col), 1);
        }
    }
}

TEST(Matrix, functionalDotProduct) {
    Matrix<int64_t> m0 = Matrix<int64_t>(2, 2, {1, 2, 3, 4}) * Matrix<int64_t>::identity(2);
    EXPECT_TRUE(m0);
    EXPECT_EQ(m0.getNumRows(), 2);
    EXPECT_EQ(m0.getNumColumns(), 2);
    for (size_t row = 0; row < m0.getNumRows(); row++) {
        for (size_t col = 0; col < m0.getNumColumns(); col++) {
            EXPECT_EQ(m0.get(row, col), row * m0.getNumColumns() + col + 1);
        }
    }
}
```

## Matrix x Matrix multiplication

Last, but not least, the `Matrix` class still requires the most important operation: matrix multiplication.

Add the `multiply` function, outside of the `Matrix` class, but in the `MatComp` namespace in `include/Matrix/Matrix.h`. `multiply` checks that the arguments have compatible dimensions, then creates an uninitialized matrix that holds the multiplication result and eventually performs the multiplication to fill it.

```cpp
/// Matrix multiplication.
template <typename Ty>
Matrix<Ty> multiply(const Matrix<Ty> &lhs, const Matrix<Ty> &rhs) {
    if (lhs.getNumColumns() != rhs.getNumRows())
        die(__FILE__, __LINE__, "Incompatible dimensions in Matrix multiply");
    Matrix<Ty> result(lhs.getNumRows(), rhs.getNumColumns());
    for (size_t row = 0; row < result.getNumRows(); row++) {
        for (size_t col = 0; col < result.getNumColumns(); col++) {
            Ty acc = Ty(0);
            for (size_t t = 0; t < lhs.getNumColumns(); t++)
                acc += lhs.get(row, t) * rhs.get(t, col);
            result.get(row, col) = acc;
        }
    }
    return result;
}
```

Add some tests in `tests/Matrix.cpp`:

```cpp
TEST(Matrix, multiplication) {
    // Test matrix x matrix.
    Matrix<double> m0 = multiply(Matrix<double>(2, 2, {1., 2., 3., 4.}),
                                  Matrix<double>(2, 2, {0., 1., 2., 3.}));
    EXPECT_TRUE(m0);
    const Matrix<double> m0_exp(2, 2, {4., 7., 8., 15.});
    EXPECT_EQ(m0.getNumColumns(), m0_exp.getNumColumns());
    EXPECT_EQ(m0.getNumRows(), m0_exp.getNumRows());
    EXPECT_EQ(m0.getNumElements(), m0_exp.getNumElements());
    for (size_t row = 0; row < m0_exp.getNumRows(); row++) {
        for (size_t col = 0; col < m0_exp.getNumColumns(); col++) {
            EXPECT_DOUBLE_EQ(m0.get(row, col), m0_exp.get(row, col));
        }
    }

    // Test matrix x vector
    Matrix<float> m1 = multiply(Matrix<float>(3, 3, {1., 2., 3., 4., 5., 6., 7., 8., 9.}),
                                 Matrix<float>(3, 1, {1., 2., 3.}));
    EXPECT_TRUE(m1);
    const Matrix<float> m1_exp(3, 1, {14., 32., 50.});
    EXPECT_EQ(m1.getNumColumns(), m1_exp.getNumColumns());
    EXPECT_EQ(m1.getNumRows(), m1_exp.getNumRows());
    EXPECT_EQ(m1.getNumElements(), m1_exp.getNumElements());
    for (size_t row = 0; row < m1_exp.getNumRows(); row++) {
        for (size_t col = 0; col < m1_exp.getNumColumns(); col++) {
            EXPECT_FLOAT_EQ(m1.get(row, col), m1_exp.get(row, col));
        }
    }
}
```

Check the code compiles and passes the checks:

```bash
cd build
ninja
ninja check
```

## Other useful operations

A complete matrix processing library requires some other important operations:

- Type conversion: the library only supports operations with the same type for now, and this is enforced by the compiler. Avoiding silent conversions is desirable, but there are cases where data have different types and *explicit* conversion is necessary.
- Transposition, that is swapping the rows and columns, is another important operation, with some room for optimization in the implementation. Can you think of different ways to implement the `transpose` operator for the `Matrix` class?
- Broadcasting: these are useful shortcuts for the `Matrix` users, as they allow to transform a scalar or a row-vector/column-vector by replicating it to make it look like a 2-D matrix. Again, on top of syntactic sugar for the users, this provides a large performance improvement, especially in memory usage, as one can play with how the elements are accessed rather than duplicating the matrix content.
- Resize: to be able to dynamically change a matrix dimensions.
- Extract: to be able to extract part of a matrix.

## Optimization

The code written so far is relatively high level and allows the compiler to perform a large number of optimizations, from propagating constants to unrolling loops to name but a few most basic ones.

The `applyEltWiseUnaryOp` and `applyEltWiseBinaryOp` helper routines from the Matrix library process one element at a time. The compiler may make use of Arm specific SIMD (Single Instruction, Multiple Data) instructions to process several elements at a time with one instruction. This is an optimization named vectorization, that can either be done automatically by the compiler (this is named *autovectorization*) or it can be done manually by the developer with the use of *intrinsics functions*. You can learn more about the compiler’s autovectorization capabilities with the [Learn about Autovectorization](https://learn.arm.com/learning-paths/cross-platform/loop-reflowing/) learning path and about other vectorization tricks with the [Optimize SIMD code with vectorization-friendly data layout](https://learn.arm.com/learning-paths/cross-platform/vectorization-friendly-data-layout/) learning path.

You can also learn how to [Accelerate Matrix Multiplication Performance with SME2](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/).

## What have you achieved so far?

At this stage, the code structure looks like:

```
Matrix/
├── CMakeLists.txt
├── build/
...
├── external/
│   └── CMakeLists.txt
├── include/
│   └── Matrix/
│       ├── Matrix.h
│       └── Operators.h
├── lib/
│   └── Matrix/
│       └── Matrix.cpp
├── src/
│   ├── getVersion.cpp
│   └── howdy.cpp
└── tests/
    ├── Matrix.cpp
    ├── Operators.cpp
    ├── Version.cpp
    └── main.cpp
```

Congratulations, you now have a minimalistic, yet fully-functional, matrix processing library, with some level of regression testing, that can be easily built and used.

The testing could - and *should* - go much deeper, as a number of corner cases have not been covered.

You can continue to add more functions, and more tests.

You can refer to this chapter source code in `code-examples/learning-paths/cross-platform/matrix/chapter-4` in the archive that you have downloaded earlier.