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#endif
17
18#include <Eigen/Dense>
19
20#if defined(__GNUC__) && !defined(__clang__)
21#pragma GCC diagnostic pop
22#endif
23#endif
24
25namespace ns3
26{
27
28#ifdef HAVE_EIGEN3
29template <class T>
30using EigenMatrix = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
31template <class T>
32using ConstEigenMatrix = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
33#endif
34
35template <class T>
36MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, size_t numPages)
37 : ValArray<T>(numRows, numCols, numPages)
38{
39}
40
41template <class T>
42MatrixArray<T>::MatrixArray(const std::valarray<T>& values)
43 : ValArray<T>(values)
44{
45}
46
47template <class T>
48MatrixArray<T>::MatrixArray(std::valarray<T>&& values)
49 : ValArray<T>(std::move(values))
50{
51}
52
53template <class T>
54MatrixArray<T>::MatrixArray(const std::vector<T>& values)
55 : ValArray<T>(values)
56{
57}
58
59template <class T>
60MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, const std::valarray<T>& values)
61 : ValArray<T>(numRows, numCols, values)
62{
63}
64
65template <class T>
66MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, std::valarray<T>&& values)
67 : ValArray<T>(numRows, numCols, std::move(values))
68{
69}
70
71template <class T>
73 size_t numCols,
74 size_t numPages,
75 const std::valarray<T>& values)
76 : ValArray<T>(numRows, numCols, numPages, values)
77{
78}
79
80template <class T>
82 size_t numCols,
83 size_t numPages,
84 std::valarray<T>&& values)
85 : ValArray<T>(numRows, numCols, numPages, std::move(values))
86{
87}
88
89template <class T>
92{
93 NS_ASSERT_MSG(m_numPages == rhs.m_numPages, "MatrixArrays have different numbers of matrices.");
94 NS_ASSERT_MSG(m_numCols == rhs.m_numRows, "Inner dimensions of matrices mismatch.");
95
96 MatrixArray<T> res{m_numRows, rhs.m_numCols, m_numPages};
97
98 for (size_t page = 0; page < res.m_numPages; ++page)
99 {
100#ifdef HAVE_EIGEN3 // Eigen found and enabled Eigen optimizations
101
102 ConstEigenMatrix<T> lhsEigenMatrix(GetPagePtr(page), m_numRows, m_numCols);
103 ConstEigenMatrix<T> rhsEigenMatrix(rhs.GetPagePtr(page), rhs.m_numRows, rhs.m_numCols);
104 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
105 resEigenMatrix = lhsEigenMatrix * rhsEigenMatrix;
106
107#else // Eigen not found or Eigen optimizations not enabled
108
109 size_t matrixOffset = page * m_numRows * m_numCols;
110 size_t rhsMatrixOffset = page * rhs.m_numRows * rhs.m_numCols;
111 for (size_t i = 0; i < res.m_numRows; ++i)
112 {
113 for (size_t j = 0; j < res.m_numCols; ++j)
114 {
115 res(i, j, page) = (m_values[std::slice(matrixOffset + i, m_numCols, m_numRows)] *
116 rhs.m_values[std::slice(rhsMatrixOffset + j * rhs.m_numRows,
117 rhs.m_numRows,
118 1)])
119 .sum();
120 }
121 }
122
123#endif
124 }
125 return res;
126}
128template <class T>
131{
132 // Create the matrix where m_numRows = this.m_numCols, m_numCols = this.m_numRows,
133 // m_numPages = this.m_numPages
134 MatrixArray<T> res{m_numCols, m_numRows, m_numPages};
135
136 for (size_t page = 0; page < m_numPages; ++page)
137 {
138#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
139
140 ConstEigenMatrix<T> thisMatrix(GetPagePtr(page), m_numRows, m_numCols);
141 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
142 resEigenMatrix = thisMatrix.transpose();
143
144#else // Eigen not found or Eigen optimizations not enabled
145
146 size_t matrixIndex = page * m_numRows * m_numCols;
147 for (size_t i = 0; i < m_numRows; ++i)
148 {
149 res.m_values[std::slice(matrixIndex + i * res.m_numRows, res.m_numRows, 1)] =
150 m_values[std::slice(matrixIndex + i, m_numCols, m_numRows)];
151 }
152
153#endif
154 }
155 return res;
156}
157
158template <class T>
159MatrixArray<T>
161{
162 MatrixArray<T> res{1, 1, m_numPages};
163 NS_ASSERT_MSG(m_numRows == m_numCols, "Matrix is not square");
164 // In case of small matrices, we use a fast path
165 if (m_numRows == 1)
166 {
167 return *this;
168 }
169 // Calculate determinant for each matrix
170 for (size_t page = 0; page < m_numPages; ++page)
171 {
172 res(0, 0, page) = 0;
173 auto pageValues = GetPagePtr(page);
174
175 // Fast path for 2x2 matrices
176 if (m_numRows == 2)
177 {
178 res(0, 0, page) = pageValues[0] * pageValues[3] - pageValues[1] * pageValues[2];
179 continue;
181 for (size_t detN = 0; detN < m_numRows; detN++)
182 {
183 auto partDetP = T{0} + 1.0;
184 auto partDetN = T{0} + 1.0;
185 for (size_t row = 0; row < m_numRows; row++)
186 {
187 // Wraparound not to have to extend the matrix
188 // Positive determinant
189 size_t col = (row + detN) % m_numCols;
190 partDetP *= pageValues[row * m_numCols + col];
191
192 // Negative determinant
193 col = m_numCols - 1 - (row + detN) % m_numCols;
194 partDetN *= pageValues[row * m_numCols + col];
195 }
196 res(0, 0, page) += partDetP - partDetN;
197 }
198 }
199 return res;
200}
201
202template <class T>
203MatrixArray<T>
205{
206 MatrixArray<T> res{1, 1, m_numPages};
207 for (size_t page = 0; page < m_numPages; ++page)
208 {
209 // Calculate the sum of squared absolute values of each matrix page
210 res[page] = 0;
211 auto pagePtr = this->GetPagePtr(page);
212 for (size_t i = 0; i < m_numRows * m_numCols; i++)
214 auto absVal = std::abs(pagePtr[i]);
215 res[page] += absVal * absVal;
216 }
217 res[page] = sqrt(res[page]);
218 }
219 return res;
220}
221
222template <class T>
223MatrixArray<T>
225 const MatrixArray<T>& rMatrix) const
226{
227 NS_ASSERT_MSG(lMatrix.m_numPages == 1 && rMatrix.m_numPages == 1,
228 "The left and right MatrixArray should have only one page.");
229 NS_ASSERT_MSG(lMatrix.m_numCols == m_numRows,
230 "Left vector numCols and this MatrixArray numRows mismatch.");
231 NS_ASSERT_MSG(m_numCols == rMatrix.m_numRows,
232 "Right vector numRows and this MatrixArray numCols mismatch.");
233
234 MatrixArray<T> res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages};
235
236#ifdef HAVE_EIGEN3
238 ConstEigenMatrix<T> lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols);
239 ConstEigenMatrix<T> rMatrixEigen(rMatrix.GetPagePtr(0), rMatrix.m_numRows, rMatrix.m_numCols);
240#endif
241
242 for (size_t page = 0; page < m_numPages; ++page)
244#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
245
246 ConstEigenMatrix<T> matrixEigen(GetPagePtr(page), m_numRows, m_numCols);
247 EigenMatrix<T> resEigenMap(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
248
249 resEigenMap = lMatrixEigen * matrixEigen * rMatrixEigen;
250
251#else // Eigen not found or Eigen optimizations not enabled
252
253 size_t matrixOffset = page * m_numRows * m_numCols;
254 for (size_t resRow = 0; resRow < res.m_numRows; ++resRow)
255 {
256 for (size_t resCol = 0; resCol < res.m_numCols; ++resCol)
257 {
258 // create intermediate row result, a multiply of resRow row of lMatrix and each
259 // column of this matrix
260 std::valarray<T> interRes(m_numCols);
261 for (size_t thisCol = 0; thisCol < m_numCols; ++thisCol)
262 {
263 interRes[thisCol] =
264 (lMatrix
265 .m_values[std::slice(resRow, lMatrix.m_numCols, lMatrix.m_numRows)] *
266 m_values[std::slice(matrixOffset + thisCol * m_numRows, m_numRows, 1)])
267 .sum();
268 }
269 // multiply intermediate results and resCol column of the rMatrix
270 res(resRow, resCol, page) =
271 (interRes *
272 rMatrix.m_values[std::slice(resCol * rMatrix.m_numRows, rMatrix.m_numRows, 1)])
273 .sum();
274 }
275 }
276#endif
277 }
278 return res;
279}
280
281template <class T>
282template <bool EnableBool, typename>
283MatrixArray<T>
285{
286 MatrixArray<std::complex<double>> retMatrix = this->Transpose();
287
288 for (size_t index = 0; index < this->GetSize(); ++index)
289 {
290 retMatrix.m_values[index] = std::conj(retMatrix.m_values[index]);
291 }
292 return retMatrix;
293}
294
295template <class T>
297MatrixArray<T>::MakeNCopies(size_t nCopies) const
298{
299 NS_ASSERT_MSG(m_numPages == 1, "The MatrixArray should have only one page to be copied.");
300 auto copiedMatrix = MatrixArray<T>{m_numRows, m_numCols, nCopies};
301 for (size_t copy = 0; copy < nCopies; copy++)
302 {
303 for (size_t i = 0; i < m_numRows * m_numCols; i++)
304 {
305 copiedMatrix.GetPagePtr(copy)[i] = m_values[i];
306 }
307 }
308 return copiedMatrix;
309}
310
311template <class T>
314{
315 NS_ASSERT_MSG(page < m_numPages, "The page to extract from the MatrixArray is out of bounds.");
316 auto extractedPage = MatrixArray<T>{m_numRows, m_numCols, 1};
317
318 for (size_t i = 0; i < m_numRows * m_numCols; ++i)
319 {
320 extractedPage.m_values[i] = GetPagePtr(page)[i];
321 }
322 return extractedPage;
323}
324
325template <class T>
328{
329 auto jointMatrix =
330 MatrixArray<T>{pages.front().GetNumRows(), pages.front().GetNumCols(), pages.size()};
331 for (size_t page = 0; page < jointMatrix.GetNumPages(); page++)
332 {
333 NS_ASSERT_MSG(pages[page].GetNumRows() == jointMatrix.GetNumRows(),
334 "All page matrices should have the same number of rows");
335 NS_ASSERT_MSG(pages[page].GetNumCols() == jointMatrix.GetNumCols(),
336 "All page matrices should have the same number of columns");
337 NS_ASSERT_MSG(pages[page].GetNumPages() == 1,
338 "All page matrices should have a single page");
339
340 size_t i = 0;
341 for (auto a : pages[page].GetValues())
342 {
343 jointMatrix.GetPagePtr(page)[i] = a;
344 i++;
345 }
346 }
347 return jointMatrix;
348}
349
350template <class T>
352MatrixArray<T>::IdentityMatrix(const size_t size, const size_t pages)
353{
354 auto identityMatrix = MatrixArray<T>{size, size, pages};
355 for (std::size_t page = 0; page < pages; page++)
356 {
357 for (std::size_t i = 0; i < size; i++)
358 {
359 identityMatrix(i, i, page) = 1.0;
360 }
361 }
362 return identityMatrix;
363}
364
365template <class T>
368{
369 NS_ASSERT_MSG(likeme.GetNumRows() == likeme.GetNumCols(), "Template array is not square.");
370 return IdentityMatrix(likeme.GetNumRows(), likeme.GetNumPages());
371}
372
374 const;
376template class MatrixArray<double>;
377template class MatrixArray<int>;
378
379} // 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.