Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MatrixOperations.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MatrixOperations.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Defines the global functions needed for matrix operations.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/LinearAlgebra/MatrixOperations.hpp>
37 #include <LibUtilities/LinearAlgebra/ExplicitInstantiation.h>
38 
39 namespace Nektar
40 {
41  template<typename ResultDataType, typename LhsDataType, typename LhsMatrixType>
42  void Multiply(NekMatrix<ResultDataType, StandardMatrixTag>& result,
43  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
44  const ResultDataType& rhs)
45  {
46  // TODO - optimize for the different matrix types.
47  int n = lhs.GetRows();
48  int m = lhs.GetColumns();
49  for(unsigned int i = 0; i < n; ++i)
50  {
51  for(unsigned int j = 0; j < m; ++j)
52  {
53  result(i,j) = lhs(i,j)*rhs;
54  }
55  }
56  }
57 
58  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (1, (const NekDouble&)))
59 
60  template<typename DataType, typename LhsDataType, typename LhsMatrixType>
61  NekMatrix<typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType, StandardMatrixTag>
62  Multiply(const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
63  const DataType& rhs)
64  {
65  typedef typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType ResultDataType;
66  NekMatrix<ResultDataType, StandardMatrixTag> result(lhs.GetRows(), lhs.GetColumns());
67  Multiply(result, lhs, rhs);
68  return result;
69  }
70 
71  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (DNekMat)), (0, ()), (1, (const NekDouble&)))
72 
73  template<typename RhsDataType, typename RhsMatrixType, typename ResultDataType>
74  void Multiply(NekMatrix<ResultDataType, StandardMatrixTag>& result,
75  const ResultDataType& lhs,
76  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
77 
78  {
79  Multiply(result, rhs, lhs);
80  }
81 
82  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (2, (DNekMat&, const NekDouble&)), (0, ()))
83 
84  template<typename DataType, typename RhsDataType, typename RhsMatrixType>
85  NekMatrix<typename NekMatrix<RhsDataType, RhsMatrixType>::NumberType, StandardMatrixTag>
86  Multiply(const DataType& lhs,
87  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
88 
89  {
90  return Multiply(rhs, lhs);
91  }
92 
93  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(Multiply, NEKTAR_ALL_MATRIX_TYPES, (1, (DNekMat)), (1, (const NekDouble&)), (0, ()))
94 
95  template<typename LhsDataType>
96  void MultiplyEqual(NekMatrix<LhsDataType, StandardMatrixTag>& lhs,
97  typename boost::call_traits<LhsDataType>::const_reference rhs)
98  {
99  typedef typename NekMatrix<LhsDataType, StandardMatrixTag>::iterator iterator;
100  for( iterator iter = lhs.begin(); iter != lhs.end(); ++iter)
101  {
102  *iter *= rhs;
103  }
104  }
105 
106  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(MultiplyEqual, (1, (DNekMat&)), (1, (void)), (0, ()), (1, (const NekDouble&)))
107 
108  template<typename ResultType, typename LhsDataType, typename RhsDataType,
109  typename LhsMatrixType, typename RhsMatrixType>
110  void NekMultiplyDefaultImpl(NekMatrix<ResultType, StandardMatrixTag>& result,
111  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
112  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
113  {
114  ASSERTL1(lhs.GetColumns() == rhs.GetRows(), std::string("A left side matrix with column count ") +
115  boost::lexical_cast<std::string>(lhs.GetColumns()) +
116  std::string(" and a right side matrix with row count ") +
117  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be multiplied."));
118 
119  for(unsigned int i = 0; i < result.GetRows(); ++i)
120  {
121  for(unsigned int j = 0; j < result.GetColumns(); ++j)
122  {
123  ResultType t = ResultType(0);
124 
125  // Set the result(i,j) element.
126  for(unsigned int k = 0; k < lhs.GetColumns(); ++k)
127  {
128  t += lhs(i,k)*rhs(k,j);
129  }
130  result(i,j) = t;
131  }
132  }
133  }
134 
135  template<typename ResultType, typename LhsDataType, typename RhsDataType,
136  typename LhsMatrixType, typename RhsMatrixType>
137  void NekMultiplyFullMatrixFullMatrix(NekMatrix<ResultType, StandardMatrixTag>& result,
138  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
139  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
140  {
141  NekMultiplyDefaultImpl(result, lhs, rhs);
142  }
143 
144  template<typename LhsDataType, typename RhsDataType, typename DataType,
145  typename LhsMatrixType, typename RhsMatrixType>
146  void Multiply(NekMatrix<DataType, StandardMatrixTag>& result,
147  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
148  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
149  {
150  ASSERTL1(lhs.GetColumns() == rhs.GetRows(), std::string("A left side matrix with column count ") +
151  boost::lexical_cast<std::string>(lhs.GetColumns()) +
152  std::string(" and a right side matrix with row count ") +
153  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be multiplied."));
154 
155  result.SetSize(lhs.GetRows(), rhs.GetColumns());
156  if( lhs.GetType() == eFULL && rhs.GetType() == eFULL)
157  {
158  NekMultiplyFullMatrixFullMatrix(result, lhs, rhs);
159  }
160  else
161  {
162  NekMultiplyDefaultImpl(result, lhs, rhs);
163  }
164  }
165 
166  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(Multiply, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
167 
168  template<typename DataType, typename RhsDataType, typename RhsMatrixType>
169  void AddEqual(NekMatrix<DataType, StandardMatrixTag>& result,
170  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
171  {
172  ASSERTL1(result.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
173  boost::lexical_cast<std::string>(result.GetRows()) + std::string(" and ") +
174  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be added."));
175  ASSERTL1(result.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
176  boost::lexical_cast<std::string>(result.GetColumns()) + std::string(" and ") +
177  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be added."));
178 
179  for(unsigned int i = 0; i < rhs.GetRows(); ++i)
180  {
181  for(unsigned int j = 0; j < rhs.GetColumns(); ++j)
182  {
183  result(i,j) += rhs(i,j);
184  }
185  }
186  }
187 
188 
189  template<typename DataType, typename RhsDataType, typename RhsMatrixType>
190  void AddEqualNegatedLhs(NekMatrix<DataType, StandardMatrixTag>& result,
191  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
192  {
193  ASSERTL1(result.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
194  boost::lexical_cast<std::string>(result.GetRows()) + std::string(" and ") +
195  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be added."));
196  ASSERTL1(result.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
197  boost::lexical_cast<std::string>(result.GetColumns()) + std::string(" and ") +
198  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be added."));
199 
200  for(unsigned int i = 0; i < rhs.GetRows(); ++i)
201  {
202  for(unsigned int j = 0; j < rhs.GetColumns(); ++j)
203  {
204  result(i,j) = -result(i,j) + rhs(i,j);
205  }
206  }
207  }
208 
209  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(AddEqual, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
210  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(AddEqualNegatedLhs, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
211 
212 
213  template<typename DataType, typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
214  void Add(NekMatrix<DataType, StandardMatrixTag>& result,
215  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
216  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
217  {
218  ASSERTL1(lhs.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
219  boost::lexical_cast<std::string>(lhs.GetRows()) + std::string(" and ") +
220  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be added."));
221  ASSERTL1(lhs.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
222  boost::lexical_cast<std::string>(lhs.GetColumns()) + std::string(" and ") +
223  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be added."));
224 
225  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
226  {
227  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
228  {
229  result(i,j) = lhs(i,j) + rhs(i,j);
230  }
231  }
232  }
233 
234  template<typename DataType, typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
235  void AddNegatedLhs(NekMatrix<DataType, StandardMatrixTag>& result,
236  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
237  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
238  {
239  ASSERTL1(lhs.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
240  boost::lexical_cast<std::string>(lhs.GetRows()) + std::string(" and ") +
241  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be added."));
242  ASSERTL1(lhs.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
243  boost::lexical_cast<std::string>(lhs.GetColumns()) + std::string(" and ") +
244  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be added."));
245 
246  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
247  {
248  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
249  {
250  result(i,j) = -lhs(i,j) + rhs(i,j);
251  }
252  }
253  }
254 
255  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(Add, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
256  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(AddNegatedLhs, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
257 
258 
259  template<typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
260  NekMatrix<typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType, StandardMatrixTag>
261  Add(const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
262  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
263  {
264  typedef typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType NumberType;
265  NekMatrix<NumberType, StandardMatrixTag> result(lhs.GetRows(), lhs.GetColumns());
266  Add(result, lhs, rhs);
267  return result;
268  }
269 
270  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(Add, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (DNekMat)), (0, ()), (0, ()))
271 
272  template<typename DataType, typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
273  void Subtract(NekMatrix<DataType, StandardMatrixTag>& result,
274  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
275  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
276  {
277  ASSERTL1(lhs.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
278  boost::lexical_cast<std::string>(lhs.GetRows()) + std::string(" and ") +
279  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be subtracted."));
280  ASSERTL1(lhs.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
281  boost::lexical_cast<std::string>(lhs.GetColumns()) + std::string(" and ") +
282  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be subtracted."));
283 
284  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
285  {
286  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
287  {
288  result(i,j) = lhs(i,j) - rhs(i,j);
289  }
290  }
291  }
292 
293  template<typename DataType, typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
294  void SubtractNegatedLhs(NekMatrix<DataType, StandardMatrixTag>& result,
295  const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
296  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
297  {
298  ASSERTL1(lhs.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
299  boost::lexical_cast<std::string>(lhs.GetRows()) + std::string(" and ") +
300  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be subtracted."));
301  ASSERTL1(lhs.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
302  boost::lexical_cast<std::string>(lhs.GetColumns()) + std::string(" and ") +
303  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be subtracted."));
304 
305  for(unsigned int i = 0; i < lhs.GetRows(); ++i)
306  {
307  for(unsigned int j = 0; j < lhs.GetColumns(); ++j)
308  {
309  result(i,j) = -lhs(i,j) - rhs(i,j);
310  }
311  }
312  }
313 
314  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(Subtract, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
315  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(SubtractNegatedLhs, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
316 
317 
318 
319  template<typename DataType, typename RhsDataType, typename RhsMatrixType>
320  void SubtractEqual(NekMatrix<DataType, StandardMatrixTag>& result,
321  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
322  {
323  ASSERTL1(result.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
324  boost::lexical_cast<std::string>(result.GetRows()) + std::string(" and ") +
325  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be subtracted."));
326  ASSERTL1(result.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
327  boost::lexical_cast<std::string>(result.GetColumns()) + std::string(" and ") +
328  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be subtracted."));
329 
330  for(unsigned int i = 0; i < rhs.GetRows(); ++i)
331  {
332  for(unsigned int j = 0; j < rhs.GetColumns(); ++j)
333  {
334  result(i,j) -= rhs(i,j);
335  }
336  }
337  }
338 
339  template<typename DataType, typename RhsDataType, typename RhsMatrixType>
340  void SubtractEqualNegatedLhs(NekMatrix<DataType, StandardMatrixTag>& result,
341  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
342  {
343  ASSERTL1(result.GetRows() == rhs.GetRows(), std::string("Matrices with different row counts ") +
344  boost::lexical_cast<std::string>(result.GetRows()) + std::string(" and ") +
345  boost::lexical_cast<std::string>(rhs.GetRows()) + std::string(" can't be subtracted."));
346  ASSERTL1(result.GetColumns() == rhs.GetColumns(), std::string("Matrices with different column counts ") +
347  boost::lexical_cast<std::string>(result.GetColumns()) + std::string(" and ") +
348  boost::lexical_cast<std::string>(rhs.GetColumns()) + std::string(" can't be subtracted."));
349 
350  for(unsigned int i = 0; i < rhs.GetRows(); ++i)
351  {
352  for(unsigned int j = 0; j < rhs.GetColumns(); ++j)
353  {
354  result(i,j) = -result(i,j) - rhs(i,j);
355  }
356  }
357  }
358 
359  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(SubtractEqual, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
360  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_SINGLE_MATRIX(SubtractEqualNegatedLhs, NEKTAR_ALL_MATRIX_TYPES, (1, (void)), (1, (DNekMat&)), (0, ()))
361 
362  template<typename LhsDataType, typename LhsMatrixType, typename RhsDataType, typename RhsMatrixType>
363  NekMatrix<typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType, StandardMatrixTag>
364  Subtract(const NekMatrix<LhsDataType, LhsMatrixType>& lhs,
365  const NekMatrix<RhsDataType, RhsMatrixType>& rhs)
366  {
367  typedef typename NekMatrix<LhsDataType, LhsMatrixType>::NumberType NumberType;
368  NekMatrix<NumberType, StandardMatrixTag> result(lhs.GetRows(), lhs.GetColumns());
369  Subtract(result, lhs, rhs);
370  return result;
371  }
372 
373  NEKTAR_GENERATE_EXPLICIT_FUNCTION_INSTANTIATION_TWO_MATRICES(Subtract, NEKTAR_ALL_MATRIX_TYPES, NEKTAR_ALL_MATRIX_TYPES, (1, (DNekMat)), (0, ()), (0, ()))
374 }
375 
376 
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2