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
This commit is contained in:
Biljana Bojovic
committed by
Biljana Bojovic
Biljana Bojovic
parent
abcee41e3d
commit
b7c3ba83fb
@@ -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
|
||||
|
||||
227
src/core/model/matrix-array.cc
Normal file
227
src/core/model/matrix-array.cc
Normal file
@@ -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 <bbojovic@cttc.es>
|
||||
*/
|
||||
|
||||
#include "matrix-array.h"
|
||||
|
||||
#ifdef HAVE_EIGEN3
|
||||
#include <Eigen/Dense>
|
||||
#endif
|
||||
|
||||
namespace ns3
|
||||
{
|
||||
|
||||
#ifdef HAVE_EIGEN3
|
||||
template <class T>
|
||||
using EigenMatrix = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
|
||||
template <class T>
|
||||
using ConstEigenMatrix = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
|
||||
#endif
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, uint16_t numPages)
|
||||
: ValArray<T>(numRows, numCols, numPages){};
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(const std::valarray<T>& values)
|
||||
: ValArray<T>(values)
|
||||
{
|
||||
}
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(std::valarray<T>&& values)
|
||||
: ValArray<T>(std::move(values))
|
||||
{
|
||||
}
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(const std::vector<T>& values)
|
||||
: ValArray<T>(values)
|
||||
{
|
||||
}
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, const std::valarray<T>& values)
|
||||
: ValArray<T>(numRows, numCols, values){};
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, std::valarray<T>&& values)
|
||||
: ValArray<T>(numRows, numCols, std::move(values)){};
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(uint16_t numRows,
|
||||
uint16_t numCols,
|
||||
uint16_t numPages,
|
||||
const std::valarray<T>& values)
|
||||
: ValArray<T>(numRows, numCols, numPages, values){};
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>::MatrixArray(uint16_t numRows,
|
||||
uint16_t numCols,
|
||||
uint16_t numPages,
|
||||
std::valarray<T>&& values)
|
||||
: ValArray<T>(numRows, numCols, numPages, std::move(values)){};
|
||||
|
||||
template <class T>
|
||||
MatrixArray<T>
|
||||
MatrixArray<T>::operator*(const MatrixArray<T>& 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<T> 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<T> lhsEigenMatrix(GetPagePtr(page), m_numRows, m_numCols);
|
||||
ConstEigenMatrix<T> rhsEigenMatrix(rhs.GetPagePtr(page), rhs.m_numRows, rhs.m_numCols);
|
||||
EigenMatrix<T> 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 <class T>
|
||||
MatrixArray<T>
|
||||
MatrixArray<T>::Transpose() const
|
||||
{
|
||||
MatrixArray<T> 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<T> thisMatrix(GetPagePtr(page), m_numRows, m_numCols);
|
||||
EigenMatrix<T> 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 <class T>
|
||||
MatrixArray<T>
|
||||
MatrixArray<T>::MultiplyByLeftAndRightMatrix(const MatrixArray<T>& lMatrix,
|
||||
const MatrixArray<T>& 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<T> res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages};
|
||||
|
||||
#ifdef HAVE_EIGEN3
|
||||
|
||||
ConstEigenMatrix<T> lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols);
|
||||
ConstEigenMatrix<T> 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<T> matrixEigen(GetPagePtr(page), m_numRows, m_numCols);
|
||||
EigenMatrix<T> 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<T> 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 <class T>
|
||||
template <bool EnableBool, typename>
|
||||
MatrixArray<T>
|
||||
MatrixArray<T>::HermitianTranspose() const
|
||||
{
|
||||
MatrixArray<std::complex<double>> retMatrix = this->Transpose();
|
||||
|
||||
for (size_t index = 0; index < this->GetSize(); ++index)
|
||||
{
|
||||
retMatrix[index] = std::conj(m_values[index]);
|
||||
}
|
||||
return retMatrix;
|
||||
}
|
||||
|
||||
template MatrixArray<std::complex<double>> MatrixArray<std::complex<double>>::HermitianTranspose()
|
||||
const;
|
||||
template class MatrixArray<std::complex<double>>;
|
||||
template class MatrixArray<double>;
|
||||
template class MatrixArray<int>;
|
||||
|
||||
} // namespace ns3
|
||||
308
src/core/model/matrix-array.h
Normal file
308
src/core/model/matrix-array.h
Normal file
@@ -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 <bbojovic@cttc.es>
|
||||
*/
|
||||
|
||||
#ifndef MATRIX_ARRAY_H
|
||||
#define MATRIX_ARRAY_H
|
||||
|
||||
#include "val-array.h"
|
||||
|
||||
#include <valarray>
|
||||
|
||||
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 T>
|
||||
class MatrixArray : public ValArray<T>
|
||||
{
|
||||
public:
|
||||
// instruct the compiler to generate the default constructor
|
||||
MatrixArray<T>() = 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<T>(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<T> values to initialize the elements.
|
||||
* \param values std::valarray<T> that will be used to initialize elements of 1D array
|
||||
*/
|
||||
explicit MatrixArray<T>(const std::valarray<T>& values);
|
||||
/**
|
||||
* \brief Constructor creates a single array of values.size () elements and 1 column,
|
||||
* and moves std::valarray<T> values to initialize the elements.
|
||||
* \param values std::valarray<T> that will be moved to initialize elements of 1D array
|
||||
*/
|
||||
MatrixArray<T>(std::valarray<T>&& values);
|
||||
/**
|
||||
* \brief Constructor creates a single array of values.size () elements and 1 column,
|
||||
* and uses std::vector<T> values to initialize the elements.
|
||||
* \param values std::vector<T> that will be used to initialize elements of 1D array
|
||||
*/
|
||||
explicit MatrixArray<T>(const std::vector<T>& values);
|
||||
/**
|
||||
* \brief Constructor creates a single matrix of numRows and numCols, and uses
|
||||
* std::valarray<T> values to initialize the elements.
|
||||
* \param numRows the number of rows
|
||||
* \param numCols the number of columns
|
||||
* \param values std::valarray<T> that will be used to initialize elements of 3D array
|
||||
*/
|
||||
MatrixArray<T>(uint16_t numRows, uint16_t numCols, const std::valarray<T>& values);
|
||||
/**
|
||||
* \brief Constructor creates a single matrix of numRows and numCols, and moves
|
||||
* std::valarray<T> values to initialize the elements.
|
||||
* \param numRows the number of rows
|
||||
* \param numCols the number of columns
|
||||
* \param values std::valarray<T> that will be moved to initialize elements of 3D array
|
||||
*/
|
||||
MatrixArray<T>(uint16_t numRows, uint16_t numCols, std::valarray<T>&& values);
|
||||
/**
|
||||
* \brief Constructor creates the array of numPages matrices of numRows x numCols dimensions,
|
||||
* and uses std::valarray<T> 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<T> that will be used to initialize elements of 3D array
|
||||
*/
|
||||
MatrixArray<T>(uint16_t numRows,
|
||||
uint16_t numCols,
|
||||
uint16_t numPages,
|
||||
const std::valarray<T>& values);
|
||||
/**
|
||||
* \brief Constructor creates the array of numPages matrices of numRows x numCols dimensions,
|
||||
* and moves std::valarray<T> 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<T> that will be used to initialize elements of 3D array
|
||||
*/
|
||||
MatrixArray<T>(uint16_t numRows,
|
||||
uint16_t numCols,
|
||||
uint16_t numPages,
|
||||
std::valarray<T>&& values);
|
||||
/** instruct the compiler to generate the implicitly declared destructor*/
|
||||
~MatrixArray<T>() override = default;
|
||||
/** instruct the compiler to generate the implicitly declared copy constructor*/
|
||||
MatrixArray<T>(const MatrixArray<T>&) = 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<T>& operator=(const MatrixArray<T>&) = default;
|
||||
/** instruct the compiler to generate the implicitly declared move constructor*/
|
||||
MatrixArray<T>(MatrixArray<T>&&) 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<T>& operator=(MatrixArray<T>&&) 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<T> operator*(const T& rhs) const;
|
||||
/**
|
||||
* \brief operator+ definition for MatrixArray<T>.
|
||||
* \param rhs The rhs MatrixArray object
|
||||
* \returns The result of operator+ of this MatrixArray and rhs MatrixAray
|
||||
*/
|
||||
MatrixArray<T> operator+(const MatrixArray<T>& rhs) const;
|
||||
/**
|
||||
* \brief binary operator- definition for MatrixArray<T>.
|
||||
* \param rhs The rhs MatrixArray object
|
||||
* \returns The result of operator- of this MatrixArray and rhs MatrixAray
|
||||
*/
|
||||
MatrixArray<T> operator-(const MatrixArray<T>& rhs) const;
|
||||
/**
|
||||
* \brief unary operator- definition for MatrixArray<T>.
|
||||
* \return the result of the operator-
|
||||
*/
|
||||
MatrixArray<T> 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<T> operator*(const MatrixArray<T>& 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<T>.
|
||||
* \return The resulting MatrixArray composed of the array of transposed matrices.
|
||||
*/
|
||||
MatrixArray<T> 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<T>
|
||||
* 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<T> MultiplyByLeftAndRightMatrix(const MatrixArray<T>& lMatrix,
|
||||
const MatrixArray<T>& rMatrix) const;
|
||||
|
||||
using ValArray<T>::GetPagePtr;
|
||||
using ValArray<T>::EqualDims;
|
||||
using ValArray<T>::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 <std::complex<double>> specialization of MatrixArray.
|
||||
* \return Returns a new matrix that is the result of the Hermitian transpose
|
||||
*/
|
||||
template <bool EnableBool = true,
|
||||
typename = typename std::enable_if<(std::is_same<T, std::complex<double>>::value &&
|
||||
EnableBool)>::type>
|
||||
MatrixArray<T> HermitianTranspose() const;
|
||||
|
||||
private:
|
||||
// To simplify functions in MatrixArray that are using members from the template base class
|
||||
using ValArray<T>::m_numRows;
|
||||
using ValArray<T>::m_numCols;
|
||||
using ValArray<T>::m_numPages;
|
||||
using ValArray<T>::m_values;
|
||||
};
|
||||
|
||||
using IntMatrixArray = MatrixArray<int>; //!< Create an alias for MatrixArray using int type
|
||||
using DoubleMatrixArray =
|
||||
MatrixArray<double>; //!< Create an alias for MatrixArray using double type
|
||||
using ComplexMatrixArray =
|
||||
MatrixArray<std::complex<double>>; //!< Create an alias for MatrixArray using complex type
|
||||
|
||||
/*************************************************
|
||||
** Class MatrixArray inline implementations
|
||||
************************************************/
|
||||
template <class T>
|
||||
inline MatrixArray<T>
|
||||
MatrixArray<T>::operator*(const T& rhs) const
|
||||
{
|
||||
return MatrixArray<T>(m_numRows,
|
||||
m_numCols,
|
||||
m_numPages,
|
||||
m_values * std::valarray<T>(rhs, m_numRows * m_numCols * m_numPages));
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline MatrixArray<T>
|
||||
MatrixArray<T>::operator+(const MatrixArray<T>& rhs) const
|
||||
{
|
||||
AssertEqualDims(rhs);
|
||||
return MatrixArray<T>(m_numRows, m_numCols, m_numPages, m_values + rhs.m_values);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline MatrixArray<T>
|
||||
MatrixArray<T>::operator-(const MatrixArray<T>& rhs) const
|
||||
{
|
||||
AssertEqualDims(rhs);
|
||||
return MatrixArray<T>(m_numRows, m_numCols, m_numPages, m_values - rhs.m_values);
|
||||
}
|
||||
|
||||
template <class T>
|
||||
inline MatrixArray<T>
|
||||
MatrixArray<T>::operator-() const
|
||||
{
|
||||
return MatrixArray<T>(m_numRows, m_numCols, m_numPages, -m_values);
|
||||
}
|
||||
|
||||
} // namespace ns3
|
||||
|
||||
#endif // MATRIX_ARRAY_H
|
||||
369
src/core/test/matrix-array-test-suite.cc
Normal file
369
src/core/test/matrix-array-test-suite.cc
Normal file
@@ -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 <bbojovic@cttc.es>
|
||||
*/
|
||||
|
||||
#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 T>
|
||||
class MatrixArrayTestCase : public TestCase
|
||||
{
|
||||
public:
|
||||
MatrixArrayTestCase<T>() = default;
|
||||
/**
|
||||
* Constructor
|
||||
*
|
||||
* \param [in] name reference name
|
||||
*/
|
||||
MatrixArrayTestCase<T>(const std::string name);
|
||||
|
||||
/** Destructor. */
|
||||
~MatrixArrayTestCase<T>() override;
|
||||
/**
|
||||
* \brief Copy constructor.
|
||||
* Instruct the compiler to generate the implicitly declared copy constructor
|
||||
*/
|
||||
MatrixArrayTestCase<T>(const MatrixArrayTestCase<T>&) = default;
|
||||
/**
|
||||
* \brief Copy assignment operator.
|
||||
* Instruct the compiler to generate the implicitly declared copy assignment operator.
|
||||
* \return A reference to this MatrixArrayTestCase
|
||||
*/
|
||||
MatrixArrayTestCase<T>& operator=(const MatrixArrayTestCase<T>&) = default;
|
||||
/**
|
||||
* \brief Move constructor.
|
||||
* Instruct the compiler to generate the implicitly declared move constructor
|
||||
*/
|
||||
MatrixArrayTestCase<T>(MatrixArrayTestCase<T>&&) noexcept = default;
|
||||
/**
|
||||
* \brief Move assignment operator.
|
||||
* Instruct the compiler to generate the implicitly declared copy constructor
|
||||
* \return A reference to this MatrixArrayTestCase
|
||||
*/
|
||||
MatrixArrayTestCase<T>& operator=(MatrixArrayTestCase<T>&&) noexcept = default;
|
||||
|
||||
protected:
|
||||
private:
|
||||
void DoRun() override;
|
||||
};
|
||||
|
||||
template <class T>
|
||||
MatrixArrayTestCase<T>::MatrixArrayTestCase(const std::string name)
|
||||
: TestCase(name)
|
||||
{
|
||||
}
|
||||
|
||||
template <class T>
|
||||
MatrixArrayTestCase<T>::~MatrixArrayTestCase<T>()
|
||||
{
|
||||
}
|
||||
|
||||
template <class T>
|
||||
void
|
||||
MatrixArrayTestCase<T>::DoRun()
|
||||
{
|
||||
// test multiplication of matrices (MatrixArray containing only 1 matrix)
|
||||
MatrixArray<T> m1 = MatrixArray<T>(2, 3);
|
||||
MatrixArray<T> m2 = MatrixArray<T>(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<T> 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<T> m4 = m3 * (static_cast<T>(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<T>(5.0)),
|
||||
m4(i, j),
|
||||
"The values are not equal");
|
||||
}
|
||||
}
|
||||
NS_LOG_INFO("m4 = m3 * 5:" << m4);
|
||||
|
||||
// test multiplication of arrays of matrices
|
||||
MatrixArray<T> m5 = MatrixArray<T>(2, 3, 2);
|
||||
MatrixArray<T> m6 = MatrixArray<T>(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<T> 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<T> 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<T> m9 = MatrixArray<T>(std::vector<T>({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<T> m10 =
|
||||
MatrixArray<T>(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<int> a{0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5};
|
||||
std::valarray<int> b{0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
|
||||
std::valarray<int> c{2, 3, 4, 6, 2, 3, 4, 6};
|
||||
std::valarray<T> aCasted(a.size());
|
||||
std::valarray<T> bCasted(b.size());
|
||||
std::valarray<T> cCasted(c.size());
|
||||
|
||||
for (size_t i = 0; i < a.size(); ++i)
|
||||
{
|
||||
aCasted[i] = static_cast<T>(a[i]);
|
||||
}
|
||||
for (size_t i = 0; i < b.size(); ++i)
|
||||
{
|
||||
bCasted[i] = static_cast<T>(b[i]);
|
||||
}
|
||||
for (size_t i = 0; i < c.size(); ++i)
|
||||
{
|
||||
cCasted[i] = static_cast<T>(c[i]);
|
||||
}
|
||||
MatrixArray<T> m11 = MatrixArray<T>(2, 3, 2, aCasted);
|
||||
MatrixArray<T> m12 = MatrixArray<T>(3, 2, 2, bCasted);
|
||||
MatrixArray<T> m13 = m11 * m12;
|
||||
MatrixArray<T> m14 = MatrixArray<T>(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<int> d{1, 1, 1};
|
||||
std::valarray<int> e{1, 1};
|
||||
std::valarray<int> f{1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3};
|
||||
std::valarray<int> g{12, 12};
|
||||
std::valarray<T> dCasted(d.size());
|
||||
std::valarray<T> eCasted(e.size());
|
||||
std::valarray<T> fCasted(f.size());
|
||||
std::valarray<T> gCasted(g.size());
|
||||
for (size_t i = 0; i < d.size(); ++i)
|
||||
{
|
||||
dCasted[i] = static_cast<T>(d[i]);
|
||||
}
|
||||
for (size_t i = 0; i < e.size(); ++i)
|
||||
{
|
||||
eCasted[i] = static_cast<T>(e[i]);
|
||||
}
|
||||
for (size_t i = 0; i < f.size(); ++i)
|
||||
{
|
||||
fCasted[i] = static_cast<T>(f[i]);
|
||||
}
|
||||
for (size_t i = 0; i < g.size(); ++i)
|
||||
{
|
||||
gCasted[i] = static_cast<T>(g[i]);
|
||||
}
|
||||
MatrixArray<T> m15 = MatrixArray<T>(1, 3, dCasted);
|
||||
MatrixArray<T> m16 = MatrixArray<T>(2, 1, eCasted);
|
||||
MatrixArray<T> m17 = MatrixArray<T>(3, 2, 2, fCasted);
|
||||
MatrixArray<T> m18 = MatrixArray<T>(1, 1, 2, gCasted);
|
||||
MatrixArray<T> m19 = m17.MultiplyByLeftAndRightMatrix(m15, m16);
|
||||
NS_TEST_ASSERT_MSG_EQ(m19, m18, "The matrices should be equal.");
|
||||
|
||||
// test MultiplyByLeftAndRightMatrix
|
||||
std::valarray<int> h{1, 3, 2, 2, 4, 0};
|
||||
std::valarray<int> j{2, 2, 3, 4, 1, 3, 0, 5};
|
||||
std::valarray<int> 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<int> l{144, 132, 128, 104, 144, 132, 128, 104};
|
||||
std::valarray<T> hCasted(h.size());
|
||||
std::valarray<T> jCasted(j.size());
|
||||
std::valarray<T> kCasted(k.size());
|
||||
std::valarray<T> lCasted(l.size());
|
||||
for (size_t i = 0; i < h.size(); ++i)
|
||||
{
|
||||
hCasted[i] = static_cast<T>(h[i]);
|
||||
}
|
||||
for (size_t i = 0; i < j.size(); ++i)
|
||||
{
|
||||
jCasted[i] = static_cast<T>(j[i]);
|
||||
}
|
||||
for (size_t i = 0; i < k.size(); ++i)
|
||||
{
|
||||
kCasted[i] = static_cast<T>(k[i]);
|
||||
}
|
||||
for (size_t i = 0; i < l.size(); ++i)
|
||||
{
|
||||
lCasted[i] = static_cast<T>(l[i]);
|
||||
}
|
||||
MatrixArray<T> m20 = MatrixArray<T>(2, 3, hCasted);
|
||||
MatrixArray<T> m21 = MatrixArray<T>(4, 2, jCasted);
|
||||
MatrixArray<T> m22 = MatrixArray<T>(3, 4, 2, kCasted);
|
||||
MatrixArray<T> m23 = MatrixArray<T>(2, 2, 2, lCasted);
|
||||
MatrixArray<T> 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<T> m25 = MatrixArray<T>(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<T> m26 = MatrixArray<T>(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<T> m27 = MatrixArray<T>(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<double>("Test MatrixArray<double>"));
|
||||
AddTestCase(
|
||||
new MatrixArrayTestCase<std::complex<double>>("Test MatrixArray<std::complex<double>>"));
|
||||
AddTestCase(new MatrixArrayTestCase<int>("Test MatrixArray<int>"));
|
||||
}
|
||||
|
||||
/**
|
||||
* \ingroup matrixArray-tests
|
||||
* MatrixArrayTestSuite instance variable.
|
||||
*/
|
||||
static MatrixArrayTestSuite g_matrixArrayTestSuite;
|
||||
|
||||
} // namespace tests
|
||||
} // namespace ns3
|
||||
Reference in New Issue
Block a user