Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StandardMatrix.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // For more information, please see: http://www.nektar.info
4 //
5 // The MIT License
6 //
7 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
8 // Department of Aeronautics, Imperial College London (UK), and Scientific
9 // Computing and Imaging Institute, University of Utah (USA).
10 //
11 // License for the specific language governing rights and limitations under
12 // Permission is hereby granted, free of charge, to any person obtaining a
13 // copy of this software and associated documentation files (the "Software"),
14 // to deal in the Software without restriction, including without limitation
15 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 // and/or sell copies of the Software, and to permit persons to whom the
17 // Software is furnished to do so, subject to the following conditions:
18 //
19 // The above copyright notice and this permission notice shall be included
20 // in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 // DEALINGS IN THE SOFTWARE.
29 //
30 ///////////////////////////////////////////////////////////////////////////////
31 
32 #include <LibUtilities/LinearAlgebra/StandardMatrix.hpp>
33 
34 
35 namespace Nektar
36 {
37 
38 
39  template<typename DataType>
40  NekMatrix<DataType, StandardMatrixTag>::NekMatrix() :
41  Matrix<DataType>(0, 0),
42  m_data(),
43  m_wrapperType(eCopy),
44  m_storagePolicy(eFULL),
45  m_numberOfSuperDiagonals(std::numeric_limits<unsigned int>::max()),
46  m_numberOfSubDiagonals(std::numeric_limits<unsigned int>::max()),
47  m_tempSpace()
48  {
49  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
50  }
51 
52  template<typename DataType>
53  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, MatrixStorage policy,
54  unsigned int subDiagonals,
55  unsigned int superDiagonals) :
56  Matrix<DataType>(rows, columns),
57  m_data(),
58  m_wrapperType(eCopy),
59  m_storagePolicy(policy),
60  m_numberOfSuperDiagonals(superDiagonals),
61  m_numberOfSubDiagonals(subDiagonals),
62  m_tempSpace()
63  {
64  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
65  }
66 
67  template<typename DataType>
68  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns,
69  MatrixStorage policy,
70  unsigned int subDiagonals,
71  unsigned int superDiagonals,
72  unsigned int capacity) :
73  Matrix<DataType>(rows, columns),
74  m_data(),
75  m_wrapperType(eCopy),
76  m_storagePolicy(policy),
77  m_numberOfSuperDiagonals(superDiagonals),
78  m_numberOfSubDiagonals(subDiagonals),
79  m_tempSpace()
80  {
81  unsigned int requiredStorage = this->GetRequiredStorageSize();
82  unsigned int actualSize = std::max(requiredStorage, capacity);
83  m_data = Array<OneD, DataType>(actualSize);
84  }
85 
86  template<typename DataType>
87  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, typename boost::call_traits<DataType>::const_reference initValue,
88  MatrixStorage policy,
89  unsigned int subDiagonals,
90  unsigned int superDiagonals) :
91  Matrix<DataType>(rows, columns),
92  m_data(),
93  m_wrapperType(eCopy),
94  m_storagePolicy(policy),
95  m_numberOfSuperDiagonals(superDiagonals),
96  m_numberOfSubDiagonals(subDiagonals),
97  m_tempSpace()
98  {
99  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
100  std::fill(m_data.begin(), m_data.end(), initValue);
101  }
102 
103  template<typename DataType>
104  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const DataType* data,
105  MatrixStorage policy,
106  unsigned int subDiagonals,
107  unsigned int superDiagonals) :
108  Matrix<DataType>(rows, columns),
109  m_data(),
110  m_wrapperType(eCopy),
111  m_storagePolicy(policy),
112  m_numberOfSuperDiagonals(superDiagonals),
113  m_numberOfSubDiagonals(subDiagonals),
114  m_tempSpace()
115  {
116  unsigned int size = GetRequiredStorageSize();
117  m_data = Array<OneD, DataType>(size);
118  std::copy(data, data + size, m_data.begin());
119  }
120 
121  template<typename DataType>
122  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const Array<OneD, const DataType>& d,
123  MatrixStorage policy,
124  unsigned int subDiagonals,
125  unsigned int superDiagonals) :
126  Matrix<DataType>(rows, columns),
127  m_data(),
128  m_wrapperType(eCopy),
129  m_storagePolicy(policy),
130  m_numberOfSuperDiagonals(superDiagonals),
131  m_numberOfSubDiagonals(subDiagonals),
132  m_tempSpace()
133  {
134  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
135  CopyArrayN(d, m_data, m_data.num_elements());
136  }
137 
138  template<typename DataType>
139  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const Array<OneD, DataType>& d, PointerWrapper wrapperType,
140  MatrixStorage policy,
141  unsigned int subDiagonals,
142  unsigned int superDiagonals) :
143  Matrix<DataType>(rows, columns),
144  m_data(),
145  m_wrapperType(wrapperType),
146  m_storagePolicy(policy),
147  m_numberOfSuperDiagonals(superDiagonals),
148  m_numberOfSubDiagonals(subDiagonals),
149  m_tempSpace()
150  {
151  if( wrapperType == eWrapper )
152  {
153  m_data = Array<OneD, DataType>(d, eVECTOR_WRAPPER);
154  }
155  else
156  {
157  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
158  CopyArrayN(d, m_data, m_data.num_elements());
159  }
160  }
161 
162  template<typename DataType>
163  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(const NekMatrix<DataType, StandardMatrixTag>& rhs) :
164  Matrix<DataType>(rhs),
165  m_data(),
166  m_wrapperType(rhs.m_wrapperType),
167  m_storagePolicy(rhs.m_storagePolicy),
168  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
169  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
170  m_tempSpace()
171  {
172  PerformCopyConstruction(rhs);
173  }
174 
175 
176  template<typename DataType>
177  MatrixStorage NekMatrix<DataType, StandardMatrixTag>::GetType() const { return m_storagePolicy; }
178 
179  template<typename DataType>
180  NekMatrix<DataType, StandardMatrixTag>& NekMatrix<DataType, StandardMatrixTag>::operator=(const NekMatrix<DataType, StandardMatrixTag>& rhs)
181  {
182  if( this == &rhs )
183  {
184  return *this;
185  }
186 
187  Matrix<DataType>::operator=(rhs);
188  m_storagePolicy = rhs.m_storagePolicy;
189  m_numberOfSubDiagonals = rhs.m_numberOfSubDiagonals;
190  m_numberOfSuperDiagonals = rhs.m_numberOfSuperDiagonals;
191 
192  ResizeDataArrayIfNeeded();
193 
194  unsigned int requiredStorageSize = GetRequiredStorageSize();
195  std::copy(rhs.m_data.data(), rhs.m_data.data() + requiredStorageSize, m_data.data());
196 
197  return *this;
198  }
199 
200 
201  template<typename DataType>
202  typename NekMatrix<DataType, StandardMatrixTag>::ConstGetValueType NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column) const
203  {
204  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
205  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
206  std::string(" rows"));
207  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
208  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
209  std::string(" columns"));
210 
211  return this->GetValue(row, column, this->GetTransposeFlag());
212  }
213 
214  template<typename DataType>
215  typename NekMatrix<DataType, StandardMatrixTag>::ConstGetValueType NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column, char transpose) const
216  {
217  return this->GetValue(row, column, transpose);
218  }
219 
220  template<typename DataType>
221  unsigned int NekMatrix<DataType, StandardMatrixTag>::GetRequiredStorageSize() const
222  {
223  return Matrix<DataType>::GetRequiredStorageSize(this->GetStorageType(),
224  this->GetRows(), this->GetColumns(),
225  this->GetNumberOfSubDiagonals(), this->GetNumberOfSuperDiagonals());
226  }
227 
228  template<typename DataType>
229  unsigned int NekMatrix<DataType, StandardMatrixTag>::CalculateIndex(unsigned int row, unsigned int col, const char transpose) const
230  {
231  unsigned int numRows = this->GetSize()[0];
232  unsigned int numColumns = this->GetSize()[1];
233  return Matrix<DataType>::CalculateIndex(this->GetStorageType(),
234  row, col, numRows, numColumns, transpose,
235  m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
236  }
237 
238  template<typename DataType>
239  typename boost::call_traits<DataType>::const_reference NekMatrix<DataType, StandardMatrixTag>::GetValue(unsigned int row, unsigned int column) const
240  {
241  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
242  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
243  std::string(" rows"));
244  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
245  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
246  std::string(" columns"));
247 
248  return GetValue(row, column, this->GetTransposeFlag());
249  }
250 
251  template<typename DataType>
252  typename NekMatrix<DataType, StandardMatrixTag>::ConstGetValueType NekMatrix<DataType, StandardMatrixTag>::GetValue(unsigned int row, unsigned int column, char transpose) const
253  {
254  static DataType defaultReturnValue;
255  unsigned int index = CalculateIndex(row, column, transpose);
256  if( index != std::numeric_limits<unsigned int>::max() )
257  {
258  return m_data[index];
259  }
260  else
261  {
262  return defaultReturnValue;
263  }
264  }
265 
266  template<typename DataType>
267  const Array<OneD, const DataType>& NekMatrix<DataType, StandardMatrixTag>::GetPtr() const
268  {
269  return m_data;
270  }
271 
272  template<typename DataType>
273  DataType NekMatrix<DataType, StandardMatrixTag>::Scale() const
274  {
275  return DataType(1);
276  }
277 
278  template<typename DataType>
279  const DataType* NekMatrix<DataType, StandardMatrixTag>::GetRawPtr() const
280  {
281  return m_data.data();
282  }
283 
284  template<typename DataType>
285  typename NekMatrix<DataType, StandardMatrixTag>::const_iterator NekMatrix<DataType, StandardMatrixTag>::begin() const
286  {
287  return begin(this->GetTransposeFlag());
288  }
289 
290  template<typename DataType>
291  typename NekMatrix<DataType, StandardMatrixTag>::const_iterator NekMatrix<DataType, StandardMatrixTag>::begin(char transpose) const
292  {
293  if( transpose == 'N' )
294  {
295  return const_iterator(m_data.data(), m_data.data() + m_data.num_elements());
296  }
297  else
298  {
299  return const_iterator(this, transpose);
300  }
301  }
302 
303  template<typename DataType>
304  typename NekMatrix<DataType, StandardMatrixTag>::const_iterator NekMatrix<DataType, StandardMatrixTag>::end() const
305  {
306  return end(this->GetTransposeFlag());
307  }
308 
309  template<typename DataType>
310  typename NekMatrix<DataType, StandardMatrixTag>::const_iterator NekMatrix<DataType, StandardMatrixTag>::end(char transpose) const
311  {
312  if( transpose == 'N' )
313  {
314  return const_iterator(m_data.data(), m_data.data() + m_data.num_elements(), true);
315  }
316  else
317  {
318  return const_iterator(this, transpose, true);
319  }
320  }
321 
322  template<typename DataType>
323  unsigned int NekMatrix<DataType, StandardMatrixTag>::GetStorageSize() const
324  {
325  return m_data.num_elements();
326  }
327 
328  template<typename DataType>
329  unsigned int NekMatrix<DataType, StandardMatrixTag>::GetNumberOfSubDiagonals() const
330  {
331  if( m_numberOfSubDiagonals != std::numeric_limits<unsigned int>::max() )
332  {
333  return m_numberOfSubDiagonals;
334  }
335  else if( this->GetRows() > 0 )
336  {
337  return this->GetRows()-1;
338  }
339  else
340  {
341  return 0;
342  }
343  }
344 
345  template<typename DataType>
346  unsigned int NekMatrix<DataType, StandardMatrixTag>::GetNumberOfSuperDiagonals() const
347  {
348  if( m_numberOfSuperDiagonals != std::numeric_limits<unsigned int>::max() )
349  {
350  return m_numberOfSuperDiagonals;
351  }
352  else if( this->GetRows() > 0 )
353  {
354  return this->GetRows()-1;
355  }
356  else
357  {
358  return 0;
359  }
360  }
361 
362  template<typename DataType>
363  unsigned int NekMatrix<DataType, StandardMatrixTag>::CalculateNumberOfRows() const
364  {
365  return GetNumberOfSubDiagonals() + GetNumberOfSuperDiagonals() + 1;
366  }
367 
368  template<typename DataType>
369  bool NekMatrix<DataType, StandardMatrixTag>::operator==(const NekMatrix<DataType, StandardMatrixTag>& rhs) const
370  {
371  if( this->GetRows() != rhs.GetRows() ||
372  this->GetColumns() != rhs.GetColumns() )
373  {
374  return false;
375  }
376 
377  if( this->GetTransposeFlag() == rhs.GetTransposeFlag() )
378  {
379  return std::equal(begin(), end(), rhs.begin());
380  }
381  else
382  {
383  for(unsigned int i = 0; i < this->GetRows(); ++i)
384  {
385  for(unsigned int j = 0; j < this->GetColumns(); ++j)
386  {
387  if( (*this)(i,j) != rhs(i,j) )
388  {
389  return false;
390  }
391  }
392  }
393  }
394 
395  return true;
396  }
397 
398  template<typename DataType>
399  PointerWrapper NekMatrix<DataType, StandardMatrixTag>::GetWrapperType() const { return m_wrapperType; }
400 
401  template<typename DataType>
402  char NekMatrix<DataType, StandardMatrixTag>::GetTransposeFlag() const
403  {
404  return this->GetRawTransposeFlag();
405  }
406 
407  template<typename DataType>
408  boost::tuples::tuple<unsigned int, unsigned int>
409  NekMatrix<DataType, StandardMatrixTag>::Advance(unsigned int curRow, unsigned int curColumn) const
410  {
411  return Advance(curRow, curColumn, this->GetTransposeFlag());
412  }
413 
414  template<typename DataType>
415  boost::tuples::tuple<unsigned int, unsigned int>
416  NekMatrix<DataType, StandardMatrixTag>::Advance(unsigned int curRow, unsigned int curColumn, char transpose) const
417  {
418  unsigned int numRows = this->GetTransposedRows(transpose);
419  unsigned int numColumns = this->GetTransposedColumns(transpose);
420 
421  switch(m_storagePolicy)
422  {
423  case eFULL:
424  return FullMatrixFuncs::Advance(
425  numRows, numColumns, curRow, curColumn);
426  break;
427  case eDIAGONAL:
428  return DiagonalMatrixFuncs::Advance(
429  numRows, numColumns, curRow, curColumn);
430  break;
431  case eUPPER_TRIANGULAR:
432  return UpperTriangularMatrixFuncs::Advance(
433  numRows, numColumns, curRow, curColumn);
434  break;
435 
436  case eLOWER_TRIANGULAR:
437  return LowerTriangularMatrixFuncs::Advance(
438  numRows, numColumns, curRow, curColumn);
439  break;
440 
441  case eSYMMETRIC:
442  case ePOSITIVE_DEFINITE_SYMMETRIC:
443  return SymmetricMatrixFuncs::Advance(
444  numRows, numColumns, curRow, curColumn);
445  break;
446  case eBANDED:
447  return BandedMatrixFuncs::Advance(
448  numRows, numColumns, curRow, curColumn);
449  break;
450  case eSYMMETRIC_BANDED:
451  case ePOSITIVE_DEFINITE_SYMMETRIC_BANDED:
452  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
453  break;
454  case eUPPER_TRIANGULAR_BANDED:
455  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
456  break;
457  case eLOWER_TRIANGULAR_BANDED:
458  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
459  break;
460 
461  default:
462  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
463  }
464  return boost::tuples::tuple<unsigned int, unsigned int>(curRow, curColumn);
465  }
466 
467  template<typename DataType>
468  NekMatrix<DataType, StandardMatrixTag> NekMatrix<DataType, StandardMatrixTag>::CreateWrapper(const NekMatrix<DataType, StandardMatrixTag>& rhs)
469  {
470  return NekMatrix<DataType, StandardMatrixTag>(rhs, eWrapper);
471  }
472 
473  template<typename DataType>
474  boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> > NekMatrix<DataType, StandardMatrixTag>::CreateWrapper(const boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >& rhs)
475  {
476  return boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >(new NekMatrix<DataType, StandardMatrixTag>(*rhs, eWrapper));
477  }
478 
479  template<typename DataType>
480  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(const NekMatrix<DataType, StandardMatrixTag>& rhs, PointerWrapper wrapperType) :
481  BaseType(rhs),
482  m_data(),
483  m_wrapperType(wrapperType),
484  m_storagePolicy(rhs.m_storagePolicy),
485  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
486  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals)
487  {
488  PerformCopyConstruction(rhs);
489  }
490 
491 
492  template<typename DataType>
493  Array<OneD, DataType>& NekMatrix<DataType, StandardMatrixTag>::GetData() { return m_data; }
494 
495  template<typename DataType>
496  void NekMatrix<DataType, StandardMatrixTag>::RemoveExcessCapacity()
497  {
498  if( m_wrapperType == eCopy )
499  {
500  unsigned int requiredStorageSize = GetRequiredStorageSize();
501  if( m_data.num_elements() > requiredStorageSize )
502  {
503  Array<OneD, DataType> newArray(requiredStorageSize);
504  CopyArrayN(m_data, newArray, requiredStorageSize);
505  m_data = newArray;
506  }
507  }
508  else if( m_wrapperType == eWrapper )
509  {
510  ASSERTL0(true, "Can't call RemoveExcessCapacity on a wrapped matrix.");
511  }
512  }
513 
514  template<typename DataType>
515  void NekMatrix<DataType, StandardMatrixTag>::ResizeDataArrayIfNeeded(unsigned int requiredStorageSize)
516  {
517  if( m_wrapperType == eCopy )
518  {
519  if( m_data.num_elements() < requiredStorageSize )
520  {
521  Array<OneD, DataType> newData(requiredStorageSize);
522  std::copy(m_data.data(), m_data.data() + m_data.num_elements(), newData.data());
523  m_data = newData;
524  }
525  }
526  else if( m_wrapperType == eWrapper )
527  {
528  // If the current matrix is wrapped, then just copy over the top,
529  // but the sizes of the two matrices must be the same.
530  ASSERTL0(m_data.num_elements() >= requiredStorageSize, "Wrapped NekMatrices must have the same dimension in operator=");
531  }
532  }
533 
534  template<typename DataType>
535  void NekMatrix<DataType, StandardMatrixTag>::ResizeDataArrayIfNeeded()
536  {
537  unsigned int requiredStorageSize = GetRequiredStorageSize();
538  ResizeDataArrayIfNeeded(requiredStorageSize);
539  }
540 
541  template<typename DataType>
542  void NekMatrix<DataType, StandardMatrixTag>::PerformCopyConstruction(const ThisType& rhs)
543  {
544  if( m_wrapperType == eWrapper )
545  {
546  m_data = rhs.m_data;
547  }
548  else
549  {
550  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
551  CopyArrayN(rhs.m_data, m_data, m_data.num_elements());
552  }
553  }
554 
555  template<typename DataType>
556  typename boost::call_traits<DataType>::value_type NekMatrix<DataType, StandardMatrixTag>::v_GetValue(unsigned int row, unsigned int column) const
557  {
558  return NekMatrix<DataType, StandardMatrixTag>::operator()(row, column);
559  }
560 
561  template<typename DataType>
562  unsigned int NekMatrix<DataType, StandardMatrixTag>::v_GetStorageSize() const
563  {
564  return NekMatrix<DataType, StandardMatrixTag>::GetStorageSize();
565  }
566 
567  template<typename DataType>
568  MatrixStorage NekMatrix<DataType, StandardMatrixTag>::v_GetStorageType() const
569  {
570  return NekMatrix<DataType, StandardMatrixTag>::GetType();
571  }
572 
573  template<typename DataType>
574  void NekMatrix<DataType, StandardMatrixTag>::v_SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d)
575  {
576  return NekMatrix<DataType, StandardMatrixTag>::SetValue(row, column, d);
577  }
578 
579 
580 
581  template<typename DataType>
582  void NekMatrix<DataType, StandardMatrixTag>::SetSize(unsigned int rows, unsigned int cols)
583  {
584  this->Resize(rows, cols);
585 
586  // Some places in Nektar++ access the matrix data array and
587  // use num_elements() to see how big it is. When using
588  // expression templates, the data array's capacity is often larger
589  // than the actual number of elements, so this statement is
590  // required to report the correct number of elements.
591  this->GetData().ChangeSize(this->GetRequiredStorageSize());
592  ASSERTL0(this->GetRequiredStorageSize() <= this->GetData().num_elements(), "Can't resize matrices if there is not enough capacity.");
593  }
594 
595 
596  template<typename DataType>
597  typename NekMatrix<DataType, StandardMatrixTag>::Proxy NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column)
598  {
599  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
600  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
601  std::string(" rows"));
602  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
603  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
604  std::string(" columns"));
605 
606  return (*this)(row, column, this->GetTransposeFlag());
607  }
608 
609  template<typename DataType>
610  typename NekMatrix<DataType, StandardMatrixTag>::Proxy NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column, char transpose)
611  {
612  unsigned int index = this->CalculateIndex(row, column, transpose);
613  if( index != std::numeric_limits<unsigned int>::max() )
614  {
615  return Proxy(this->GetData()[index]);
616  }
617  else
618  {
619  return Proxy();
620  }
621 
622  }
623 
624  template<typename DataType>
625  void NekMatrix<DataType, StandardMatrixTag>::SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d)
626  {
627  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
628  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
629  std::string(" rows"));
630  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
631  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
632  std::string(" columns"));
633  SetValue(row, column, d, this->GetTransposeFlag());
634  }
635 
636  template<typename DataType>
637  void NekMatrix<DataType, StandardMatrixTag>::SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d, char transpose)
638  {
639  unsigned int index = this->CalculateIndex(row, column, transpose);
640  if( index != std::numeric_limits<unsigned int>::max() )
641  {
642  this->GetData()[index] = d;
643  }
644  else
645  {
646  NEKERROR(ErrorUtil::efatal, "Can't assign values into zeroed elements of a special array.");
647  }
648  }
649 
650  template<typename DataType>
651  Array<OneD, DataType>& NekMatrix<DataType, StandardMatrixTag>::GetPtr()
652  {
653  return this->GetData();
654  }
655 
656  template<typename DataType>
657  DataType* NekMatrix<DataType, StandardMatrixTag>::GetRawPtr()
658  {
659  return this->GetData().data();
660  }
661 
662  template<typename DataType>
663  typename NekMatrix<DataType, StandardMatrixTag>::iterator NekMatrix<DataType, StandardMatrixTag>::begin()
664  {
665  return begin(this->GetTransposeFlag());
666  }
667 
668  template<typename DataType>
669  typename NekMatrix<DataType, StandardMatrixTag>::iterator NekMatrix<DataType, StandardMatrixTag>::begin(char transpose)
670  {
671  if( transpose == 'N' )
672  {
673  return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements());
674  }
675  else
676  {
677  return iterator(this, transpose);
678  }
679  }
680 
681  template<typename DataType>
682  typename NekMatrix<DataType, StandardMatrixTag>::iterator NekMatrix<DataType, StandardMatrixTag>::end()
683  {
684  return end(this->GetTransposeFlag());
685  }
686 
687  template<typename DataType>
688  typename NekMatrix<DataType, StandardMatrixTag>::iterator NekMatrix<DataType, StandardMatrixTag>::end(char transpose)
689  {
690  if( transpose == 'N' )
691  {
692  return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements(), true);
693  }
694  else
695  {
696  return iterator(this, transpose, true);
697  }
698  }
699 
700  template<typename DataType>
701  NekDouble NekMatrix<DataType, StandardMatrixTag>::AbsMaxtoMinEigenValueRatio(void)
702  {
703  NekDouble returnval;
704  int nvals = this->GetColumns();
705  Array<OneD, NekDouble> EigValReal(nvals);
706  Array<OneD, NekDouble> EigValImag(nvals);
707 
708  EigenSolve(EigValReal,EigValImag);
709 
710  Vmath::Vmul(nvals,EigValReal,1,EigValReal,1,EigValReal,1);
711  Vmath::Vmul(nvals,EigValImag,1,EigValImag,1,EigValImag,1);
712  Vmath::Vadd(nvals,EigValReal,1,EigValImag,1,EigValReal,1);
713 
714  returnval = sqrt(Vmath::Vmax(nvals,EigValReal,1)/Vmath::Vmin(nvals,EigValReal,1));
715 
716  return returnval;
717  }
718 
719  template<typename DataType>
720  void NekMatrix<DataType, StandardMatrixTag>::EigenSolve(Array<OneD, NekDouble> &EigValReal,
721  Array<OneD, NekDouble> &EigValImag,
722  Array<OneD, NekDouble> &EigVecs)
723  {
724  ASSERTL0(this->GetRows()==this->GetColumns(), "Only square matrices can be called");
725 
726  switch(this->GetType())
727  {
728  case eFULL:
729  FullMatrixFuncs::EigenSolve(this->GetRows(),
730  this->GetData(), EigValReal,
731  EigValImag, EigVecs);
732  break;
733  case eDIAGONAL:
734  Vmath::Vcopy(this->GetRows(),&(this->GetData())[0],1, &EigValReal[0],1);
735  Vmath::Zero(this->GetRows(),&EigValImag[0],1);
736  break;
737  case eUPPER_TRIANGULAR:
738  case eLOWER_TRIANGULAR:
739  case eSYMMETRIC:
740  case eBANDED:
741  case eSYMMETRIC_BANDED:
742  case eUPPER_TRIANGULAR_BANDED:
743  case eLOWER_TRIANGULAR_BANDED:
744  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
745  break;
746 
747  default:
748  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
749  }
750  }
751 
752  template<typename DataType>
753  void NekMatrix<DataType, StandardMatrixTag>::Invert()
754  {
755  ASSERTL0(this->GetRows()==this->GetColumns(), "Only square matrices can be inverted.");
756  ASSERTL0(this->GetTransposeFlag()=='N', "Only untransposed matrices may be inverted.");
757 
758  switch(this->GetType())
759  {
760  case eFULL:
761  FullMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
762  this->GetData(), this->GetTransposeFlag());
763  break;
764  case eDIAGONAL:
765  DiagonalMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
766  this->GetData());
767  break;
768  case eUPPER_TRIANGULAR:
769  case eLOWER_TRIANGULAR:
770  case eSYMMETRIC:
771  case eBANDED:
772  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
773  break;
774  case eSYMMETRIC_BANDED:
775  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
776  break;
777  case eUPPER_TRIANGULAR_BANDED:
778  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
779  break;
780  case eLOWER_TRIANGULAR_BANDED:
781  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
782  break;
783 
784  default:
785  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
786  }
787  }
788 
789  template<typename DataType>
790  Array<OneD, DataType>& NekMatrix<DataType, StandardMatrixTag>::GetTempSpace()
791  {
792  if( m_tempSpace.capacity() == 0 )
793  {
794  m_tempSpace = Array<OneD, DataType>(this->GetData().capacity());
795  }
796  return m_tempSpace;
797  }
798 
799  template<typename DataType>
800  void NekMatrix<DataType, StandardMatrixTag>::SwapTempAndDataBuffers()
801  {
802  std::swap(m_tempSpace, this->GetData());
803  }
804 
805  #ifndef NEKTAR_USE_EXPRESSION_TEMPLATES
806  template<typename DataType>
807  NekMatrix<DataType, StandardMatrixTag> NekMatrix<DataType, StandardMatrixTag>::operator-() const
808  {
809  NekMatrix<DataType, StandardMatrixTag> result(*this);
810  NegateInPlace(result);
811  return result;
812  }
813  #endif
814 
815  template<typename DataType>
816  NekMatrix<DataType, StandardMatrixTag>& NekMatrix<DataType, StandardMatrixTag>::operator*=(const DataType& s)
817  {
818  for(unsigned int i = 0; i < this->GetPtr().num_elements(); ++i)
819  {
820  this->GetPtr()[i] *= s;
821  }
822  return *this;
823  }
824 
825 
826  template<typename DataType>
827  NekMatrix<DataType, StandardMatrixTag>
828  Transpose(NekMatrix<DataType, StandardMatrixTag>& rhs)
829  {
830  NekMatrix<DataType, StandardMatrixTag> result(rhs.GetRows(), rhs.GetColumns(),
831  rhs.GetPtr(), eWrapper, rhs.GetType(), rhs.GetNumberOfSubDiagonals(),
832  rhs.GetNumberOfSuperDiagonals());
833  result.Transpose();
834  return result;
835  }
836 
837  template LIB_UTILITIES_EXPORT NekMatrix<NekDouble, StandardMatrixTag> Transpose(NekMatrix<NekDouble, StandardMatrixTag>& rhs);
838 
839  template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>;
840 
841  template<typename DataType>
842  void NegateInPlace(NekMatrix<DataType, StandardMatrixTag>& m)
843  {
844  for(unsigned int i = 0; i < m.GetRows(); ++i)
845  {
846  for(unsigned int j = 0; j < m.GetColumns(); ++j)
847  {
848  m(i,j) *= -1.0;
849  }
850  }
851  }
852 
853  template LIB_UTILITIES_EXPORT void NegateInPlace(NekMatrix<double, StandardMatrixTag>& v);
854 
855 // template<typename DataType>
856 // template<typename T, typename MatrixType>
857 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::iterator_impl(typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::pointer d, typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::pointer e, bool isEnd) :
858 // m_data(d),
859 // m_end(e),
860 // m_curRow(std::numeric_limits<unsigned int>::max()),
861 // m_curColumn(std::numeric_limits<unsigned int>::max()),
862 // m_matrix(NULL),
863 // m_curIndex(std::numeric_limits<unsigned int>::max()),
864 // m_transpose('N')
865 // {
866 // if( isEnd )
867 // {
868 // m_data = m_end;
869 // }
870 // }
871 
872 // template<typename DataType>
873 // template<typename T, typename MatrixType>
874 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::iterator_impl(MatrixType* m, char transpose, bool isEnd) :
875 // m_data(NULL),
876 // m_end(NULL),
877 // m_curRow(0),
878 // m_curColumn(0),
879 // m_matrix(m),
880 // m_curIndex(0),
881 // m_transpose(transpose)
882 // {
883 // if( isEnd )
884 // {
885 // m_curRow = std::numeric_limits<unsigned int>::max();
886 // m_curColumn = std::numeric_limits<unsigned int>::max();
887 // m_curIndex = std::numeric_limits<unsigned int>::max();
888 // }
889 // }
890 
891 // template<typename DataType>
892 // template<typename T, typename MatrixType>
893 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::iterator_impl(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs) :
894 // m_data(rhs.m_data),
895 // m_end(rhs.m_end),
896 // m_curRow(rhs.m_curRow),
897 // m_curColumn(rhs.m_curColumn),
898 // m_matrix(rhs.m_matrix),
899 // m_curIndex(rhs.m_curIndex),
900 // m_transpose(rhs.m_transpose)
901 // {
902 // }
903 
904 // template<typename DataType>
905 // template<typename T, typename MatrixType>
906 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>&
907 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator=(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
908 // {
909 // m_data = rhs.m_data;
910 // m_end = rhs.m_end;
911 // m_curRow = rhs.m_curRow;
912 // m_curColumn = rhs.m_curColumn;
913 // m_matrix = rhs.m_matrix;
914 // m_curIndex = rhs.m_curIndex;
915 // m_transpose = rhs.m_transpose;
916 // return *this;
917 // }
918 
919 // template<typename DataType>
920 // template<typename T, typename MatrixType>
921 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::reference NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator*()
922 // {
923 // if( m_data )
924 // {
925 // ASSERTL1(m_data < m_end, "Attempt to dereference matrix iterator after its end.");
926 // return *m_data;
927 // }
928 // else
929 // {
930 // return m_matrix->GetPtr()[m_curIndex];
931 // }
932 // }
933 
934 // template<typename DataType>
935 // template<typename T, typename MatrixType>
936 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::const_reference NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator*() const
937 // {
938 // if( m_data )
939 // {
940 // ASSERTL1(m_data < m_end, "Attempt to dereference matrix iterator after its end.");
941 // return *m_data;
942 // }
943 // else
944 // {
945 // return m_matrix->GetPtr()[m_curIndex];
946 // }
947 // }
948 
949 // template<typename DataType>
950 // template<typename T, typename MatrixType>
951 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator++()
952 // {
953 // if( m_data )
954 // {
955 // ++m_data;
956 // }
957 // else
958 // {
959 // boost::tie(m_curRow, m_curColumn) =
960 // m_matrix->Advance(m_curRow, m_curColumn, m_transpose);
961 // if( m_curRow == std::numeric_limits<unsigned int>::max() )
962 // {
963 // m_curIndex = m_curRow;
964 // }
965 // else
966 // {
967 // m_curIndex = m_matrix->CalculateIndex(m_curRow, m_curColumn, m_transpose);
968 // }
969 // }
970 // return *this;
971 // }
972 
973 // template<typename DataType>
974 // template<typename T, typename MatrixType>
975 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType> NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator++(int)
976 // {
977 // iterator_impl<T, MatrixType> result = *this;
978 // ++(*this);
979 // return result;
980 // }
981 
982 // template<typename DataType>
983 // template<typename T, typename MatrixType>
984 // bool NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator==(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
985 // {
986 // return m_data == rhs.m_data &&
987 // m_end == rhs.m_end &&
988 // m_curRow == rhs.m_curRow &&
989 // m_curColumn == rhs.m_curColumn &&
990 // m_matrix == rhs.m_matrix &&
991 // m_curIndex == rhs.m_curIndex &&
992 // m_transpose == rhs.m_transpose;
993 // }
994 
995 // template<typename DataType>
996 // template<typename T, typename MatrixType>
997 // bool NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator!=(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
998 // {
999 // return !(*this == rhs);
1000 // }
1001 
1002 // template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>::iterator_impl<const NekDouble, const NekMatrix<NekDouble, StandardMatrixTag> >;
1003 // template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>::iterator_impl<NekDouble, NekMatrix<NekDouble, StandardMatrixTag> >;
1004 
1005 }
1006 
1007 
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2