MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SmallMatrix.h
Go to the documentation of this file.
1 /*
2  This file forms part of hpGEM. This package has been developed over a number of years by various people at the University of Twente and a full list of contributors can be found at
3  http://hpgem.org/about-the-code/team
4 
5  This code is distributed using BSD 3-Clause License. A copy of which can found below.
6 
7 
8  Copyright (c) 2014, University of Twente
9  All rights reserved.
10 
11  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
12 
13  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
14 
15  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
16 
17  3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
18 
19  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
20  */
21 
22 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
23 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
24 //
25 //Redistribution and use in source and binary forms, with or without
26 //modification, are permitted provided that the following conditions are met:
27 // * Redistributions of source code must retain the above copyright
28 // notice, this list of conditions and the following disclaimer.
29 // * Redistributions in binary form must reproduce the above copyright
30 // notice, this list of conditions and the following disclaimer in the
31 // documentation and/or other materials provided with the distribution.
32 // * Neither the name MercuryDPM nor the
33 // names of its contributors may be used to endorse or promote products
34 // derived from this software without specific prior written permission.
35 //
36 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
37 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
38 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
39 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
40 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
41 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
42 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
43 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
44 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
45 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
46 
47 //Note: This code is copied and adapted from hpGEM (see license above), version 22th of January 2016. It has been
48 //integrated into MercuryDPM at 16th of March 2017.
49 
50 #ifndef MERCURY_SMALLMATRIX_H
51 #define MERCURY_SMALLMATRIX_H
52 
53 #include "GeneralDefine.h"
54 #include "Logger.h"
55 #include "SmallVector.h"
56 
66 template<unsigned int numberOfRows, unsigned int numberOfColumns>
68 {
69 public:
70 
73  : data_()
74  {
75  }
76 
78  : data_()
79  {
80  logger.assert(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  }
83 
85  SmallMatrix(const Mdouble& c)
86  : data_()
87  {
88  data_.fill(c);
89  }
90 
91  SmallMatrix(const std::initializer_list<SmallVector<numberOfRows>>& entries)
92  : data_()
93  {
94  logger.assert(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  }
107 
109  SmallMatrix(const SmallMatrix& other)
110  : data_()
111  {
112  std::copy(other.data_.begin(), other.data_.end(), data_.begin());
113  }
114 
116  SmallMatrix(std::array<SmallVector<numberOfRows>, numberOfColumns> entries)
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  }
127 
130  : data_(std::move(other.data_))
131  {
132  }
133 
135  Mdouble& operator()(unsigned int n, unsigned int m)
136  {
137  logger.assert(n < numberOfRows, "Requested row number % for a matrix with only % rows", n, numberOfRows);
138  logger.assert(m < numberOfColumns, "Requested column number % for a matrix with only % columns", m,
139  numberOfColumns);
140  return data_[n + m * numberOfRows];
141  }
142 
144  const Mdouble& operator()(unsigned int n, unsigned int m) const
145  {
146  logger.assert(n < numberOfRows, "Requested row number % for a matrix with only % rows", n, numberOfRows);
147  logger.assert(m < numberOfColumns, "Requested column number % for a matrix with only % columns", m,
148  numberOfColumns);
149  return data_[n + m * numberOfRows];
150  }
151 
153  Mdouble& operator[](const unsigned int n)
154  {
155  logger.assert(n < numberOfRows * numberOfColumns, "Requested entry % for a matrix with only % entries", n,
156  numberOfRows * numberOfColumns);
157  return data_[n];
158  }
159 
160  const Mdouble& operator[](const unsigned int n) const
161  {
162  logger.assert(n < numberOfRows * numberOfColumns, "Requested entry % for a matrix with only % entries", n,
163  numberOfRows * numberOfColumns);
164  return data_[n];
165  }
166 
169 
171 
173  SmallMatrix operator*(const Mdouble& right) const
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  }
180 
182  template<unsigned int K>
184 
185  template<unsigned int K>
187 
189  {
190  std::transform(data_.begin(), data_.end(), other.data_.begin(), data_.begin(), std::plus<Mdouble>());
191  return *this;
192  }
193 
195  {
196  std::transform(data_.begin(), data_.end(), other.data_.begin(), data_.begin(), std::minus<Mdouble>());
197  return *this;
198  }
199 
200  SmallMatrix operator+(const SmallMatrix& other) const
201  {
202  SmallMatrix result;
203  std::transform(data_.begin(), data_.end(), other.data_.begin(), result.data_.begin(), std::plus<Mdouble>());
204  return result;
205  }
206 
207  SmallMatrix operator-(const SmallMatrix& other) const
208  {
209  SmallMatrix result;
210  std::transform(data_.begin(), data_.end(), other.data_.begin(), result.data_.begin(), std::minus<Mdouble>());
211  return result;
212  }
213 
215  {
216  return *this * -1.;
217  }
218 
221  {
222  std::transform(data_.begin(), data_.end(), data_.begin(),
223  std::bind(std::multiplies<Mdouble>(), std::placeholders::_1, scalar));
224  return *this;
225  }
226 
230 
233  {
234  std::transform(data_.begin(), data_.end(), data_.begin(),
235  std::bind(std::divides<Mdouble>(), std::placeholders::_1, scalar));
236  return *this;
237  }
238 
240  SmallMatrix operator/(const Mdouble& scalar) const
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  }
247 
250  {
251  std::copy(right.data_.begin(), right.data_.end(), data_.begin());
252  return *this;
253  }
254 
257  {
258  std::move(right.data_.begin(), right.data_.end(), data_.begin());
259  return *this;
260  }
261 
264 
266  void axpy(Mdouble a, const SmallMatrix& x)
267  {
268  for (unsigned int i = 0; i < numberOfRows * numberOfColumns; ++i)
269  {
270  data_[i] += a * x[i];
271  }
272  }
273 
275  unsigned int size() const
276  {
277  return numberOfRows * numberOfColumns;
278  }
279 
281  unsigned int getNumberOfRows() const
282  {
283  return numberOfRows;
284  }
285 
287  unsigned int getNRows() const
288  {
289  return getNumberOfRows();
290  }
291 
293  unsigned int getNumberOfColumns() const
294  {
295  return numberOfColumns;
296  }
297 
298  unsigned int getNCols() const
299  {
300  return getNumberOfColumns();
301  }
302 
304  SmallVector<numberOfRows> getColumn(unsigned int j) const
305  {
306  logger.assert(j < numberOfColumns, "Asked for column %, but there are only % columns", j, numberOfColumns);
307  return SmallVector<numberOfRows>(data() + j * numberOfRows);
308  }
309 
312  {
313  logger.assert(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  }
321 
324 
325  Mdouble determinant() const;
326 
328  SmallMatrix inverse() const;
329 
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  }
342 
344  template<unsigned int numberOfRightHandSideColumns>
346 
349  void solve(SmallVector<numberOfRows>& b) const;
350 
352  {
353  return data_.data();
354  }
355 
356  const Mdouble* data() const
357  {
358  return data_.data();
359  }
360 
361 private:
363  std::array<Mdouble, numberOfRows * numberOfColumns> data_;
364 };
365 
367 template<unsigned int numberOfRows, unsigned int numberOfColumns>
368 std::ostream& operator<<(std::ostream& os, const SmallMatrix<numberOfRows, numberOfColumns>& A)
369 {
370  for (unsigned int i = 0; i < numberOfRows; ++i)
371  {
372  os << A.getRow(i) << std::endl;
373  }
374  return os;
375 }
376 
378 template<unsigned int numberOfRows, unsigned int numberOfColumns>
381 {
382  return mat * d;
383 }
384 
386 template<unsigned int numberOfRows, unsigned int numberOfColumns>
388 
389 #include "SmallMatrix_impl.h"
390 
391 
392 #endif //MERCURY_SMALLMATRIX_H
SmallMatrix(const SmallVector< numberOfRows > &other)
Definition: SmallMatrix.h:77
SmallMatrix operator*(const Mdouble &right) const
Does matrix A_ij=scalar*B_ij.
Definition: SmallMatrix.h:173
Mdouble * data()
Definition: SmallMatrix.h:351
SmallMatrix & operator=(SmallMatrix &&right)
Assigns one matrix to another.
Definition: SmallMatrix.h:256
const Mdouble & operator()(unsigned int n, unsigned int m) const
defines the operator(n,m) to access the element on row n and column m
Definition: SmallMatrix.h:144
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
SmallMatrix(SmallMatrix &&other)
Move Matrix from another Matrix.
Definition: SmallMatrix.h:129
const Mdouble * data() const
Definition: SmallMatrix.h:356
SmallMatrix(std::array< SmallVector< numberOfRows >, numberOfColumns > entries)
Glues one or more vectors with the same number of rows together.
Definition: SmallMatrix.h:116
SmallMatrix operator-(const SmallMatrix &other) const
Definition: SmallMatrix.h:207
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:330
unsigned int size() const
Get total number of Matrix entries.
Definition: SmallMatrix.h:275
Mdouble & operator[](const unsigned int n)
Access the n linear element in the matrix.
Definition: SmallMatrix.h:153
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
Mdouble & operator()(unsigned int n, unsigned int m)
defines the operator(n,m) to access the element on row n and column m
Definition: SmallMatrix.h:135
void axpy(Mdouble a, const SmallMatrix &x)
Applies the matrix y=ax + y, where x is another matrix and a is a scalar.
Definition: SmallMatrix.h:266
unsigned int getNumberOfColumns() const
Get the number of columns.
Definition: SmallMatrix.h:293
unsigned int getNCols() const
Definition: SmallMatrix.h:298
SmallVector< numberOfColumns > getRow(unsigned int i) const
get the i^th row
Definition: SmallMatrix.h:311
unsigned int getNumberOfRows() const
Get the number of rows.
Definition: SmallMatrix.h:281
std::array< Mdouble, numberOfRows *numberOfColumns > data_
The actually data of the matrix class.
Definition: SmallMatrix.h:363
SmallMatrix & operator/=(const Mdouble &scalar)
Does matrix A_ij=scalar*A_ij.
Definition: SmallMatrix.h:232
Mdouble determinant() const
const Mdouble & operator[](const unsigned int n) const
Definition: SmallMatrix.h:160
SmallMatrix & operator+=(const SmallMatrix &other)
Definition: SmallMatrix.h:188
const Mdouble * data() const
Definition: SmallVector.h:238
unsigned int getNRows() const
Definition: SmallMatrix.h:287
SmallVector< numberOfRows > getColumn(unsigned int j) const
get the j^th column
Definition: SmallMatrix.h:304
SmallMatrix< numberOfRows, numberOfColumns > operator*(const Mdouble d, const SmallMatrix< numberOfRows, numberOfColumns > &mat)
Multiplies a matrix with a Mdouble.
Definition: SmallMatrix.h:380
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...
SmallMatrix operator-() const
Definition: SmallMatrix.h:214
SmallMatrix & operator=(const SmallMatrix &right)
Assigns one matrix to another.
Definition: SmallMatrix.h:249
SmallMatrix inverse() const
return the inverse in the vector result. The size of result matches the matrix.
SmallMatrix()
Constructs a matrix of size n-rows by m-columns.
Definition: SmallMatrix.h:72
SmallMatrix & operator-=(const SmallMatrix &other)
Definition: SmallMatrix.h:194
SmallVector< numberOfRows > computeWedgeStuffVector() const
computeWedgeStuffVector.
SmallMatrix LUfactorisation() const
Return the LUfactorisation of the matrix.
SmallMatrix(const SmallMatrix &other)
Construct and copy Matrix from another Matrix i.e. B(A) where B and A are both matrices.
Definition: SmallMatrix.h:109
SmallMatrix & operator*=(const Mdouble &scalar)
Does matrix A_ij=scalar*A_ij.
Definition: SmallMatrix.h:220
SmallMatrix(const Mdouble &c)
Constructs a matrix of size n-rows by m-columns and initialises all entry to a constant.
Definition: SmallMatrix.h:85
Data type for small dense matrix.
Definition: SmallMatrix.h:67
SmallMatrix operator+(const SmallMatrix &other) const
Definition: SmallMatrix.h:200
SmallMatrix(const std::initializer_list< SmallVector< numberOfRows >> &entries)
Definition: SmallMatrix.h:91
SmallVector< numberOfRows > operator*(SmallVector< numberOfColumns > &right)
Defines Matrix A times vector B and return vector C i.e. C_,j= A_ij B_,j.
SmallMatrix operator/(const Mdouble &scalar) const
this does element by divided by a scalar
Definition: SmallMatrix.h:240