A Discrete-Event Network Simulator
API
Loading...
Searching...
No Matches
matrix-array.cc
Go to the documentation of this file.
1/*
2 * Copyright (c) 2022 Centre Tecnologic de Telecomunicacions de Catalunya (CTTC)
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License version 2 as
6 * published by the Free Software Foundation;
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program; if not, write to the Free Software
15 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 *
17 * Author: Biljana Bojovic <bbojovic@cttc.es>
18 */
19
20#include "matrix-array.h"
21
22#ifdef HAVE_EIGEN3
23#include <Eigen/Dense>
24#endif
25
26namespace ns3
27{
28
29#ifdef HAVE_EIGEN3
30template <class T>
31using EigenMatrix = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
32template <class T>
33using ConstEigenMatrix = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
34#endif
35
36template <class T>
37MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, uint16_t numPages)
38 : ValArray<T>(numRows, numCols, numPages)
39{
40}
41
42template <class T>
43MatrixArray<T>::MatrixArray(const std::valarray<T>& values)
44 : ValArray<T>(values)
45{
46}
47
48template <class T>
49MatrixArray<T>::MatrixArray(std::valarray<T>&& values)
50 : ValArray<T>(std::move(values))
51{
52}
53
54template <class T>
55MatrixArray<T>::MatrixArray(const std::vector<T>& values)
56 : ValArray<T>(values)
57{
58}
59
60template <class T>
61MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, const std::valarray<T>& values)
62 : ValArray<T>(numRows, numCols, values)
63{
64}
65
66template <class T>
67MatrixArray<T>::MatrixArray(uint16_t numRows, uint16_t numCols, std::valarray<T>&& values)
68 : ValArray<T>(numRows, numCols, std::move(values))
69{
70}
71
72template <class T>
74 uint16_t numCols,
75 uint16_t numPages,
76 const std::valarray<T>& values)
77 : ValArray<T>(numRows, numCols, numPages, values)
78{
79}
80
81template <class T>
83 uint16_t numCols,
84 uint16_t numPages,
85 std::valarray<T>&& values)
86 : ValArray<T>(numRows, numCols, numPages, std::move(values))
87{
88}
89
90template <class T>
93{
94 NS_ASSERT_MSG(m_numPages == rhs.m_numPages, "MatrixArrays have different numbers of matrices.");
95 NS_ASSERT_MSG(m_numCols == rhs.m_numRows, "Inner dimensions of matrices mismatch.");
96
97 MatrixArray<T> res{m_numRows, rhs.m_numCols, m_numPages};
98
99 for (uint16_t page = 0; page < res.m_numPages; ++page)
100 {
101#ifdef HAVE_EIGEN3 // Eigen found and enabled Eigen optimizations
102
103 ConstEigenMatrix<T> lhsEigenMatrix(GetPagePtr(page), m_numRows, m_numCols);
104 ConstEigenMatrix<T> rhsEigenMatrix(rhs.GetPagePtr(page), rhs.m_numRows, rhs.m_numCols);
105 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
106 resEigenMatrix = lhsEigenMatrix * rhsEigenMatrix;
107
108#else // Eigen not found or Eigen optimizations not enabled
109
110 size_t matrixOffset = page * m_numRows * m_numCols;
111 for (uint16_t i = 0; i < res.m_numRows; ++i)
112 {
113 for (uint16_t j = 0; j < res.m_numCols; ++j)
114 {
115 res(i, j, page) =
116 (m_values[std::slice(matrixOffset + i, m_numCols, m_numRows)] *
117 rhs.m_values[std::slice(matrixOffset + j * rhs.m_numRows, rhs.m_numRows, 1)])
118 .sum();
119 }
120 }
121
122#endif
123 }
124 return res;
125}
126
127template <class T>
130{
131 // Create the matrix where m_numRows = this.m_numCols, m_numCols = this.m_numRows,
132 // m_numPages = this.m_numPages
133 MatrixArray<T> res{m_numCols, m_numRows, m_numPages};
134
135 for (uint16_t page = 0; page < m_numPages; ++page)
136 {
137#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
138
139 ConstEigenMatrix<T> thisMatrix(GetPagePtr(page), m_numRows, m_numCols);
140 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
141 resEigenMatrix = thisMatrix.transpose();
142
143#else // Eigen not found or Eigen optimizations not enabled
144
145 size_t matrixIndex = page * m_numRows * m_numCols;
146 for (uint16_t i = 0; i < m_numRows; ++i)
147 {
148 res.m_values[std::slice(matrixIndex + i * res.m_numRows, res.m_numRows, 1)] =
149 m_values[std::slice(matrixIndex + i, m_numCols, m_numRows)];
150 }
151
152#endif
153 }
154 return res;
155}
156
157template <class T>
160 const MatrixArray<T>& rMatrix) const
161{
162 NS_ASSERT_MSG(lMatrix.m_numPages == 1 && rMatrix.m_numPages == 1,
163 "The left and right MatrixArray should have only one page.");
164 NS_ASSERT_MSG(lMatrix.m_numCols == m_numRows,
165 "Left vector numCols and this MatrixArray numRows mismatch.");
166 NS_ASSERT_MSG(m_numCols == rMatrix.m_numRows,
167 "Right vector numRows and this MatrixArray numCols mismatch.");
168
169 MatrixArray<T> res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages};
170
171#ifdef HAVE_EIGEN3
172
173 ConstEigenMatrix<T> lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols);
174 ConstEigenMatrix<T> rMatrixEigen(rMatrix.GetPagePtr(0), rMatrix.m_numRows, rMatrix.m_numCols);
175#endif
176
177 for (uint16_t page = 0; page < m_numPages; ++page)
178 {
179#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
180
181 ConstEigenMatrix<T> matrixEigen(GetPagePtr(page), m_numRows, m_numCols);
182 EigenMatrix<T> resEigenMap(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
183
184 resEigenMap = lMatrixEigen * matrixEigen * rMatrixEigen;
185
186#else // Eigen not found or Eigen optimizations not enabled
187
188 size_t matrixOffset = page * m_numRows * m_numCols;
189 for (uint16_t resRow = 0; resRow < res.m_numRows; ++resRow)
190 {
191 for (uint16_t resCol = 0; resCol < res.m_numCols; ++resCol)
192 {
193 // create intermediate row result, a multiply of resRow row of lMatrix and each
194 // column of this matrix
195 std::valarray<T> interRes(m_numCols);
196 for (uint16_t thisCol = 0; thisCol < m_numCols; ++thisCol)
197 {
198 interRes[thisCol] =
199 (lMatrix
200 .m_values[std::slice(resRow, lMatrix.m_numCols, lMatrix.m_numRows)] *
201 m_values[std::slice(matrixOffset + thisCol * m_numRows, m_numRows, 1)])
202 .sum();
203 }
204 // multiply intermediate results and resCol column of the rMatrix
205 res(resRow, resCol, page) =
206 (interRes *
207 rMatrix.m_values[std::slice(resCol * rMatrix.m_numRows, rMatrix.m_numRows, 1)])
208 .sum();
209 }
210 }
211#endif
212 }
213 return res;
214}
215
216template <class T>
217template <bool EnableBool, typename>
220{
221 MatrixArray<std::complex<double>> retMatrix = this->Transpose();
222
223 for (size_t index = 0; index < this->GetSize(); ++index)
224 {
225 retMatrix[index] = std::conj(m_values[index]);
226 }
227 return retMatrix;
228}
229
231 const;
233template class MatrixArray<double>;
234template class MatrixArray<int>;
235
236} // namespace ns3
MatrixArray class inherits ValArray class and provides additional interfaces to ValArray which enable...
Definition: matrix-array.h:83
MatrixArray< T > HermitianTranspose() const
Function that performs the Hermitian transpose of this MatrixArray and returns a new matrix that is t...
MatrixArray< T > operator*(const T &rhs) const
Element-wise multiplication with a scalar value.
Definition: matrix-array.h:278
MatrixArray< T > Transpose() const
This operator interprets the 3D array as an array of matrices, and performs a linear algebra operatio...
MatrixArray()=default
MatrixArray< T > MultiplyByLeftAndRightMatrix(const MatrixArray< T > &lMatrix, const MatrixArray< T > &rMatrix) const
Multiply each matrix in the array by the left and the right matrix.
ValArray is a class to efficiently store 3D array.
Definition: val-array.h:80
uint16_t m_numRows
The size of the first dimension, i.e., the number of rows of each 2D array.
Definition: val-array.h:368
T * GetPagePtr(uint16_t pageIndex)
Get a data pointer to a specific 2D array for use in linear algebra libraries.
Definition: val-array.h:525
uint16_t m_numPages
The size of the third dimension, i.e., the number of 2D arrays.
Definition: val-array.h:372
std::valarray< T > m_values
The data values.
Definition: val-array.h:373
uint16_t m_numCols
The size of the second dimension, i.e., the number of columns of each 2D array.
Definition: val-array.h:370
#define NS_ASSERT_MSG(condition, message)
At runtime, in debugging builds, if this condition is not true, the program prints the message to out...
Definition: assert.h:86
Every class exported by the ns3 library is enclosed in the ns3 namespace.
uint32_t GetSize(Ptr< const Packet > packet, const WifiMacHeader *hdr, bool isAmpdu)
Return the total size of the packet after WifiMacHeader and FCS trailer have been added.
Definition: wifi-utils.cc:132
STL namespace.