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 * SPDX-License-Identifier: GPL-2.0-only
5 *
6 * Author: Biljana Bojovic <bbojovic@cttc.es>
7 */
8
9#include "matrix-array.h"
10
11#ifdef HAVE_EIGEN3
12
13#if defined(__GNUC__) && !defined(__clang__)
14#pragma GCC diagnostic push
15#pragma GCC diagnostic ignored "-Wclass-memaccess"
16#pragma GCC diagnostic ignored "-Wunused-variable"
17#endif
18#if defined(__clang__)
19#pragma clang diagnostic push
20#pragma clang diagnostic ignored "-Wdeprecated-declarations"
21#endif
22
23#include <Eigen/Dense>
24
25#if defined(__clang__)
26#pragma clang diagnostic pop
27#endif
28#if defined(__GNUC__) && !defined(__clang__)
29#pragma GCC diagnostic pop
30#endif
31#endif
32
33namespace ns3
34{
35
36#ifdef HAVE_EIGEN3
37template <class T>
38using EigenMatrix = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
39template <class T>
40using ConstEigenMatrix = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
41#endif
42
43template <class T>
44MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, size_t numPages)
45 : ValArray<T>(numRows, numCols, numPages)
46{
47}
48
49template <class T>
50MatrixArray<T>::MatrixArray(const std::valarray<T>& values)
51 : ValArray<T>(values)
52{
53}
54
55template <class T>
56MatrixArray<T>::MatrixArray(std::valarray<T>&& values)
57 : ValArray<T>(std::move(values))
58{
59}
60
61template <class T>
62MatrixArray<T>::MatrixArray(const std::vector<T>& values)
63 : ValArray<T>(values)
64{
65}
66
67template <class T>
68MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, const std::valarray<T>& values)
69 : ValArray<T>(numRows, numCols, values)
70{
71}
72
73template <class T>
74MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, std::valarray<T>&& values)
75 : ValArray<T>(numRows, numCols, std::move(values))
76{
77}
78
79template <class T>
81 size_t numCols,
82 size_t numPages,
83 const std::valarray<T>& values)
84 : ValArray<T>(numRows, numCols, numPages, values)
85{
86}
87
88template <class T>
90 size_t numCols,
91 size_t numPages,
92 std::valarray<T>&& values)
93 : ValArray<T>(numRows, numCols, numPages, std::move(values))
94{
95}
97template <class T>
100{
101 NS_ASSERT_MSG(m_numPages == rhs.m_numPages, "MatrixArrays have different numbers of matrices.");
102 NS_ASSERT_MSG(m_numCols == rhs.m_numRows, "Inner dimensions of matrices mismatch.");
103
104 MatrixArray<T> res{m_numRows, rhs.m_numCols, m_numPages};
105
106 for (size_t page = 0; page < res.m_numPages; ++page)
107 {
108#ifdef HAVE_EIGEN3 // Eigen found and enabled Eigen optimizations
109
110 ConstEigenMatrix<T> lhsEigenMatrix(GetPagePtr(page), m_numRows, m_numCols);
111 ConstEigenMatrix<T> rhsEigenMatrix(rhs.GetPagePtr(page), rhs.m_numRows, rhs.m_numCols);
112 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
113 resEigenMatrix = lhsEigenMatrix * rhsEigenMatrix;
114
115#else // Eigen not found or Eigen optimizations not enabled
116
117 size_t matrixOffset = page * m_numRows * m_numCols;
118 size_t rhsMatrixOffset = page * rhs.m_numRows * rhs.m_numCols;
119 for (size_t i = 0; i < res.m_numRows; ++i)
120 {
121 for (size_t j = 0; j < res.m_numCols; ++j)
122 {
123 res(i, j, page) = (m_values[std::slice(matrixOffset + i, m_numCols, m_numRows)] *
124 rhs.m_values[std::slice(rhsMatrixOffset + j * rhs.m_numRows,
125 rhs.m_numRows,
126 1)])
127 .sum();
128 }
129 }
130
131#endif
132 }
133 return res;
134}
135
136template <class T>
139{
140 // Create the matrix where m_numRows = this.m_numCols, m_numCols = this.m_numRows,
141 // m_numPages = this.m_numPages
142 MatrixArray<T> res{m_numCols, m_numRows, m_numPages};
143
144 for (size_t page = 0; page < m_numPages; ++page)
145 {
146#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
147
148 ConstEigenMatrix<T> thisMatrix(GetPagePtr(page), m_numRows, m_numCols);
149 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
150 resEigenMatrix = thisMatrix.transpose();
151
152#else // Eigen not found or Eigen optimizations not enabled
153
154 size_t matrixIndex = page * m_numRows * m_numCols;
155 for (size_t i = 0; i < m_numRows; ++i)
156 {
157 res.m_values[std::slice(matrixIndex + i * res.m_numRows, res.m_numRows, 1)] =
158 m_values[std::slice(matrixIndex + i, m_numCols, m_numRows)];
159 }
160
161#endif
162 }
163 return res;
164}
165
166template <class T>
169{
170 MatrixArray<T> res{1, 1, m_numPages};
171 NS_ASSERT_MSG(m_numRows == m_numCols, "Matrix is not square");
172 // In case of small matrices, we use a fast path
173 if (m_numRows == 1)
174 {
175 return *this;
176 }
177 // Calculate determinant for each matrix
178 for (size_t page = 0; page < m_numPages; ++page)
179 {
180 res(0, 0, page) = 0;
181 auto pageValues = GetPagePtr(page);
182
183 // Fast path for 2x2 matrices
184 if (m_numRows == 2)
186 res(0, 0, page) = pageValues[0] * pageValues[3] - pageValues[1] * pageValues[2];
187 continue;
188 }
189 for (size_t detN = 0; detN < m_numRows; detN++)
191 auto partDetP = T{0} + 1.0;
192 auto partDetN = T{0} + 1.0;
193 for (size_t row = 0; row < m_numRows; row++)
194 {
195 // Wraparound not to have to extend the matrix
196 // Positive determinant
197 size_t col = (row + detN) % m_numCols;
198 partDetP *= pageValues[row * m_numCols + col];
199
200 // Negative determinant
201 col = m_numCols - 1 - (row + detN) % m_numCols;
202 partDetN *= pageValues[row * m_numCols + col];
203 }
204 res(0, 0, page) += partDetP - partDetN;
205 }
206 }
207 return res;
208}
209
210template <class T>
211MatrixArray<T>
214 MatrixArray<T> res{1, 1, m_numPages};
215 for (size_t page = 0; page < m_numPages; ++page)
216 {
217 // Calculate the sum of squared absolute values of each matrix page
218 res[page] = 0;
219 auto pagePtr = this->GetPagePtr(page);
220 for (size_t i = 0; i < m_numRows * m_numCols; i++)
221 {
222 auto absVal = std::abs(pagePtr[i]);
223 res[page] += absVal * absVal;
224 }
225 res[page] = sqrt(res[page]);
226 }
227 return res;
228}
230template <class T>
233 const MatrixArray<T>& rMatrix) const
234{
235 NS_ASSERT_MSG(lMatrix.m_numPages == 1 && rMatrix.m_numPages == 1,
236 "The left and right MatrixArray should have only one page.");
237 NS_ASSERT_MSG(lMatrix.m_numCols == m_numRows,
238 "Left vector numCols and this MatrixArray numRows mismatch.");
239 NS_ASSERT_MSG(m_numCols == rMatrix.m_numRows,
240 "Right vector numRows and this MatrixArray numCols mismatch.");
241
242 MatrixArray<T> res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages};
244#ifdef HAVE_EIGEN3
245
246 ConstEigenMatrix<T> lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols);
247 ConstEigenMatrix<T> rMatrixEigen(rMatrix.GetPagePtr(0), rMatrix.m_numRows, rMatrix.m_numCols);
248#endif
250 for (size_t page = 0; page < m_numPages; ++page)
251 {
252#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
253
254 ConstEigenMatrix<T> matrixEigen(GetPagePtr(page), m_numRows, m_numCols);
255 EigenMatrix<T> resEigenMap(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
257 resEigenMap = lMatrixEigen * matrixEigen * rMatrixEigen;
258
259#else // Eigen not found or Eigen optimizations not enabled
260
261 size_t matrixOffset = page * m_numRows * m_numCols;
262 for (size_t resRow = 0; resRow < res.m_numRows; ++resRow)
264 for (size_t resCol = 0; resCol < res.m_numCols; ++resCol)
265 {
266 // create intermediate row result, a multiply of resRow row of lMatrix and each
267 // column of this matrix
268 std::valarray<T> interRes(m_numCols);
269 for (size_t thisCol = 0; thisCol < m_numCols; ++thisCol)
270 {
271 interRes[thisCol] =
272 (lMatrix
273 .m_values[std::slice(resRow, lMatrix.m_numCols, lMatrix.m_numRows)] *
274 m_values[std::slice(matrixOffset + thisCol * m_numRows, m_numRows, 1)])
275 .sum();
276 }
277 // multiply intermediate results and resCol column of the rMatrix
278 res(resRow, resCol, page) =
279 (interRes *
280 rMatrix.m_values[std::slice(resCol * rMatrix.m_numRows, rMatrix.m_numRows, 1)])
281 .sum();
282 }
283 }
284#endif
285 }
286 return res;
287}
288
289template <class T>
290template <bool EnableBool, typename>
291MatrixArray<T>
293{
294 MatrixArray<std::complex<double>> retMatrix = this->Transpose();
295
296 for (size_t index = 0; index < this->GetSize(); ++index)
297 {
298 retMatrix.m_values[index] = std::conj(retMatrix.m_values[index]);
299 }
300 return retMatrix;
301}
302
303template <class T>
305MatrixArray<T>::MakeNCopies(size_t nCopies) const
306{
307 NS_ASSERT_MSG(m_numPages == 1, "The MatrixArray should have only one page to be copied.");
308 auto copiedMatrix = MatrixArray<T>{m_numRows, m_numCols, nCopies};
309 for (size_t copy = 0; copy < nCopies; copy++)
310 {
311 for (size_t i = 0; i < m_numRows * m_numCols; i++)
312 {
313 copiedMatrix.GetPagePtr(copy)[i] = m_values[i];
314 }
315 }
316 return copiedMatrix;
317}
318
319template <class T>
322{
323 NS_ASSERT_MSG(page < m_numPages, "The page to extract from the MatrixArray is out of bounds.");
324 auto extractedPage = MatrixArray<T>{m_numRows, m_numCols, 1};
325
326 for (size_t i = 0; i < m_numRows * m_numCols; ++i)
327 {
328 extractedPage.m_values[i] = GetPagePtr(page)[i];
329 }
330 return extractedPage;
331}
332
333template <class T>
336{
337 auto jointMatrix =
338 MatrixArray<T>{pages.front().GetNumRows(), pages.front().GetNumCols(), pages.size()};
339 for (size_t page = 0; page < jointMatrix.GetNumPages(); page++)
340 {
341 NS_ASSERT_MSG(pages[page].GetNumRows() == jointMatrix.GetNumRows(),
342 "All page matrices should have the same number of rows");
343 NS_ASSERT_MSG(pages[page].GetNumCols() == jointMatrix.GetNumCols(),
344 "All page matrices should have the same number of columns");
345 NS_ASSERT_MSG(pages[page].GetNumPages() == 1,
346 "All page matrices should have a single page");
347
348 size_t i = 0;
349 for (auto a : pages[page].GetValues())
350 {
351 jointMatrix.GetPagePtr(page)[i] = a;
352 i++;
353 }
354 }
355 return jointMatrix;
356}
357
358template <class T>
360MatrixArray<T>::IdentityMatrix(const size_t size, const size_t pages)
361{
362 auto identityMatrix = MatrixArray<T>{size, size, pages};
363 for (std::size_t page = 0; page < pages; page++)
364 {
365 for (std::size_t i = 0; i < size; i++)
366 {
367 identityMatrix(i, i, page) = 1.0;
368 }
369 }
370 return identityMatrix;
371}
372
373template <class T>
376{
377 NS_ASSERT_MSG(likeme.GetNumRows() == likeme.GetNumCols(), "Template array is not square.");
378 return IdentityMatrix(likeme.GetNumRows(), likeme.GetNumPages());
379}
380
382 const;
384template class MatrixArray<double>;
385template class MatrixArray<int>;
386
387} // namespace ns3
MatrixArray class inherits ValArray class and provides additional interfaces to ValArray which enable...
MatrixArray< T > HermitianTranspose() const
Function that performs the Hermitian transpose of this MatrixArray and returns a new matrix that is t...
MatrixArray operator*(const T &rhs) const
Element-wise multiplication with a scalar value.
MatrixArray Determinant() const
This operator calculates a vector o determinants, one for each page.
static MatrixArray< T > IdentityMatrix(const size_t size, const size_t pages=1)
Function produces an identity MatrixArray with the specified size.
static MatrixArray< T > JoinPages(const std::vector< MatrixArray< T > > &pages)
Function joins multiple pages into a single MatrixArray.
MatrixArray Transpose() const
This operator interprets the 3D array as an array of matrices, and performs a linear algebra operatio...
MatrixArray< T > MakeNCopies(size_t nCopies) const
Function that copies the current 1-page matrix into a new matrix with n copies of the original matrix...
MatrixArray()=default
MatrixArray< T > ExtractPage(size_t page) const
Function extracts a page from a MatrixArray.
MatrixArray MultiplyByLeftAndRightMatrix(const MatrixArray< T > &lMatrix, const MatrixArray< T > &rMatrix) const
Multiply each matrix in the array by the left and the right matrix.
MatrixArray FrobeniusNorm() const
This operator calculates a vector of Frobenius norm, one for each page.
ValArray is a class to efficiently store 3D array.
Definition val-array.h:74
T * GetPagePtr(size_t pageIndex)
Get a data pointer to a specific 2D array for use in linear algebra libraries.
Definition val-array.h:516
size_t GetNumPages() const
Definition val-array.h:387
size_t m_numCols
The size of the second dimension, i.e., the number of columns of each 2D array.
Definition val-array.h:361
std::valarray< T > m_values
The data values.
Definition val-array.h:364
size_t GetNumRows() const
Definition val-array.h:373
size_t m_numRows
The size of the first dimension, i.e., the number of rows of each 2D array.
Definition val-array.h:359
size_t GetNumCols() const
Definition val-array.h:380
size_t m_numPages
The size of the third dimension, i.e., the number of 2D arrays.
Definition val-array.h:363
#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:75
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.
STL namespace.