From b7c3ba83fb1167d6103b85a27be9003c34f63ef5 Mon Sep 17 00:00:00 2001 From: Biljana Bojovic Date: Thu, 15 Dec 2022 15:32:24 +0100 Subject: [PATCH] core: Add MatrixArray and MatrixArrayTestSuite MatrixArray class is introduced for efficient operations on arrays of matrices core: Allow using GetPagePtr as public function of MatrixArray core: Remove Eigen include from matrix-array.h core: Fix MatrixArray Doxygen core: Fix MatrixArrayTestSuite Doxygen core: Make HermitianTranspose member of MatrixArray fix matrix-array spell check core: Improve constructor Doxygen description core: Move out MatrixArray inline function implementation from the class declaration (Thanks to Peter Barnes) Use auto in MatrixArrayTestCase to be consistent --- src/core/CMakeLists.txt | 3 + src/core/model/matrix-array.cc | 227 ++++++++++++++ src/core/model/matrix-array.h | 308 +++++++++++++++++++ src/core/test/matrix-array-test-suite.cc | 369 +++++++++++++++++++++++ 4 files changed, 907 insertions(+) create mode 100644 src/core/model/matrix-array.cc create mode 100644 src/core/model/matrix-array.h create mode 100644 src/core/test/matrix-array-test-suite.cc diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index e2f731cb9..ad693e31b 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -209,6 +209,7 @@ set(source_files model/trickle-timer.cc model/realtime-simulator-impl.cc model/wall-clock-synchronizer.cc + model/matrix-array.cc ) # Define core lib headers @@ -314,6 +315,7 @@ set(header_files model/realtime-simulator-impl.h model/wall-clock-synchronizer.h model/val-array.h + model/matrix-array.h ) set(test_sources @@ -350,6 +352,7 @@ set(test_sources test/type-traits-test-suite.cc test/watchdog-test-suite.cc test/val-array-test-suite.cc + test/matrix-array-test-suite.cc ) # Build core lib diff --git a/src/core/model/matrix-array.cc b/src/core/model/matrix-array.cc new file mode 100644 index 000000000..6d63591cf --- /dev/null +++ b/src/core/model/matrix-array.cc @@ -0,0 +1,227 @@ +/* + * Copyright (c) 2022 Centre Tecnologic de Telecomunicacions de Catalunya (CTTC) + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License version 2 as + * published by the Free Software Foundation; + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * Author: Biljana Bojovic + */ + +#include "matrix-array.h" + +#ifdef HAVE_EIGEN3 +#include +#endif + +namespace ns3 +{ + +#ifdef HAVE_EIGEN3 +template +using EigenMatrix = Eigen::Map>; +template +using ConstEigenMatrix = Eigen::Map>; +#endif + +template +MatrixArray::MatrixArray(uint16_t numRows, uint16_t numCols, uint16_t numPages) + : ValArray(numRows, numCols, numPages){}; + +template +MatrixArray::MatrixArray(const std::valarray& values) + : ValArray(values) +{ +} + +template +MatrixArray::MatrixArray(std::valarray&& values) + : ValArray(std::move(values)) +{ +} + +template +MatrixArray::MatrixArray(const std::vector& values) + : ValArray(values) +{ +} + +template +MatrixArray::MatrixArray(uint16_t numRows, uint16_t numCols, const std::valarray& values) + : ValArray(numRows, numCols, values){}; + +template +MatrixArray::MatrixArray(uint16_t numRows, uint16_t numCols, std::valarray&& values) + : ValArray(numRows, numCols, std::move(values)){}; + +template +MatrixArray::MatrixArray(uint16_t numRows, + uint16_t numCols, + uint16_t numPages, + const std::valarray& values) + : ValArray(numRows, numCols, numPages, values){}; + +template +MatrixArray::MatrixArray(uint16_t numRows, + uint16_t numCols, + uint16_t numPages, + std::valarray&& values) + : ValArray(numRows, numCols, numPages, std::move(values)){}; + +template +MatrixArray +MatrixArray::operator*(const MatrixArray& rhs) const +{ + NS_ASSERT_MSG(m_numPages == rhs.m_numPages, "MatrixArrays have different numbers of matrices."); + NS_ASSERT_MSG(m_numCols == rhs.m_numRows, "Inner dimensions of matrices mismatch."); + + MatrixArray res{m_numRows, rhs.m_numCols, m_numPages}; + + for (auto page = 0; page < res.m_numPages; ++page) + { +#ifdef HAVE_EIGEN3 // Eigen found and enabled Eigen optimizations + + ConstEigenMatrix lhsEigenMatrix(GetPagePtr(page), m_numRows, m_numCols); + ConstEigenMatrix rhsEigenMatrix(rhs.GetPagePtr(page), rhs.m_numRows, rhs.m_numCols); + EigenMatrix resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols); + resEigenMatrix = lhsEigenMatrix * rhsEigenMatrix; + +#else // Eigen not found or Eigen optimizations not enabled + + size_t matrixOffset = page * m_numRows * m_numCols; + for (auto i = 0; i < res.m_numRows; ++i) + { + for (auto j = 0; j < res.m_numCols; ++j) + { + res(i, j, page) = + (m_values[std::slice(matrixOffset + i, m_numCols, m_numRows)] * + rhs.m_values[std::slice(matrixOffset + j * rhs.m_numRows, rhs.m_numRows, 1)]) + .sum(); + } + } + +#endif + } + return res; +} + +template +MatrixArray +MatrixArray::Transpose() const +{ + MatrixArray res{m_numCols, + m_numRows, + m_numPages}; // create the matrix where m_numRows = this.m_numCols, m_numCols + // = this.m_numRows, m_numPages = this.m_numPages + + for (auto page = 0; page < m_numPages; ++page) + { +#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled + + ConstEigenMatrix thisMatrix(GetPagePtr(page), m_numRows, m_numCols); + EigenMatrix resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols); + resEigenMatrix = thisMatrix.transpose(); + +#else // Eigen not found or Eigen optimizations not enabled + + size_t matrixIndex = page * m_numRows * m_numCols; + for (auto i = 0; i < m_numRows; ++i) + { + res.m_values[std::slice(matrixIndex + i * res.m_numRows, res.m_numRows, 1)] = + m_values[std::slice(matrixIndex + i, m_numCols, m_numRows)]; + } + +#endif + } + return res; +} + +template +MatrixArray +MatrixArray::MultiplyByLeftAndRightMatrix(const MatrixArray& lMatrix, + const MatrixArray& rMatrix) const +{ + NS_ASSERT_MSG(lMatrix.m_numPages == 1 && rMatrix.m_numPages == 1, + "The left and right MatrixArray should have only one page."); + NS_ASSERT_MSG(lMatrix.m_numCols == m_numRows, + "Left vector numCols and this MatrixArray numRows mismatch."); + NS_ASSERT_MSG(m_numCols == rMatrix.m_numRows, + "Right vector numRows and this MatrixArray numCols mismatch."); + + MatrixArray res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages}; + +#ifdef HAVE_EIGEN3 + + ConstEigenMatrix lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols); + ConstEigenMatrix rMatrixEigen(rMatrix.GetPagePtr(0), rMatrix.m_numRows, rMatrix.m_numCols); +#endif + + for (auto page = 0; page < m_numPages; ++page) + { +#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled + + ConstEigenMatrix matrixEigen(GetPagePtr(page), m_numRows, m_numCols); + EigenMatrix resEigenMap(res.GetPagePtr(page), res.m_numRows, res.m_numCols); + + resEigenMap = lMatrixEigen * matrixEigen * rMatrixEigen; + +#else // Eigen not found or Eigen optimizations not enabled + + size_t matrixOffset = page * m_numRows * m_numCols; + for (auto resRow = 0; resRow < res.m_numRows; ++resRow) + { + for (auto resCol = 0; resCol < res.m_numCols; ++resCol) + { + // create intermediate row result, a multiply of resRow row of lMatrix and each + // column of this matrix + std::valarray interRes(m_numCols); + for (auto thisCol = 0; thisCol < m_numCols; ++thisCol) + { + interRes[thisCol] = + (lMatrix + .m_values[std::slice(resRow, lMatrix.m_numCols, lMatrix.m_numRows)] * + m_values[std::slice(matrixOffset + thisCol * m_numRows, m_numRows, 1)]) + .sum(); + } + // multiply intermediate results and resCol column of the rMatrix + res(resRow, resCol, page) = + (interRes * + rMatrix.m_values[std::slice(resCol * rMatrix.m_numRows, rMatrix.m_numRows, 1)]) + .sum(); + } + } +#endif + } + return res; +} + +template +template +MatrixArray +MatrixArray::HermitianTranspose() const +{ + MatrixArray> retMatrix = this->Transpose(); + + for (size_t index = 0; index < this->GetSize(); ++index) + { + retMatrix[index] = std::conj(m_values[index]); + } + return retMatrix; +} + +template MatrixArray> MatrixArray>::HermitianTranspose() + const; +template class MatrixArray>; +template class MatrixArray; +template class MatrixArray; + +} // namespace ns3 diff --git a/src/core/model/matrix-array.h b/src/core/model/matrix-array.h new file mode 100644 index 000000000..b5628d5ae --- /dev/null +++ b/src/core/model/matrix-array.h @@ -0,0 +1,308 @@ +/* + * Copyright (c) 2022 Centre Tecnologic de Telecomunicacions de Catalunya (CTTC) + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License version 2 as + * published by the Free Software Foundation; + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * Author: Biljana Bojovic + */ + +#ifndef MATRIX_ARRAY_H +#define MATRIX_ARRAY_H + +#include "val-array.h" + +#include + +namespace ns3 +{ + +/** + * \ingroup Matrices + * + * \brief MatrixArray class inherits ValArray class and provides additional + * interfaces to ValArray which enable page-wise linear algebra operations for + * arrays of matrices. While ValArray class is focused on efficient storage, + * access and basic algebra element-wise operations on 3D numerical arrays, + * MatrixArray is focused on the efficient algebra operations on arrays of + * matrices. + * + * Each page (2-dimensional array) within a 3D array is considered to + * be a matrix, and operations are performed separately for each page/matrix. + * For example, a multiplication of two MatrixArrays means that a matrix + * multiplication is performed between each page in the first MatrixArray with + * the corresponding page in the second MatrixArray. Some of the operations can + * use the Eigen3 matrix library for speedup. In addition, this class can be + * extended to provide interfaces to the Eigen3 matrix library and enable + * additional linear algebra operations. The objective of this class is to two-fold: + * + * - First, it provides simple interfaces to perform matrix operations in parallel + * on all pages without the need to write an outer loop over all the pages. + * + * - Second, it can improve the performance of such parallel matrix operations, + * as this class uses a single memory allocation and one common set of size + * variables (rows, cols) for all pages, which can be an advantage compared to + * having separate memory allocations and size variables for each page/matrix, + * especially when those matrices are very small. + * + * Important characteristics: + * + * - Matrices are in column-major order. + * + * - Matrices can be 1-dimensional, in the sense that either the number of rows + * or columns is equal to 1. + * + * - The array of matrices can contain just one matrix or more matrices. + * + * - For binary operators, both MatrixArrays must have an equal number of pages, + * as well compatible row and column dimensions between the matrices. + * + * - Binary operators are performed among matrices with the same index from lhs + * and rhs matrix. + * + * - MatrixArray improves computationally complex matrix operations by using + * Eigen library when it is available at compile time. MatrixArray class makes + * use of Eigen Map class, which allows re-using already existing portion of + * memory as an Eigen type, e.g. ValArray can return a pointer to a memory + * where are placed elements of a specific matrix, and this can be passed to + * Eigen Map class which does not do a copy of the elements, but is using these + * elements directly in linear algebra operations. + */ +template +class MatrixArray : public ValArray +{ + public: + // instruct the compiler to generate the default constructor + MatrixArray() = default; + /** + * \brief Constructor "pages" number of matrices that are of size "numRows"x"numCols", + * and are initialized with all-zero elements. + * If only 1 parameter, numRows, is provided then a single 1D array is being created. + * \param numRows the number of rows + * \param numCols the number of columns + * \param numPages the number of pages + */ + MatrixArray(uint16_t numRows, uint16_t numCols = 1, uint16_t numPages = 1); + /** + * \brief Constructor creates a single array of values.size () elements and 1 column, + * and uses std::valarray values to initialize the elements. + * \param values std::valarray that will be used to initialize elements of 1D array + */ + explicit MatrixArray(const std::valarray& values); + /** + * \brief Constructor creates a single array of values.size () elements and 1 column, + * and moves std::valarray values to initialize the elements. + * \param values std::valarray that will be moved to initialize elements of 1D array + */ + MatrixArray(std::valarray&& values); + /** + * \brief Constructor creates a single array of values.size () elements and 1 column, + * and uses std::vector values to initialize the elements. + * \param values std::vector that will be used to initialize elements of 1D array + */ + explicit MatrixArray(const std::vector& values); + /** + * \brief Constructor creates a single matrix of numRows and numCols, and uses + * std::valarray values to initialize the elements. + * \param numRows the number of rows + * \param numCols the number of columns + * \param values std::valarray that will be used to initialize elements of 3D array + */ + MatrixArray(uint16_t numRows, uint16_t numCols, const std::valarray& values); + /** + * \brief Constructor creates a single matrix of numRows and numCols, and moves + * std::valarray values to initialize the elements. + * \param numRows the number of rows + * \param numCols the number of columns + * \param values std::valarray that will be moved to initialize elements of 3D array + */ + MatrixArray(uint16_t numRows, uint16_t numCols, std::valarray&& values); + /** + * \brief Constructor creates the array of numPages matrices of numRows x numCols dimensions, + * and uses std::valarray values to initialize all the elements. + * \param numRows the number of rows + * \param numCols the number of columns + * \param numPages the number of pages + * \param values std::valarray that will be used to initialize elements of 3D array + */ + MatrixArray(uint16_t numRows, + uint16_t numCols, + uint16_t numPages, + const std::valarray& values); + /** + * \brief Constructor creates the array of numPages matrices of numRows x numCols dimensions, + * and moves std::valarray values to initialize all the elements. + * \param numRows the number of rows + * \param numCols the number of columns + * \param numPages the number of pages + * \param values std::valarray that will be used to initialize elements of 3D array + */ + MatrixArray(uint16_t numRows, + uint16_t numCols, + uint16_t numPages, + std::valarray&& values); + /** instruct the compiler to generate the implicitly declared destructor*/ + ~MatrixArray() override = default; + /** instruct the compiler to generate the implicitly declared copy constructor*/ + MatrixArray(const MatrixArray&) = default; + /** + * \brief Copy assignment operator. + * Instruct the compiler to generate the implicitly declared copy assignment operator. + * \return a reference to the result of the assignment + */ + MatrixArray& operator=(const MatrixArray&) = default; + /** instruct the compiler to generate the implicitly declared move constructor*/ + MatrixArray(MatrixArray&&) noexcept = default; + /** + * \brief Move assignment operator. + * Instruct the compiler to generate the implicitly declared move assignment operator. + * \return a reference to the result of the assignment + */ + MatrixArray& operator=(MatrixArray&&) noexcept = default; + /** + * \brief Element-wise multiplication with a scalar value. + * \param rhs is a scalar value of type T + * \returns The matrix in which each element has been multiplied by the given + * scalar value. + */ + MatrixArray operator*(const T& rhs) const; + /** + * \brief operator+ definition for MatrixArray. + * \param rhs The rhs MatrixArray object + * \returns The result of operator+ of this MatrixArray and rhs MatrixAray + */ + MatrixArray operator+(const MatrixArray& rhs) const; + /** + * \brief binary operator- definition for MatrixArray. + * \param rhs The rhs MatrixArray object + * \returns The result of operator- of this MatrixArray and rhs MatrixAray + */ + MatrixArray operator-(const MatrixArray& rhs) const; + /** + * \brief unary operator- definition for MatrixArray. + * \return the result of the operator- + */ + MatrixArray operator-() const; + /** + * \brief Page-wise matrix multiplication. + * This operator interprets the 3D array as an array of matrices, + * and performs a linear algebra operation on each of the matrices (pages), + * i.e., matrices from this MatrixArray and rhs with the same page indices are + * multiplied. + * The number of columns of this MatrixArray must be equal to the number of + * rows in rhs, and rhs must have the same number of pages as this MatrixArray. + * \param rhs is another MatrixArray instance + * \returns The array of results of matrix multiplications. + */ + MatrixArray operator*(const MatrixArray& rhs) const; + /** + * \brief This operator interprets the 3D array as an array of matrices, + * and performs a linear algebra operation on each of the matrices (pages), + * i.e., transposes the matrix or array of matrices definition for MatrixArray. + * \return The resulting MatrixArray composed of the array of transposed matrices. + */ + MatrixArray Transpose() const; + /** + *\brief Multiply each matrix in the array by the left and the right matrix. + * For each page of this MatrixArray the operation performed is + * lMatrix * matrix(pageIndex) * rMatrix, and the resulting MatrixArray + * contains the array of the results per page. If "this" has dimensions + * M x N x P, lMatrix must have dimensions M number of columns, and a single + * page, and rMatrix must have N number of rows, and also a single page. + * + * The result of this function will have the dimensions J x K x P. Where J is + * the number of rows of the lMatrix, and K is the number of + * columns of rMatrix. Dimensions are: + * + * lMatrix(JxMx1) * this (MxNxP) * rMatrix(NxKx1) -> result(JxKxP) + * + * This operation is not possible when using the multiplication operator because + * the number of pages does not match. + * + * \param lMatrix the left matrix in the multiplication + * \param rMatrix the right matrix in the multiplication + * \returns Returns the result of the multiplication which is a 3D MatrixArray + * with dimensions J x K x P as explained previously. + */ + MatrixArray MultiplyByLeftAndRightMatrix(const MatrixArray& lMatrix, + const MatrixArray& rMatrix) const; + + using ValArray::GetPagePtr; + using ValArray::EqualDims; + using ValArray::AssertEqualDims; + + /** + * \brief Function that performs the Hermitian transpose of this MatrixArray + * and returns a new matrix that is the result of the operation. + * This function assumes that the matrix contains complex numbers, because of that + * it is only available for the > specialization of MatrixArray. + * \return Returns a new matrix that is the result of the Hermitian transpose + */ + template >::value && + EnableBool)>::type> + MatrixArray HermitianTranspose() const; + + private: + // To simplify functions in MatrixArray that are using members from the template base class + using ValArray::m_numRows; + using ValArray::m_numCols; + using ValArray::m_numPages; + using ValArray::m_values; +}; + +using IntMatrixArray = MatrixArray; //!< Create an alias for MatrixArray using int type +using DoubleMatrixArray = + MatrixArray; //!< Create an alias for MatrixArray using double type +using ComplexMatrixArray = + MatrixArray>; //!< Create an alias for MatrixArray using complex type + +/************************************************* + ** Class MatrixArray inline implementations + ************************************************/ +template +inline MatrixArray +MatrixArray::operator*(const T& rhs) const +{ + return MatrixArray(m_numRows, + m_numCols, + m_numPages, + m_values * std::valarray(rhs, m_numRows * m_numCols * m_numPages)); +} + +template +inline MatrixArray +MatrixArray::operator+(const MatrixArray& rhs) const +{ + AssertEqualDims(rhs); + return MatrixArray(m_numRows, m_numCols, m_numPages, m_values + rhs.m_values); +} + +template +inline MatrixArray +MatrixArray::operator-(const MatrixArray& rhs) const +{ + AssertEqualDims(rhs); + return MatrixArray(m_numRows, m_numCols, m_numPages, m_values - rhs.m_values); +} + +template +inline MatrixArray +MatrixArray::operator-() const +{ + return MatrixArray(m_numRows, m_numCols, m_numPages, -m_values); +} + +} // namespace ns3 + +#endif // MATRIX_ARRAY_H diff --git a/src/core/test/matrix-array-test-suite.cc b/src/core/test/matrix-array-test-suite.cc new file mode 100644 index 000000000..c47cbe9dd --- /dev/null +++ b/src/core/test/matrix-array-test-suite.cc @@ -0,0 +1,369 @@ +/* + * Copyright (c) 2022 Centre Tecnologic de Telecomunicacions de Catalunya (CTTC) + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License version 2 as + * published by the Free Software Foundation; + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * + * Author: Biljana Bojovic + */ + +#include "ns3/log.h" +#include "ns3/matrix-array.h" +#include "ns3/test.h" + +/** + * \file + * \ingroup core-tests + * \ingroup matrixArray + * \ingroup matrixArray-tests + * MatrixArray test suite + */ +namespace ns3 +{ + +namespace tests +{ + +NS_LOG_COMPONENT_DEFINE("MatrixArrayTest"); + +/** + * \ingroup matrixArray-tests + * MatrixArray test for testing constructors + */ +template +class MatrixArrayTestCase : public TestCase +{ + public: + MatrixArrayTestCase() = default; + /** + * Constructor + * + * \param [in] name reference name + */ + MatrixArrayTestCase(const std::string name); + + /** Destructor. */ + ~MatrixArrayTestCase() override; + /** + * \brief Copy constructor. + * Instruct the compiler to generate the implicitly declared copy constructor + */ + MatrixArrayTestCase(const MatrixArrayTestCase&) = default; + /** + * \brief Copy assignment operator. + * Instruct the compiler to generate the implicitly declared copy assignment operator. + * \return A reference to this MatrixArrayTestCase + */ + MatrixArrayTestCase& operator=(const MatrixArrayTestCase&) = default; + /** + * \brief Move constructor. + * Instruct the compiler to generate the implicitly declared move constructor + */ + MatrixArrayTestCase(MatrixArrayTestCase&&) noexcept = default; + /** + * \brief Move assignment operator. + * Instruct the compiler to generate the implicitly declared copy constructor + * \return A reference to this MatrixArrayTestCase + */ + MatrixArrayTestCase& operator=(MatrixArrayTestCase&&) noexcept = default; + + protected: + private: + void DoRun() override; +}; + +template +MatrixArrayTestCase::MatrixArrayTestCase(const std::string name) + : TestCase(name) +{ +} + +template +MatrixArrayTestCase::~MatrixArrayTestCase() +{ +} + +template +void +MatrixArrayTestCase::DoRun() +{ + // test multiplication of matrices (MatrixArray containing only 1 matrix) + MatrixArray m1 = MatrixArray(2, 3); + MatrixArray m2 = MatrixArray(m1.GetNumCols(), m1.GetNumRows()); + for (auto i = 0; i < m1.GetNumRows(); ++i) + { + for (auto j = 0; j < m1.GetNumCols(); ++j) + { + m1(i, j) = 1; + m2(j, i) = 1; + } + } + MatrixArray m3 = m1 * m2; + NS_LOG_INFO("m1:" << m1); + NS_LOG_INFO("m2:" << m2); + NS_LOG_INFO("m3 = m1 * m2:" << m3); + NS_TEST_ASSERT_MSG_EQ(m3.GetNumRows(), + m1.GetNumRows(), + "The number of rows in resulting matrix is not correct"); + NS_TEST_ASSERT_MSG_EQ(m3.GetNumCols(), + m2.GetNumCols(), + "The number of cols in resulting matrix is not correct"); + NS_TEST_ASSERT_MSG_EQ(m3.GetNumRows(), + m3.GetNumCols(), + "The number of rows and cols should be equal"); + for (auto i = 0; i < m3.GetNumCols(); ++i) + { + for (auto j = 0; j < m3.GetNumRows(); ++j) + { + NS_TEST_ASSERT_MSG_EQ(std::real(m3(i, j)), + m1.GetNumCols(), + "The element value should be " << m1.GetNumCols()); + } + } + + // multiplication with a scalar value + MatrixArray m4 = m3 * (static_cast(5.0)); + for (auto i = 0; i < m4.GetNumCols(); ++i) + { + for (auto j = 0; j < m4.GetNumRows(); ++j) + { + NS_TEST_ASSERT_MSG_EQ(m3(i, j) * (static_cast(5.0)), + m4(i, j), + "The values are not equal"); + } + } + NS_LOG_INFO("m4 = m3 * 5:" << m4); + + // test multiplication of arrays of matrices + MatrixArray m5 = MatrixArray(2, 3, 2); + MatrixArray m6 = MatrixArray(m5.GetNumCols(), m5.GetNumRows(), m5.GetNumPages()); + for (auto p = 0; p < m5.GetNumPages(); ++p) + { + for (auto i = 0; i < m5.GetNumRows(); ++i) + { + for (auto j = 0; j < m5.GetNumCols(); ++j) + { + m5(i, j, p) = 1; + m6(j, i, p) = 1; + } + } + } + MatrixArray m7 = m5 * m6; + NS_TEST_ASSERT_MSG_EQ(m7.GetNumRows(), + m5.GetNumRows(), + "The number of rows in resulting matrix is not correct"); + NS_TEST_ASSERT_MSG_EQ(m7.GetNumCols(), + m6.GetNumCols(), + "The number of cols in resulting matrix is not correct"); + NS_TEST_ASSERT_MSG_EQ(m7.GetNumRows(), + m7.GetNumCols(), + "The number of rows and cols should be equal"); + + for (auto p = 0; p < m7.GetNumPages(); ++p) + { + for (auto i = 0; i < m7.GetNumCols(); ++i) + { + for (auto j = 0; j < m7.GetNumRows(); ++j) + { + NS_TEST_ASSERT_MSG_EQ(std::real(m7(i, j, p)), + m5.GetNumCols(), + "The element value should be " << m5.GetNumCols()); + } + } + } + // test ostream operator + NS_LOG_INFO("m5:" << m5); + NS_LOG_INFO("m6:" << m6); + NS_LOG_INFO("m7 = m5 * m6:" << m7); + + // test transpose function + MatrixArray m8 = m5.Transpose(); + NS_TEST_ASSERT_MSG_EQ(m6, m8, "These two matrices should be equal"); + NS_LOG_INFO("m8 = m5.Transpose ()" << m8); + + // test 1D array creation, i.e. vector and transposing it + MatrixArray m9 = MatrixArray(std::vector({0, 1, 2, 3, 4, 5, 6, 7})); + NS_TEST_ASSERT_MSG_EQ((m9.GetNumRows() == 8) && (m9.GetNumCols() == 1) && + (m9.GetNumPages() == 1), + true, + "Creation of vector is not correct."); + + NS_LOG_INFO("Vector:" << m9); + NS_LOG_INFO("Vector after transposing:" << m9.Transpose()); + + // Test basic operators + MatrixArray m10 = + MatrixArray(m9.GetNumRows(), m9.GetNumCols(), m9.GetNumPages(), m9.GetValues()); + NS_TEST_ASSERT_MSG_EQ(m10, m9, "m10 and m9 should be equal"); + m10 -= m9; + NS_TEST_ASSERT_MSG_NE(m10, m9, "m10 and m9 should not be equal"); + m10 += m9; + NS_TEST_ASSERT_MSG_EQ(m10, m9, "m10 and m9 should be equal"); + m10 = m9; + NS_TEST_ASSERT_MSG_EQ(m10, m9, "m10 and m9 should be equal"); + m10 = m9 + m9; + NS_TEST_ASSERT_MSG_NE(m10, m9, "m10 and m9 should not be equal"); + m10 = m10 - m9; + NS_TEST_ASSERT_MSG_EQ(m10, m9, "m10 and m9 should be equal"); + + // test multiplication by using an initialization matrixArray + std::valarray a{0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5}; + std::valarray b{0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}; + std::valarray c{2, 3, 4, 6, 2, 3, 4, 6}; + std::valarray aCasted(a.size()); + std::valarray bCasted(b.size()); + std::valarray cCasted(c.size()); + + for (size_t i = 0; i < a.size(); ++i) + { + aCasted[i] = static_cast(a[i]); + } + for (size_t i = 0; i < b.size(); ++i) + { + bCasted[i] = static_cast(b[i]); + } + for (size_t i = 0; i < c.size(); ++i) + { + cCasted[i] = static_cast(c[i]); + } + MatrixArray m11 = MatrixArray(2, 3, 2, aCasted); + MatrixArray m12 = MatrixArray(3, 2, 2, bCasted); + MatrixArray m13 = m11 * m12; + MatrixArray m14 = MatrixArray(2, 2, 2, cCasted); + NS_TEST_ASSERT_MSG_EQ(m13.GetNumCols(), + m14.GetNumCols(), + "The number of columns is not as expected."); + NS_TEST_ASSERT_MSG_EQ(m13.GetNumRows(), + m14.GetNumRows(), + "The number of rows is not as expected."); + NS_TEST_ASSERT_MSG_EQ(m13, m14, "The values are not equal."); + NS_LOG_INFO("m11:" << m11); + NS_LOG_INFO("m12:" << m12); + NS_LOG_INFO("m13 = matrixArrayA * matrixArrayB:" << m13); + + // test MultiplyByLeftAndRightMatrix + std::valarray d{1, 1, 1}; + std::valarray e{1, 1}; + std::valarray f{1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3}; + std::valarray g{12, 12}; + std::valarray dCasted(d.size()); + std::valarray eCasted(e.size()); + std::valarray fCasted(f.size()); + std::valarray gCasted(g.size()); + for (size_t i = 0; i < d.size(); ++i) + { + dCasted[i] = static_cast(d[i]); + } + for (size_t i = 0; i < e.size(); ++i) + { + eCasted[i] = static_cast(e[i]); + } + for (size_t i = 0; i < f.size(); ++i) + { + fCasted[i] = static_cast(f[i]); + } + for (size_t i = 0; i < g.size(); ++i) + { + gCasted[i] = static_cast(g[i]); + } + MatrixArray m15 = MatrixArray(1, 3, dCasted); + MatrixArray m16 = MatrixArray(2, 1, eCasted); + MatrixArray m17 = MatrixArray(3, 2, 2, fCasted); + MatrixArray m18 = MatrixArray(1, 1, 2, gCasted); + MatrixArray m19 = m17.MultiplyByLeftAndRightMatrix(m15, m16); + NS_TEST_ASSERT_MSG_EQ(m19, m18, "The matrices should be equal."); + + // test MultiplyByLeftAndRightMatrix + std::valarray h{1, 3, 2, 2, 4, 0}; + std::valarray j{2, 2, 3, 4, 1, 3, 0, 5}; + std::valarray k{1, 2, 0, 0, 2, 3, 4, 1, 2, 3, 4, 1, 1, 2, 0, 0, 2, 3, 4, 1, 2, 3, 4, 1}; + std::valarray l{144, 132, 128, 104, 144, 132, 128, 104}; + std::valarray hCasted(h.size()); + std::valarray jCasted(j.size()); + std::valarray kCasted(k.size()); + std::valarray lCasted(l.size()); + for (size_t i = 0; i < h.size(); ++i) + { + hCasted[i] = static_cast(h[i]); + } + for (size_t i = 0; i < j.size(); ++i) + { + jCasted[i] = static_cast(j[i]); + } + for (size_t i = 0; i < k.size(); ++i) + { + kCasted[i] = static_cast(k[i]); + } + for (size_t i = 0; i < l.size(); ++i) + { + lCasted[i] = static_cast(l[i]); + } + MatrixArray m20 = MatrixArray(2, 3, hCasted); + MatrixArray m21 = MatrixArray(4, 2, jCasted); + MatrixArray m22 = MatrixArray(3, 4, 2, kCasted); + MatrixArray m23 = MatrixArray(2, 2, 2, lCasted); + MatrixArray m24 = m22.MultiplyByLeftAndRightMatrix(m20, m21); + NS_TEST_ASSERT_MSG_EQ(m24, m23, "The matrices should be equal."); + NS_LOG_INFO("m20:" << m20); + NS_LOG_INFO("m21:" << m21); + NS_LOG_INFO("m22:" << m22); + NS_LOG_INFO("m24 = m20 * m22 * m21" << m24); + + // test initialization with moving + NS_LOG_INFO("size() of lCasted before move: " << lCasted.size()); + MatrixArray m25 = MatrixArray(2, 2, 2, std::move(lCasted)); + NS_LOG_INFO("size() of lCasted after move: " << lCasted.size()); + NS_LOG_INFO("m25.GetSize ()" << m25.GetSize()); + NS_TEST_ASSERT_MSG_EQ(lCasted.size(), 0, "The size of lCasted should be 0."); + + NS_LOG_INFO("size() of hCasted before move: " << hCasted.size()); + MatrixArray m26 = MatrixArray(2, 3, std::move(hCasted)); + NS_LOG_INFO("size() of hCasted after move: " << hCasted.size()); + NS_LOG_INFO("m26.GetSize ()" << m26.GetSize()); + NS_TEST_ASSERT_MSG_EQ(hCasted.size(), 0, "The size of hCasted should be 0."); + + NS_LOG_INFO("size() of jCasted before move: " << jCasted.size()); + MatrixArray m27 = MatrixArray(std::move(jCasted)); + NS_LOG_INFO("size() of jCasted after move: " << jCasted.size()); + NS_LOG_INFO("m27.GetSize ()" << m27.GetSize()); + NS_TEST_ASSERT_MSG_EQ(hCasted.size(), 0, "The size of jCasted should be 0."); +} + +/** + * \ingroup matrixArray-tests + * MatrixArray test suite + */ +class MatrixArrayTestSuite : public TestSuite +{ + public: + /** Constructor. */ + MatrixArrayTestSuite(); +}; + +MatrixArrayTestSuite::MatrixArrayTestSuite() + : TestSuite("matrix-array-test") +{ + AddTestCase(new MatrixArrayTestCase("Test MatrixArray")); + AddTestCase( + new MatrixArrayTestCase>("Test MatrixArray>")); + AddTestCase(new MatrixArrayTestCase("Test MatrixArray")); +} + +/** + * \ingroup matrixArray-tests + * MatrixArrayTestSuite instance variable. + */ +static MatrixArrayTestSuite g_matrixArrayTestSuite; + +} // namespace tests +} // namespace ns3