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(size_t numRows, size_t numCols, size_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(size_t numRows, size_t numCols, const std::valarray<T>& values)
62 : ValArray<T>(numRows, numCols, values)
63{
64}
65
66template <class T>
67MatrixArray<T>::MatrixArray(size_t numRows, size_t numCols, std::valarray<T>&& values)
68 : ValArray<T>(numRows, numCols, std::move(values))
69{
70}
71
72template <class T>
74 size_t numCols,
75 size_t numPages,
76 const std::valarray<T>& values)
77 : ValArray<T>(numRows, numCols, numPages, values)
78{
79}
80
81template <class T>
83 size_t numCols,
84 size_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 (size_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 size_t rhsMatrixOffset = page * rhs.m_numRows * rhs.m_numCols;
112 for (size_t i = 0; i < res.m_numRows; ++i)
113 {
114 for (size_t j = 0; j < res.m_numCols; ++j)
115 {
116 res(i, j, page) = (m_values[std::slice(matrixOffset + i, m_numCols, m_numRows)] *
117 rhs.m_values[std::slice(rhsMatrixOffset + j * rhs.m_numRows,
118 rhs.m_numRows,
119 1)])
120 .sum();
121 }
122 }
123
124#endif
125 }
126 return res;
127}
128
129template <class T>
132{
133 // Create the matrix where m_numRows = this.m_numCols, m_numCols = this.m_numRows,
134 // m_numPages = this.m_numPages
135 MatrixArray<T> res{m_numCols, m_numRows, m_numPages};
136
137 for (size_t page = 0; page < m_numPages; ++page)
138 {
139#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
140
141 ConstEigenMatrix<T> thisMatrix(GetPagePtr(page), m_numRows, m_numCols);
142 EigenMatrix<T> resEigenMatrix(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
143 resEigenMatrix = thisMatrix.transpose();
144
145#else // Eigen not found or Eigen optimizations not enabled
146
147 size_t matrixIndex = page * m_numRows * m_numCols;
148 for (size_t i = 0; i < m_numRows; ++i)
149 {
150 res.m_values[std::slice(matrixIndex + i * res.m_numRows, res.m_numRows, 1)] =
151 m_values[std::slice(matrixIndex + i, m_numCols, m_numRows)];
152 }
153
154#endif
155 }
156 return res;
157}
158
159template <class T>
162 const MatrixArray<T>& rMatrix) const
163{
164 NS_ASSERT_MSG(lMatrix.m_numPages == 1 && rMatrix.m_numPages == 1,
165 "The left and right MatrixArray should have only one page.");
166 NS_ASSERT_MSG(lMatrix.m_numCols == m_numRows,
167 "Left vector numCols and this MatrixArray numRows mismatch.");
168 NS_ASSERT_MSG(m_numCols == rMatrix.m_numRows,
169 "Right vector numRows and this MatrixArray numCols mismatch.");
170
171 MatrixArray<T> res{lMatrix.m_numRows, rMatrix.m_numCols, m_numPages};
172
173#ifdef HAVE_EIGEN3
174
175 ConstEigenMatrix<T> lMatrixEigen(lMatrix.GetPagePtr(0), lMatrix.m_numRows, lMatrix.m_numCols);
176 ConstEigenMatrix<T> rMatrixEigen(rMatrix.GetPagePtr(0), rMatrix.m_numRows, rMatrix.m_numCols);
177#endif
178
179 for (size_t page = 0; page < m_numPages; ++page)
180 {
181#ifdef HAVE_EIGEN3 // Eigen found and Eigen optimizations enabled
182
183 ConstEigenMatrix<T> matrixEigen(GetPagePtr(page), m_numRows, m_numCols);
184 EigenMatrix<T> resEigenMap(res.GetPagePtr(page), res.m_numRows, res.m_numCols);
185
186 resEigenMap = lMatrixEigen * matrixEigen * rMatrixEigen;
187
188#else // Eigen not found or Eigen optimizations not enabled
189
190 size_t matrixOffset = page * m_numRows * m_numCols;
191 for (size_t resRow = 0; resRow < res.m_numRows; ++resRow)
192 {
193 for (size_t resCol = 0; resCol < res.m_numCols; ++resCol)
194 {
195 // create intermediate row result, a multiply of resRow row of lMatrix and each
196 // column of this matrix
197 std::valarray<T> interRes(m_numCols);
198 for (size_t thisCol = 0; thisCol < m_numCols; ++thisCol)
199 {
200 interRes[thisCol] =
201 (lMatrix
202 .m_values[std::slice(resRow, lMatrix.m_numCols, lMatrix.m_numRows)] *
203 m_values[std::slice(matrixOffset + thisCol * m_numRows, m_numRows, 1)])
204 .sum();
205 }
206 // multiply intermediate results and resCol column of the rMatrix
207 res(resRow, resCol, page) =
208 (interRes *
209 rMatrix.m_values[std::slice(resCol * rMatrix.m_numRows, rMatrix.m_numRows, 1)])
210 .sum();
211 }
212 }
213#endif
214 }
215 return res;
216}
217
218template <class T>
219template <bool EnableBool, typename>
222{
223 MatrixArray<std::complex<double>> retMatrix = this->Transpose();
224
225 for (size_t index = 0; index < this->GetSize(); ++index)
226 {
227 retMatrix.m_values[index] = std::conj(retMatrix.m_values[index]);
228 }
229 return retMatrix;
230}
231
233 const;
235template class MatrixArray<double>;
236template class MatrixArray<int>;
237
238} // 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:272
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
T * GetPagePtr(size_t pageIndex)
Get a data pointer to a specific 2D array for use in linear algebra libraries.
Definition: val-array.h:522
size_t m_numCols
The size of the second dimension, i.e., the number of columns of each 2D array.
Definition: val-array.h:367
std::valarray< T > m_values
The data values.
Definition: val-array.h:370
size_t m_numRows
The size of the first dimension, i.e., the number of rows of each 2D array.
Definition: val-array.h:365
size_t m_numPages
The size of the third dimension, i.e., the number of 2D arrays.
Definition: val-array.h:369
#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.