revision: v0.14
SmallMatrix< numberOfRows, numberOfColumns > Class Template Reference

Data type for small dense matrix. More...

#include <SmallMatrix.h>

Public Member Functions

 SmallMatrix ()
 Constructs a matrix of size n-rows by m-columns. More...
 
 SmallMatrix (const SmallVector< numberOfRows > &other)
 
 SmallMatrix (const Mdouble &c)
 Constructs a matrix of size n-rows by m-columns and initialises all entry to a constant. More...
 
 SmallMatrix (const std::initializer_list< SmallVector< numberOfRows >> &entries)
 
 SmallMatrix (const SmallMatrix &other)
 Construct and copy Matrix from another Matrix i.e. B(A) where B and A are both matrices. More...
 
 SmallMatrix (std::array< SmallVector< numberOfRows >, numberOfColumns > entries)
 Glues one or more vectors with the same number of rows together. More...
 
 SmallMatrix (SmallMatrix &&other)
 Move Matrix from another Matrix. More...
 
Mdoubleoperator() (unsigned int n, unsigned int m)
 defines the operator(n,m) to access the element on row n and column m More...
 
const Mdoubleoperator() (unsigned int n, unsigned int m) const
 defines the operator(n,m) to access the element on row n and column m More...
 
Mdoubleoperator[] (const unsigned int n)
 Access the n linear element in the matrix. More...
 
const Mdoubleoperator[] (const unsigned int n) const
 
SmallVector< numberOfRows > operator* (SmallVector< numberOfColumns > &right)
 Defines Matrix A times vector B and return vector C i.e. C_,j= A_ij B_,j. More...
 
SmallVector< numberOfRows > operator* (SmallVector< numberOfColumns > &right) const
 
SmallMatrix operator* (const Mdouble &right) const
 Does matrix A_ij=scalar*B_ij. More...
 
template<unsigned int K>
SmallMatrix< numberOfRows, K > operator* (const SmallMatrix< numberOfColumns, K > &other)
 Does matrix A_ij = B_ik * C_kj. More...
 
template<unsigned int K>
SmallMatrix< numberOfRows, K > operator* (const SmallMatrix< numberOfColumns, K > &other) const
 
SmallMatrixoperator+= (const SmallMatrix &other)
 
SmallMatrixoperator-= (const SmallMatrix &other)
 
SmallMatrix operator+ (const SmallMatrix &other) const
 
SmallMatrix operator- (const SmallMatrix &other) const
 
SmallMatrix operator- () const
 
SmallMatrixoperator*= (const Mdouble &scalar)
 Does matrix A_ij=scalar*A_ij. More...
 
SmallMatrixoperator*= (const SmallMatrix< numberOfColumns, numberOfColumns > &other)
 Does matrix A_ij = A_ik * B_kj note that other must be square because this is a fixed-size matrix. More...
 
SmallMatrixoperator/= (const Mdouble &scalar)
 Does matrix A_ij=scalar*A_ij. More...
 
SmallMatrix operator/ (const Mdouble &scalar) const
 this does element by divided by a scalar More...
 
SmallMatrixoperator= (const SmallMatrix &right)
 Assigns one matrix to another. More...
 
SmallMatrixoperator= (SmallMatrix &&right)
 Assigns one matrix to another. More...
 
SmallVector< numberOfRows > computeWedgeStuffVector () const
 computeWedgeStuffVector. More...
 
void axpy (Mdouble a, const SmallMatrix &x)
 Applies the matrix y=ax + y, where x is another matrix and a is a scalar. More...
 
unsigned int size () const
 Get total number of Matrix entries. More...
 
unsigned int getNumberOfRows () const
 Get the number of rows. More...
 
unsigned int getNRows () const
 
unsigned int getNumberOfColumns () const
 Get the number of columns. More...
 
unsigned int getNCols () const
 
SmallVector< numberOfRows > getColumn (unsigned int j) const
 get the j^th column More...
 
SmallVector< numberOfColumns > getRow (unsigned int i) const
 get the i^th row More...
 
SmallMatrix LUfactorisation () const
 Return the LUfactorisation of the matrix. More...
 
Mdouble determinant () const
 
SmallMatrix inverse () const
 return the inverse in the vector result. The size of result matches the matrix. More...
 
SmallMatrix< numberOfColumns, numberOfRows > transpose () const
 
template<unsigned int numberOfRightHandSideColumns>
void solve (SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &B) const
 solves Ax=B where A is the current matrix and B is passed in. The result is returned in B. More...
 
void solve (SmallVector< numberOfRows > &b) const
 solves Ax=b where A is the current matrix and NumericalVector b is the input parameter. The result is returned in b. More...
 
Mdoubledata ()
 
const Mdoubledata () const
 

Private Attributes

std::array< Mdouble, numberOfRows *numberOfColumns > data_
 The actually data of the matrix class. More...
 

Detailed Description

template<unsigned int numberOfRows, unsigned int numberOfColumns>
class SmallMatrix< numberOfRows, numberOfColumns >

Data type for small dense matrix.

Stores small dense matrix efficiently. It only store Mdoubles as this is the main type linear algebra is done on in hpGEM It stores the matrix in fortran style (column-major) to give quicker access to extern BLAS libraries. For example, the order they are stored in a 2x2 matrix is 0 2 1 3

Constructor & Destructor Documentation

◆ SmallMatrix() [1/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( )
inline

Constructs a matrix of size n-rows by m-columns.

73  : data_()
74  {
75  }

◆ SmallMatrix() [2/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( const SmallVector< numberOfRows > &  other)
inline
78  : data_()
79  {
80  logger.assert_debug(numberOfColumns == 1, "Trying to construct a matrix with more than 1 columns from a vector");
81  std::copy(other.data(), other.data() + numberOfRows, data_.begin());
82  }

References SmallVector< numberOfRows >::data(), SmallMatrix< numberOfRows, numberOfColumns >::data_, and logger.

◆ SmallMatrix() [3/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( const Mdouble c)
inline

Constructs a matrix of size n-rows by m-columns and initialises all entry to a constant.

86  : data_()
87  {
88  data_.fill(c);
89  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ SmallMatrix() [4/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( const std::initializer_list< SmallVector< numberOfRows >> &  entries)
inline
92  : data_()
93  {
94  logger.assert_debug(entries.size() == numberOfColumns, "expected a matrix with % "
95  "columns, but got a matrix with % columns", numberOfColumns,
96  entries.size());
97  unsigned int column = 0;
98  for (const SmallVector<numberOfRows>& entry : entries)
99  {
100  for (unsigned int i = 0; i < numberOfRows; ++i)
101  {
102  (*this)(i, column) = entry[i];
103  }
104  ++column;
105  }
106  }

References constants::i, and logger.

◆ SmallMatrix() [5/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( const SmallMatrix< numberOfRows, numberOfColumns > &  other)
inline

Construct and copy Matrix from another Matrix i.e. B(A) where B and A are both matrices.

110  : data_()
111  {
112  std::copy(other.data_.begin(), other.data_.end(), data_.begin());
113  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ SmallMatrix() [6/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( std::array< SmallVector< numberOfRows >, numberOfColumns >  entries)
inline

Glues one or more vectors with the same number of rows together.

117  : data_()
118  {
119  for (unsigned int i = 0; i < numberOfRows; ++i)
120  {
121  for (unsigned int j = 0; j < numberOfColumns; ++j)
122  {
123  (*this)(i, j) = entries[j][i];
124  }
125  }
126  }

References constants::i.

◆ SmallMatrix() [7/7]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns >::SmallMatrix ( SmallMatrix< numberOfRows, numberOfColumns > &&  other)
inline

Move Matrix from another Matrix.

130  : data_(std::move(other.data_))
131  {
132  }

Member Function Documentation

◆ axpy()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
void SmallMatrix< numberOfRows, numberOfColumns >::axpy ( Mdouble  a,
const SmallMatrix< numberOfRows, numberOfColumns > &  x 
)
inline

Applies the matrix y=ax + y, where x is another matrix and a is a scalar.

267  {
268  for (unsigned int i = 0; i < numberOfRows * numberOfColumns; ++i)
269  {
270  data_[i] += a * x[i];
271  }
272  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_, and constants::i.

Referenced by main().

◆ computeWedgeStuffVector()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallVector< numberOfRows > SmallMatrix< numberOfRows, numberOfColumns >::computeWedgeStuffVector

computeWedgeStuffVector.

202 {
203  //copied from MiddleSizeMatrix to prevent constructing a temporary MiddleSizeMatrix
204  logger.assert_debug(numberOfColumns == numberOfRows - 1,
205  "Matrix has wrong dimensions to construct the wedge stuff vector");
207 
208  switch (numberOfRows)
209  {
210  case 2:
211  result[0] = -(*this)(1, 0);
212  result[1] = +(*this)(0, 0);
213  break;
214  case 3:
215  result[0] = (*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1);
216  result[1] = (*this)(0, 1) * (*this)(2, 0) - (*this)(0, 0) * (*this)(2, 1); // includes minus sign already!
217  result[2] = (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(0, 1);
218  break;
219  case 4:
220  result[0] = (*this)(1, 0) * (-(*this)(2, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(2, 2)) +
221  (*this)(2, 0) * ((*this)(1, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(1, 2)) +
222  (*this)(3, 0) * (-(*this)(1, 1) * (*this)(2, 2) + (*this)(2, 1) * (*this)(1, 2));
223 
224  result[1] = (*this)(0, 0) * ((*this)(2, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(2, 2)) +
225  (*this)(2, 0) * (-(*this)(0, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(0, 2)) +
226  (*this)(3, 0) * ((*this)(0, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(0, 2));
227  result[2] = (*this)(0, 0) * (-(*this)(1, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(1, 2)) +
228  (*this)(1, 0) * ((*this)(0, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(0, 2)) +
229  (*this)(3, 0) * (-(*this)(0, 1) * (*this)(1, 2) + (*this)(1, 1) * (*this)(0, 2));
230  result[3] = (*this)(0, 0) * ((*this)(1, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(1, 2)) +
231  (*this)(1, 0) * (-(*this)(0, 1) * (*this)(2, 2) + (*this)(2, 1) * (*this)(0, 2)) +
232  (*this)(2, 0) * ((*this)(0, 1) * (*this)(1, 2) - (*this)(1, 1) * (*this)(0, 2));
233  break;
234  default:
235  logger(ERROR, "Wedge product not implemented for this dimension");
236  } //end switch
237 
238  return (result);
239 }

References ERROR, and logger.

◆ data() [1/2]

◆ data() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
const Mdouble* SmallMatrix< numberOfRows, numberOfColumns >::data ( ) const
inline
357  {
358  return data_.data();
359  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ determinant()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
Mdouble SmallMatrix< numberOfRows, numberOfColumns >::determinant
262 {
263  logger.assert_debug(numberOfRows == numberOfColumns, "Matrix should be square to have a determinant!");
264 
265  switch (numberOfRows)
266  {
267  case 0:
268  return 1;
269  case 1:
270  return (*this)(0, 0);
271  case 2:
272  return (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
273 
274  case 3:
275  return (*this)(0, 0) * ((*this)(1, 1) * (*this)(2, 2) - (*this)(1, 2) * (*this)(2, 1)) -
276  (*this)(0, 1) * ((*this)(1, 0) * (*this)(2, 2) - (*this)(2, 0) * (*this)(1, 2)) +
277  (*this)(0, 2) * ((*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1));
278 
279  case 4:
280  return ((*this)(3, 0) * (*this)(2, 1) * (*this)(0, 3) - (*this)(2, 0) * (*this)(3, 1) * (*this)(0, 3)) *
281  (*this)(1, 2) +
282  (-(*this)(3, 0) * (*this)(0, 3) * (*this)(2, 2) + (*this)(2, 0) * (*this)(0, 3) * (*this)(3, 2)) *
283  (*this)(1, 1) +
284  ((*this)(3, 1) * (*this)(0, 3) * (*this)(2, 2) - (*this)(2, 1) * (*this)(0, 3) * (*this)(3, 2)) *
285  (*this)(1, 0) +
286  (-(*this)(3, 0) * (*this)(2, 1) * (*this)(1, 3) + (*this)(2, 0) * (*this)(3, 1) * (*this)(1, 3) +
287  (-(*this)(2, 0) * (*this)(3, 3) + (*this)(3, 0) * (*this)(2, 3)) * (*this)(1, 1) +
288  ((*this)(2, 1) * (*this)(3, 3) - (*this)(3, 1) * (*this)(2, 3)) * (*this)(1, 0)) * (*this)(0, 2) +
289  ((*this)(3, 0) * (*this)(1, 3) * (*this)(2, 2) - (*this)(2, 0) * (*this)(1, 3) * (*this)(3, 2) +
290  ((*this)(2, 0) * (*this)(3, 3) - (*this)(3, 0) * (*this)(2, 3)) * (*this)(1, 2) +
291  (-(*this)(2, 2) * (*this)(3, 3) + (*this)(2, 3) * (*this)(3, 2)) * (*this)(1, 0)) * (*this)(0, 1) +
292  (-(*this)(3, 1) * (*this)(1, 3) * (*this)(2, 2) + (*this)(2, 1) * (*this)(1, 3) * (*this)(3, 2) +
293  ((*this)(3, 1) * (*this)(2, 3) - (*this)(2, 1) * (*this)(3, 3)) * (*this)(1, 2) +
294  (*this)(1, 1) * ((*this)(2, 2) * (*this)(3, 3) - (*this)(2, 3) * (*this)(3, 2))) * (*this)(0, 0);
295  // ... says Maple; this can possibly be done more efficiently,
296  // maybe even with LU (with pivoting, though...)
297  default:
298  logger(ERROR, "Computing the Determinant for size % is not implemented", numberOfRows);
299  break;
300  }
301  return 0;
302 }

References ERROR, and logger.

◆ getColumn()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallVector<numberOfRows> SmallMatrix< numberOfRows, numberOfColumns >::getColumn ( unsigned int  j) const
inline

get the j^th column

305  {
306  logger.assert_debug(j < numberOfColumns, "Asked for column %, but there are only % columns", j, numberOfColumns);
307  return SmallVector<numberOfRows>(data() + j * numberOfRows);
308  }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), and logger.

◆ getNCols()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
unsigned int SmallMatrix< numberOfRows, numberOfColumns >::getNCols ( ) const
inline
299  {
300  return getNumberOfColumns();
301  }

References SmallMatrix< numberOfRows, numberOfColumns >::getNumberOfColumns().

Referenced by main().

◆ getNRows()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
unsigned int SmallMatrix< numberOfRows, numberOfColumns >::getNRows ( ) const
inline
Deprecated:
Does not conform naming convention, please use getNumberOfRows instead.
288  {
289  return getNumberOfRows();
290  }

References SmallMatrix< numberOfRows, numberOfColumns >::getNumberOfRows().

◆ getNumberOfColumns()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
unsigned int SmallMatrix< numberOfRows, numberOfColumns >::getNumberOfColumns ( ) const
inline

Get the number of columns.

294  {
295  return numberOfColumns;
296  }

Referenced by SmallMatrix< numberOfRows, numberOfColumns >::getNCols().

◆ getNumberOfRows()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
unsigned int SmallMatrix< numberOfRows, numberOfColumns >::getNumberOfRows ( ) const
inline

Get the number of rows.

282  {
283  return numberOfRows;
284  }

Referenced by SmallMatrix< numberOfRows, numberOfColumns >::getNRows(), and main().

◆ getRow()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallVector<numberOfColumns> SmallMatrix< numberOfRows, numberOfColumns >::getRow ( unsigned int  i) const
inline

get the i^th row

312  {
313  logger.assert_debug(i < numberOfRows, "Asked for row %, but there are only % rows", i, numberOfRows);
315  for (unsigned int j = 0; j < numberOfColumns; ++j)
316  {
317  result[j] = (*this)(i, j);
318  }
319  return result;
320  }

References constants::i, and logger.

◆ inverse()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns > SmallMatrix< numberOfRows, numberOfColumns >::inverse

return the inverse in the vector result. The size of result matches the matrix.

306 {
307  logger.assert_debug(numberOfRows == numberOfColumns, "Cannot invert a non-square matrix");
309 
310  int nr = numberOfRows;
311  int nc = numberOfColumns;
312 
313  int nPivot = numberOfRows;
314  int iPivot[nPivot];
315 
316  int info = 0;
317 
318  dgetrf_(&nr, &nc, result.data(), &nr, iPivot, &info);
319 
320  int lwork = numberOfRows * numberOfColumns;
321  SmallMatrix work;
322  dgetri_(&nc, result.data(), &nc, iPivot, work.data(), &lwork, &info);
323 
324  return result;
325 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), dgetrf_(), dgetri_(), and logger.

◆ LUfactorisation()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns > SmallMatrix< numberOfRows, numberOfColumns >::LUfactorisation

Return the LUfactorisation of the matrix.

243 {
244  int nr = numberOfRows;
245  int nc = numberOfColumns;
246  int nPivot = std::min(numberOfRows, numberOfColumns);
247  int iPivot[nPivot];
248 
249  SmallMatrix result(*this);
250 
251  int info;
252 
253  dgetrf_(&nr, &nc, result.data(), &nr, iPivot, &info);
254 
255  return result;
256 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), and dgetrf_().

◆ operator()() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
Mdouble& SmallMatrix< numberOfRows, numberOfColumns >::operator() ( unsigned int  n,
unsigned int  m 
)
inline

defines the operator(n,m) to access the element on row n and column m

136  {
137  logger.assert_debug(n < numberOfRows, "Requested row number % for a matrix with only % rows", n, numberOfRows);
138  logger.assert_debug(m < numberOfColumns, "Requested column number % for a matrix with only % columns", m,
139  numberOfColumns);
140  return data_[n + m * numberOfRows];
141  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_, logger, and n.

◆ operator()() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
const Mdouble& SmallMatrix< numberOfRows, numberOfColumns >::operator() ( unsigned int  n,
unsigned int  m 
) const
inline

defines the operator(n,m) to access the element on row n and column m

145  {
146  logger.assert_debug(n < numberOfRows, "Requested row number % for a matrix with only % rows", n, numberOfRows);
147  logger.assert_debug(m < numberOfColumns, "Requested column number % for a matrix with only % columns", m,
148  numberOfColumns);
149  return data_[n + m * numberOfRows];
150  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_, logger, and n.

◆ operator*() [1/5]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix SmallMatrix< numberOfRows, numberOfColumns >::operator* ( const Mdouble right) const
inline

Does matrix A_ij=scalar*B_ij.

174  {
175  SmallMatrix result;
176  std::transform(data_.begin(), data_.end(), result.data_.begin(),
177  std::bind(std::multiplies<Mdouble>(), std::placeholders::_1, right));
178  return result;
179  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator*() [2/5]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
template<unsigned int K>
SmallMatrix< numberOfRows, K > SmallMatrix< numberOfRows, numberOfColumns >::operator* ( const SmallMatrix< numberOfColumns, K > &  other)

Does matrix A_ij = B_ik * C_kj.

142 {
143  int i = numberOfRows;
144  int j = numberOfColumns;
145  int k = K;
146 
147  if (numberOfColumns == 0)
148  {
149  logger(WARN, "Trying to multiply a matrix with a matrix without any columns.");
151  }
152  //The result of the matrix is left.numberOfRows, right.numberOfColumns()
154 
155  double d_one = 1.0;
156  double d_zero = 0.0;
157 
158  //Let the actual multiplication be done by Fortran
159  dgemm_("N", "N", &i, &k, &j, &d_one, this->data(), &i, const_cast<double*>(other.data()), &j, &d_zero, C.data(),
160  &i);
161 
162  return C;
163 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), dgemm_(), constants::i, logger, and WARN.

◆ operator*() [3/5]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
template<unsigned int K>
SmallMatrix< numberOfRows, K > SmallMatrix< numberOfRows, numberOfColumns >::operator* ( const SmallMatrix< numberOfColumns, K > &  other) const
169 {
170  int i = numberOfRows;
171  int j = numberOfColumns;
172  int k = K;
173 
174  if (numberOfColumns == 0)
175  {
176  logger(WARN, "Trying to multiply a matrix with a matrix without any columns.");
178  }
179  //The result of the matrix is left.Nrows, right.NCols()
181 
182  double d_one = 1.0;
183  double d_zero = 0.0;
184 
185  //Let the actual multiplication be done by Fortran
186  dgemm_("N", "N", &i, &k, &j, &d_one, const_cast<double*>(this->data()), &i, const_cast<double*>(other.data()), &j,
187  &d_zero, C.data(), &i);
188 
189  return C;
190 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), dgemm_(), constants::i, logger, and WARN.

◆ operator*() [4/5]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallVector< numberOfRows > SmallMatrix< numberOfRows, numberOfColumns >::operator* ( SmallVector< numberOfColumns > &  right)

Defines Matrix A times vector B and return vector C i.e. C_,j= A_ij B_,j.

82 {
83  if (numberOfRows == 0)
84  {
85  logger(WARN, "Trying to multiply a vector with a matrix without any rows.");
87  }
88  if (numberOfColumns == 0)
89  {
90  logger(WARN, "Trying to multiply a vector with a matrix without any columns.");
92  }
93  int nr = numberOfRows;
94  int nc = numberOfColumns;
95 
96  int i_one = 1;
97  double d_one = 1.0;
98  double d_zero = 0.0;
99 
101 
102  logger(DEBUG, "Matrix size: % x % \n Vector size: %", nr, nc, right.size());
103 
104  dgemv_("N", &nr, &nc, &d_one, this->data(), &nr, right.data(), &i_one, &d_zero, result.data(), &i_one);
105  return result;
106 }

References SmallVector< numberOfRows >::data(), DEBUG, dgemv_(), logger, SmallVector< numberOfRows >::size(), and WARN.

◆ operator*() [5/5]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallVector< numberOfRows > SmallMatrix< numberOfRows, numberOfColumns >::operator* ( SmallVector< numberOfColumns > &  right) const
111 {
112  if (numberOfRows == 0)
113  {
114  logger(WARN, "Trying to multiply a vector with a matrix without any rows.");
115  return SmallVector<numberOfRows>();
116  }
117  if (numberOfColumns == 0)
118  {
119  logger(WARN, "Trying to multiply a vector with a matrix without any columns.");
120  return SmallVector<numberOfRows>();
121  }
122  int nr = numberOfRows;
123  int nc = numberOfColumns;
124 
125  int i_one = 1;
126  double d_one = 1.0;
127  double d_zero = 0.0;
128 
130 
131  logger(DEBUG, "Matrix size: % x % \n Vector size: %", nr, nc, right.size());
132 
133  dgemv_("N", &nr, &nc, &d_one, (const_cast<double*>(this->data())), &nr, right.data(), &i_one, &d_zero,
134  result.data(), &i_one);
135  return result;
136 }

References SmallVector< numberOfRows >::data(), DEBUG, dgemv_(), logger, SmallVector< numberOfRows >::size(), and WARN.

◆ operator*=() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator*= ( const Mdouble scalar)
inline

Does matrix A_ij=scalar*A_ij.

221  {
222  std::transform(data_.begin(), data_.end(), data_.begin(),
223  std::bind(std::multiplies<Mdouble>(), std::placeholders::_1, scalar));
224  return *this;
225  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator*=() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix< numberOfRows, numberOfColumns > & SmallMatrix< numberOfRows, numberOfColumns >::operator*= ( const SmallMatrix< numberOfColumns, numberOfColumns > &  other)

Does matrix A_ij = A_ik * B_kj note that other must be square because this is a fixed-size matrix.

195 {
196  //blas does not support in-place multiply
197  return (*this) = (*this) * other;
198 }

◆ operator+()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix SmallMatrix< numberOfRows, numberOfColumns >::operator+ ( const SmallMatrix< numberOfRows, numberOfColumns > &  other) const
inline
201  {
202  SmallMatrix result;
203  std::transform(data_.begin(), data_.end(), other.data_.begin(), result.data_.begin(), std::plus<Mdouble>());
204  return result;
205  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator+=()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator+= ( const SmallMatrix< numberOfRows, numberOfColumns > &  other)
inline
189  {
190  std::transform(data_.begin(), data_.end(), other.data_.begin(), data_.begin(), std::plus<Mdouble>());
191  return *this;
192  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator-() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix SmallMatrix< numberOfRows, numberOfColumns >::operator- ( ) const
inline
215  {
216  return *this * -1.;
217  }

◆ operator-() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix SmallMatrix< numberOfRows, numberOfColumns >::operator- ( const SmallMatrix< numberOfRows, numberOfColumns > &  other) const
inline
208  {
209  SmallMatrix result;
210  std::transform(data_.begin(), data_.end(), other.data_.begin(), result.data_.begin(), std::minus<Mdouble>());
211  return result;
212  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator-=()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator-= ( const SmallMatrix< numberOfRows, numberOfColumns > &  other)
inline
195  {
196  std::transform(data_.begin(), data_.end(), other.data_.begin(), data_.begin(), std::minus<Mdouble>());
197  return *this;
198  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator/()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix SmallMatrix< numberOfRows, numberOfColumns >::operator/ ( const Mdouble scalar) const
inline

this does element by divided by a scalar

241  {
242  SmallMatrix result;
243  std::transform(data_.begin(), data_.end(), result.data_.begin(),
244  std::bind(std::divides<Mdouble>(), std::placeholders::_1, scalar));
245  return result;
246  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator/=()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator/= ( const Mdouble scalar)
inline

Does matrix A_ij=scalar*A_ij.

233  {
234  std::transform(data_.begin(), data_.end(), data_.begin(),
235  std::bind(std::divides<Mdouble>(), std::placeholders::_1, scalar));
236  return *this;
237  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator=() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator= ( const SmallMatrix< numberOfRows, numberOfColumns > &  right)
inline

Assigns one matrix to another.

250  {
251  std::copy(right.data_.begin(), right.data_.end(), data_.begin());
252  return *this;
253  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator=() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix& SmallMatrix< numberOfRows, numberOfColumns >::operator= ( SmallMatrix< numberOfRows, numberOfColumns > &&  right)
inline

Assigns one matrix to another.

257  {
258  std::move(right.data_.begin(), right.data_.end(), data_.begin());
259  return *this;
260  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_.

◆ operator[]() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
Mdouble& SmallMatrix< numberOfRows, numberOfColumns >::operator[] ( const unsigned int  n)
inline

Access the n linear element in the matrix.

154  {
155  logger.assert_debug(n < numberOfRows * numberOfColumns, "Requested entry % for a matrix with only % entries", n,
156  numberOfRows * numberOfColumns);
157  return data_[n];
158  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_, logger, and n.

◆ operator[]() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
const Mdouble& SmallMatrix< numberOfRows, numberOfColumns >::operator[] ( const unsigned int  n) const
inline
161  {
162  logger.assert_debug(n < numberOfRows * numberOfColumns, "Requested entry % for a matrix with only % entries", n,
163  numberOfRows * numberOfColumns);
164  return data_[n];
165  }

References SmallMatrix< numberOfRows, numberOfColumns >::data_, logger, and n.

◆ size()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
unsigned int SmallMatrix< numberOfRows, numberOfColumns >::size ( ) const
inline

Get total number of Matrix entries.

276  {
277  return numberOfRows * numberOfColumns;
278  }

Referenced by main().

◆ solve() [1/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
template<unsigned int numberOfRightHandSideColumns>
void SmallMatrix< numberOfRows, numberOfColumns >::solve ( SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &  B) const

solves Ax=B where A is the current matrix and B is passed in. The result is returned in B.

330 {
331  logger.assert_debug(numberOfRows == numberOfColumns, "can only solve for square matrixes");
332 
333  int n = numberOfRows;
334  int nrhs = numberOfRightHandSideColumns;
335  int info;
336 
337  int IPIV[numberOfRows];
339  dgesv_(&n, &nrhs, matThis.data(), &n, IPIV, B.data(), &n, &info);
340 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), dgesv_(), logger, and n.

Referenced by SuperQuadricParticle::computeContactPoint().

◆ solve() [2/2]

template<unsigned int numberOfRows, unsigned int numberOfColumns>
void SmallMatrix< numberOfRows, numberOfColumns >::solve ( SmallVector< numberOfRows > &  b) const

solves Ax=b where A is the current matrix and NumericalVector b is the input parameter. The result is returned in b.

344 {
345  logger.assert_debug(numberOfRows == numberOfColumns, "can only solve for square matrixes");
346 
347  int n = numberOfRows;
348  int nrhs = 1;
349  int info;
350 
351  int IPIV[numberOfRows];
352  SmallMatrix matThis = *this;
353  dgesv_(&n, &nrhs, matThis.data(), &n, IPIV, b.data(), &n, &info);
354 }

References SmallMatrix< numberOfRows, numberOfColumns >::data(), SmallVector< numberOfRows >::data(), dgesv_(), logger, and n.

◆ transpose()

template<unsigned int numberOfRows, unsigned int numberOfColumns>
SmallMatrix<numberOfColumns, numberOfRows> SmallMatrix< numberOfRows, numberOfColumns >::transpose ( ) const
inline
331  {
333  for (unsigned int i = 0; i < numberOfRows; ++i)
334  {
335  for (unsigned int j = 0; j < numberOfColumns; ++j)
336  {
337  result(j, i) = (*this)(i, j);
338  }
339  }
340  return result;
341  }

References constants::i.

Referenced by SuperQuadricParticle::getCurvature(), ShapeGradientHessianTester::testCushion(), ShapeGradientHessianTester::testEllipsoid(), and ShapeGradientHessianTester::testRoundedBeam().

Member Data Documentation

◆ data_


The documentation for this class was generated from the following files:
dgesv_
void dgesv_(int *N, int *NRHS, double *A, int *lda, int *IPIV, double *B, int *LDB, int *INFO)
This is used for solve Ax=B for x. Again this is from LAPACK.
SmallMatrix::data
Mdouble * data()
Definition: SmallMatrix.h:351
dgetrf_
void dgetrf_(int *M, int *N, double *A, int *lda, int *IPIV, int *INFO)
This is LU factorisation of the matrix A. This has been taken from LAPACK.
logger
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
dgemm_
int dgemm_(const char *transA, const char *transB, int *M, int *N, int *k, double *alpha, double *A, int *LDA, double *B, int *LDB, double *beta, double *C, int *LDC)
This is the gernal matrix multiplication from blas level 3.
SmallMatrix::getNumberOfColumns
unsigned int getNumberOfColumns() const
Get the number of columns.
Definition: SmallMatrix.h:293
SmallMatrix::data_
std::array< Mdouble, numberOfRows *numberOfColumns > data_
The actually data of the matrix class.
Definition: SmallMatrix.h:363
SmallMatrix::getNumberOfRows
unsigned int getNumberOfRows() const
Get the number of rows.
Definition: SmallMatrix.h:281
dgemv_
void dgemv_(const char *trans, int *m, int *n, double *alpha, double *A, int *LDA, double *x, int *incx, double *beta, double *y, int *incy)
This does matrix times vector and is from blas level 2.
SmallMatrix
Data type for small dense matrix.
Definition: SmallMatrix.h:68
ERROR
LL< Log::ERROR > ERROR
Error log level.
Definition: Logger.cc:53
WARN
LL< Log::WARN > WARN
Warning log level.
Definition: Logger.cc:54
dgetri_
void dgetri_(int *N, double *A, int *lda, int *IPIV, double *WORK, int *lwork, int *INFO)
This is the inverse calulation also from LAPACK. Calculates inverse if you pass it the LU factorisati...
Log::FATAL
@ FATAL
SmallVector::size
unsigned int size() const
Definition: SmallVector.h:233
SmallVector::data
const Mdouble * data() const
Definition: SmallVector.h:238
constants::i
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:32
SmallVector
Definition: SmallVector.h:62