Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 
33 
34 
35 namespace Nektar
36 {
37 
38 
39  template<typename DataType>
41  Matrix<DataType>(0, 0),
42  m_data(),
43  m_wrapperType(eCopy),
44  m_numberOfSuperDiagonals(std::numeric_limits<unsigned int>::max()),
45  m_numberOfSubDiagonals(std::numeric_limits<unsigned int>::max()),
46  m_tempSpace()
47  {
48  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
49  }
50 
51  template<typename DataType>
52  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, MatrixStorage policy,
53  unsigned int subDiagonals,
54  unsigned int superDiagonals) :
55  Matrix<DataType>(rows, columns, policy),
56  m_data(),
57  m_wrapperType(eCopy),
58  m_numberOfSuperDiagonals(superDiagonals),
59  m_numberOfSubDiagonals(subDiagonals),
60  m_tempSpace()
61  {
62  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
63  }
64 
65  template<typename DataType>
66  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns,
67  MatrixStorage policy,
68  unsigned int subDiagonals,
69  unsigned int superDiagonals,
70  unsigned int capacity) :
71  Matrix<DataType>(rows, columns, policy),
72  m_data(),
73  m_wrapperType(eCopy),
74  m_numberOfSuperDiagonals(superDiagonals),
75  m_numberOfSubDiagonals(subDiagonals),
76  m_tempSpace()
77  {
78  unsigned int requiredStorage = this->GetRequiredStorageSize();
79  unsigned int actualSize = std::max(requiredStorage, capacity);
80  m_data = Array<OneD, DataType>(actualSize);
81  }
82 
83  template<typename DataType>
84  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, typename boost::call_traits<DataType>::const_reference initValue,
85  MatrixStorage policy,
86  unsigned int subDiagonals,
87  unsigned int superDiagonals) :
88  Matrix<DataType>(rows, columns, policy),
89  m_data(),
90  m_wrapperType(eCopy),
91  m_numberOfSuperDiagonals(superDiagonals),
92  m_numberOfSubDiagonals(subDiagonals),
93  m_tempSpace()
94  {
95  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
96  std::fill(m_data.begin(), m_data.end(), initValue);
97  }
98 
99  template<typename DataType>
100  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const DataType* data,
101  MatrixStorage policy,
102  unsigned int subDiagonals,
103  unsigned int superDiagonals) :
104  Matrix<DataType>(rows, columns, policy),
105  m_data(),
106  m_wrapperType(eCopy),
107  m_numberOfSuperDiagonals(superDiagonals),
108  m_numberOfSubDiagonals(subDiagonals),
109  m_tempSpace()
110  {
111  unsigned int size = GetRequiredStorageSize();
112  m_data = Array<OneD, DataType>(size);
113  std::copy(data, data + size, m_data.begin());
114  }
115 
116  template<typename DataType>
117  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const Array<OneD, const DataType>& d,
118  MatrixStorage policy,
119  unsigned int subDiagonals,
120  unsigned int superDiagonals) :
121  Matrix<DataType>(rows, columns, policy),
122  m_data(),
123  m_wrapperType(eCopy),
124  m_numberOfSuperDiagonals(superDiagonals),
125  m_numberOfSubDiagonals(subDiagonals),
126  m_tempSpace()
127  {
128  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
129  CopyArrayN(d, m_data, m_data.num_elements());
130  }
131 
132  template<typename DataType>
133  NekMatrix<DataType, StandardMatrixTag>::NekMatrix(unsigned int rows, unsigned int columns, const Array<OneD, DataType>& d, PointerWrapper wrapperType,
134  MatrixStorage policy,
135  unsigned int subDiagonals,
136  unsigned int superDiagonals) :
137  Matrix<DataType>(rows, columns, policy),
138  m_data(),
139  m_wrapperType(wrapperType),
140  m_numberOfSuperDiagonals(superDiagonals),
141  m_numberOfSubDiagonals(subDiagonals),
142  m_tempSpace()
143  {
144  if( wrapperType == eWrapper )
145  {
147  }
148  else
149  {
150  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
151  CopyArrayN(d, m_data, m_data.num_elements());
152  }
153  }
154 
155  template<typename DataType>
157  Matrix<DataType>(rhs),
158  m_data(),
159  m_wrapperType(rhs.m_wrapperType),
160  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
161  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals),
162  m_tempSpace()
163  {
164  PerformCopyConstruction(rhs);
165  }
166 
167  template<typename DataType>
169  {
170  if( this == &rhs )
171  {
172  return *this;
173  }
174 
176  m_numberOfSubDiagonals = rhs.m_numberOfSubDiagonals;
177  m_numberOfSuperDiagonals = rhs.m_numberOfSuperDiagonals;
178 
179  ResizeDataArrayIfNeeded();
180 
181  unsigned int requiredStorageSize = GetRequiredStorageSize();
182  std::copy(rhs.m_data.data(), rhs.m_data.data() + requiredStorageSize, m_data.data());
183 
184  return *this;
185  }
186 
187 
188  template<typename DataType>
190  {
191  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
192  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
193  std::string(" rows"));
194  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
195  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
196  std::string(" columns"));
197 
198  return this->GetValue(row, column, this->GetTransposeFlag());
199  }
200 
201  template<typename DataType>
202  typename NekMatrix<DataType, StandardMatrixTag>::ConstGetValueType NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column, char transpose) const
203  {
204  return this->GetValue(row, column, transpose);
205  }
206 
207  template<typename DataType>
209  {
210  return Matrix<DataType>::GetRequiredStorageSize(this->GetStorageType(),
211  this->GetRows(), this->GetColumns(),
212  this->GetNumberOfSubDiagonals(), this->GetNumberOfSuperDiagonals());
213  }
214 
215  template<typename DataType>
216  unsigned int NekMatrix<DataType, StandardMatrixTag>::CalculateIndex(unsigned int row, unsigned int col, const char transpose) const
217  {
218  unsigned int numRows = this->GetSize()[0];
219  unsigned int numColumns = this->GetSize()[1];
220  return Matrix<DataType>::CalculateIndex(this->GetStorageType(),
221  row, col, numRows, numColumns, transpose,
222  m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
223  }
224 
225  template<typename DataType>
226  typename boost::call_traits<DataType>::const_reference NekMatrix<DataType, StandardMatrixTag>::GetValue(unsigned int row, unsigned int column) const
227  {
228  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
229  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
230  std::string(" rows"));
231  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
232  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
233  std::string(" columns"));
234 
235  return GetValue(row, column, this->GetTransposeFlag());
236  }
237 
238  template<typename DataType>
239  typename NekMatrix<DataType, StandardMatrixTag>::ConstGetValueType NekMatrix<DataType, StandardMatrixTag>::GetValue(unsigned int row, unsigned int column, char transpose) const
240  {
241  static DataType defaultReturnValue;
242  unsigned int index = CalculateIndex(row, column, transpose);
243  if( index != std::numeric_limits<unsigned int>::max() )
244  {
245  return m_data[index];
246  }
247  else
248  {
249  return defaultReturnValue;
250  }
251  }
252 
253  template<typename DataType>
255  {
256  return m_data;
257  }
258 
259  template<typename DataType>
261  {
262  return DataType(1);
263  }
264 
265  template<typename DataType>
267  {
268  return m_data.data();
269  }
270 
271  template<typename DataType>
273  {
274  return begin(this->GetTransposeFlag());
275  }
276 
277  template<typename DataType>
279  {
280  if( transpose == 'N' )
281  {
282  return const_iterator(m_data.data(), m_data.data() + m_data.num_elements());
283  }
284  else
285  {
286  return const_iterator(this, transpose);
287  }
288  }
289 
290  template<typename DataType>
292  {
293  return end(this->GetTransposeFlag());
294  }
295 
296  template<typename DataType>
298  {
299  if( transpose == 'N' )
300  {
301  return const_iterator(m_data.data(), m_data.data() + m_data.num_elements(), true);
302  }
303  else
304  {
305  return const_iterator(this, transpose, true);
306  }
307  }
308 
309  template<typename DataType>
311  {
312  return m_data.num_elements();
313  }
314 
315  template<typename DataType>
317  {
318  if( m_numberOfSubDiagonals != std::numeric_limits<unsigned int>::max() )
319  {
320  return m_numberOfSubDiagonals;
321  }
322  else if( this->GetRows() > 0 )
323  {
324  return this->GetRows()-1;
325  }
326  else
327  {
328  return 0;
329  }
330  }
331 
332  template<typename DataType>
334  {
335  if( m_numberOfSuperDiagonals != std::numeric_limits<unsigned int>::max() )
336  {
337  return m_numberOfSuperDiagonals;
338  }
339  else if( this->GetRows() > 0 )
340  {
341  return this->GetRows()-1;
342  }
343  else
344  {
345  return 0;
346  }
347  }
348 
349  template<typename DataType>
351  {
352  return GetNumberOfSubDiagonals() + GetNumberOfSuperDiagonals() + 1;
353  }
354 
355  template<typename DataType>
357  {
358  if( this->GetRows() != rhs.GetRows() ||
359  this->GetColumns() != rhs.GetColumns() )
360  {
361  return false;
362  }
363 
364  if( this->GetTransposeFlag() == rhs.GetTransposeFlag() )
365  {
366  return std::equal(begin(), end(), rhs.begin());
367  }
368  else
369  {
370  for(unsigned int i = 0; i < this->GetRows(); ++i)
371  {
372  for(unsigned int j = 0; j < this->GetColumns(); ++j)
373  {
374  if( (*this)(i,j) != rhs(i,j) )
375  {
376  return false;
377  }
378  }
379  }
380  }
381 
382  return true;
383  }
384 
385  template<typename DataType>
387 
388  template<typename DataType>
389  boost::tuples::tuple<unsigned int, unsigned int>
390  NekMatrix<DataType, StandardMatrixTag>::Advance(unsigned int curRow, unsigned int curColumn) const
391  {
392  return Advance(curRow, curColumn, this->GetTransposeFlag());
393  }
394 
395  template<typename DataType>
396  boost::tuples::tuple<unsigned int, unsigned int>
397  NekMatrix<DataType, StandardMatrixTag>::Advance(unsigned int curRow, unsigned int curColumn, char transpose) const
398  {
399  unsigned int numRows = this->GetTransposedRows(transpose);
400  unsigned int numColumns = this->GetTransposedColumns(transpose);
401 
402  switch(this->GetStorageType())
403  {
404  case eFULL:
406  numRows, numColumns, curRow, curColumn);
407  break;
408  case eDIAGONAL:
410  numRows, numColumns, curRow, curColumn);
411  break;
412  case eUPPER_TRIANGULAR:
414  numRows, numColumns, curRow, curColumn);
415  break;
416 
417  case eLOWER_TRIANGULAR:
419  numRows, numColumns, curRow, curColumn);
420  break;
421 
422  case eSYMMETRIC:
425  numRows, numColumns, curRow, curColumn);
426  break;
427  case eBANDED:
429  numRows, numColumns, curRow, curColumn);
430  break;
431  case eSYMMETRIC_BANDED:
433  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
434  break;
436  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
437  break;
439  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
440  break;
441 
442  default:
443  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
444  }
445  return boost::tuples::tuple<unsigned int, unsigned int>(curRow, curColumn);
446  }
447 
448  template<typename DataType>
450  {
452  }
453 
454  template<typename DataType>
455  boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> > NekMatrix<DataType, StandardMatrixTag>::CreateWrapper(const boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >& rhs)
456  {
457  return boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >(new NekMatrix<DataType, StandardMatrixTag>(*rhs, eWrapper));
458  }
459 
460  template<typename DataType>
462  BaseType(rhs),
463  m_data(),
464  m_wrapperType(wrapperType),
465  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
466  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals)
467  {
468  PerformCopyConstruction(rhs);
469  }
470 
471 
472  template<typename DataType>
474 
475  template<typename DataType>
477  {
478  if( m_wrapperType == eCopy )
479  {
480  unsigned int requiredStorageSize = GetRequiredStorageSize();
481  if( m_data.num_elements() > requiredStorageSize )
482  {
483  Array<OneD, DataType> newArray(requiredStorageSize);
484  CopyArrayN(m_data, newArray, requiredStorageSize);
485  m_data = newArray;
486  }
487  }
488  else if( m_wrapperType == eWrapper )
489  {
490  ASSERTL0(true, "Can't call RemoveExcessCapacity on a wrapped matrix.");
491  }
492  }
493 
494  template<typename DataType>
496  {
497  if( m_wrapperType == eCopy )
498  {
499  if( m_data.num_elements() < requiredStorageSize )
500  {
501  Array<OneD, DataType> newData(requiredStorageSize);
502  std::copy(m_data.data(), m_data.data() + m_data.num_elements(), newData.data());
503  m_data = newData;
504  }
505  }
506  else if( m_wrapperType == eWrapper )
507  {
508  // If the current matrix is wrapped, then just copy over the top,
509  // but the sizes of the two matrices must be the same.
510  ASSERTL0(m_data.num_elements() >= requiredStorageSize, "Wrapped NekMatrices must have the same dimension in operator=");
511  }
512  }
513 
514  template<typename DataType>
516  {
517  unsigned int requiredStorageSize = GetRequiredStorageSize();
518  ResizeDataArrayIfNeeded(requiredStorageSize);
519  }
520 
521  template<typename DataType>
523  {
524  if( m_wrapperType == eWrapper )
525  {
526  m_data = rhs.m_data;
527  }
528  else
529  {
530  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
531  CopyArrayN(rhs.m_data, m_data, m_data.num_elements());
532  }
533  }
534 
535  template<typename DataType>
536  typename boost::call_traits<DataType>::value_type NekMatrix<DataType, StandardMatrixTag>::v_GetValue(unsigned int row, unsigned int column) const
537  {
539  }
540 
541  template<typename DataType>
543  {
545  }
546 
547  template<typename DataType>
548  void NekMatrix<DataType, StandardMatrixTag>::v_SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d)
549  {
551  }
552 
553 
554 
555  template<typename DataType>
556  void NekMatrix<DataType, StandardMatrixTag>::SetSize(unsigned int rows, unsigned int cols)
557  {
558  this->Resize(rows, cols);
559 
560  // Some places in Nektar++ access the matrix data array and
561  // use num_elements() to see how big it is. When using
562  // expression templates, the data array's capacity is often larger
563  // than the actual number of elements, so this statement is
564  // required to report the correct number of elements.
565  this->GetData().ChangeSize(this->GetRequiredStorageSize());
566  ASSERTL0(this->GetRequiredStorageSize() <= this->GetData().num_elements(), "Can't resize matrices if there is not enough capacity.");
567  }
568 
569 
570  template<typename DataType>
572  {
573  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
574  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
575  std::string(" rows"));
576  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
577  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
578  std::string(" columns"));
579 
580  return (*this)(row, column, this->GetTransposeFlag());
581  }
582 
583  template<typename DataType>
584  typename NekMatrix<DataType, StandardMatrixTag>::Proxy NekMatrix<DataType, StandardMatrixTag>::operator()(unsigned int row, unsigned int column, char transpose)
585  {
586  unsigned int index = this->CalculateIndex(row, column, transpose);
587  if( index != std::numeric_limits<unsigned int>::max() )
588  {
589  return Proxy(this->GetData()[index]);
590  }
591  else
592  {
593  return Proxy();
594  }
595 
596  }
597 
598  template<typename DataType>
599  void NekMatrix<DataType, StandardMatrixTag>::SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d)
600  {
601  ASSERTL2(row < this->GetRows(), std::string("Row ") + boost::lexical_cast<std::string>(row) +
602  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetRows()) +
603  std::string(" rows"));
604  ASSERTL2(column < this->GetColumns(), std::string("Column ") + boost::lexical_cast<std::string>(column) +
605  std::string(" requested in a matrix with a maximum of ") + boost::lexical_cast<std::string>(this->GetColumns()) +
606  std::string(" columns"));
607  SetValue(row, column, d, this->GetTransposeFlag());
608  }
609 
610  template<typename DataType>
611  void NekMatrix<DataType, StandardMatrixTag>::SetValue(unsigned int row, unsigned int column, typename boost::call_traits<DataType>::const_reference d, char transpose)
612  {
613  unsigned int index = this->CalculateIndex(row, column, transpose);
614  if( index != std::numeric_limits<unsigned int>::max() )
615  {
616  this->GetData()[index] = d;
617  }
618  else
619  {
620  NEKERROR(ErrorUtil::efatal, "Can't assign values into zeroed elements of a special array.");
621  }
622  }
623 
624  template<typename DataType>
626  {
627  return this->GetData();
628  }
629 
630  template<typename DataType>
632  {
633  return this->GetData().data();
634  }
635 
636  template<typename DataType>
638  {
639  return begin(this->GetTransposeFlag());
640  }
641 
642  template<typename DataType>
644  {
645  if( transpose == 'N' )
646  {
647  return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements());
648  }
649  else
650  {
651  return iterator(this, transpose);
652  }
653  }
654 
655  template<typename DataType>
657  {
658  return end(this->GetTransposeFlag());
659  }
660 
661  template<typename DataType>
663  {
664  if( transpose == 'N' )
665  {
666  return iterator(this->GetData().data(), this->GetData().data() + this->GetData().num_elements(), true);
667  }
668  else
669  {
670  return iterator(this, transpose, true);
671  }
672  }
673 
674  template<typename DataType>
676  {
677  NekDouble returnval;
678  int nvals = this->GetColumns();
679  Array<OneD, NekDouble> EigValReal(nvals);
680  Array<OneD, NekDouble> EigValImag(nvals);
681 
682  EigenSolve(EigValReal,EigValImag);
683 
684  Vmath::Vmul(nvals,EigValReal,1,EigValReal,1,EigValReal,1);
685  Vmath::Vmul(nvals,EigValImag,1,EigValImag,1,EigValImag,1);
686  Vmath::Vadd(nvals,EigValReal,1,EigValImag,1,EigValReal,1);
687 
688  returnval = sqrt(Vmath::Vmax(nvals,EigValReal,1)/Vmath::Vmin(nvals,EigValReal,1));
689 
690  return returnval;
691  }
692 
693  template<typename DataType>
695  Array<OneD, NekDouble> &EigValImag,
696  Array<OneD, NekDouble> &EigVecs)
697  {
698  ASSERTL0(this->GetRows()==this->GetColumns(), "Only square matrices can be called");
699 
700  switch(this->GetType())
701  {
702  case eFULL:
703  FullMatrixFuncs::EigenSolve(this->GetRows(),
704  this->GetData(), EigValReal,
705  EigValImag, EigVecs);
706  break;
707  case eDIAGONAL:
708  Vmath::Vcopy(this->GetRows(),&(this->GetData())[0],1, &EigValReal[0],1);
709  Vmath::Zero(this->GetRows(),&EigValImag[0],1);
710  break;
711  case eUPPER_TRIANGULAR:
712  case eLOWER_TRIANGULAR:
713  case eSYMMETRIC:
714  case eBANDED:
715  case eSYMMETRIC_BANDED:
718  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
719  break;
720 
721  default:
722  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
723  }
724  }
725 
726  template<typename DataType>
728  {
729  ASSERTL0(this->GetRows()==this->GetColumns(), "Only square matrices can be inverted.");
730  ASSERTL0(this->GetTransposeFlag()=='N', "Only untransposed matrices may be inverted.");
731 
732  switch(this->GetType())
733  {
734  case eFULL:
735  FullMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
736  this->GetData(), this->GetTransposeFlag());
737  break;
738  case eDIAGONAL:
739  DiagonalMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
740  this->GetData());
741  break;
742  case eSYMMETRIC:
743  SymmetricMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
744  this->GetData());
745  break;
746  case eUPPER_TRIANGULAR:
747  case eLOWER_TRIANGULAR:
748  case eBANDED:
749  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
750  break;
751  case eSYMMETRIC_BANDED:
752  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
753  break;
755  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
756  break;
758  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
759  break;
760 
761  default:
762  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
763  }
764  }
765 
766  template<typename DataType>
768  {
769  if( m_tempSpace.capacity() == 0 )
770  {
771  m_tempSpace = Array<OneD, DataType>(this->GetData().capacity());
772  }
773  return m_tempSpace;
774  }
775 
776  template<typename DataType>
778  {
779  std::swap(m_tempSpace, this->GetData());
780  }
781 
782  #ifndef NEKTAR_USE_EXPRESSION_TEMPLATES
783  template<typename DataType>
785  {
787  NegateInPlace(result);
788  return result;
789  }
790  #endif
791 
792  template<typename DataType>
794  {
795  for(unsigned int i = 0; i < this->GetPtr().num_elements(); ++i)
796  {
797  this->GetPtr()[i] *= s;
798  }
799  return *this;
800  }
801 
802 
803  template<typename DataType>
806  {
808  rhs.GetPtr(), eWrapper, rhs.GetType(), rhs.GetNumberOfSubDiagonals(),
810  result.Transpose();
811  return result;
812  }
813 
814  template LIB_UTILITIES_EXPORT NekMatrix<NekDouble, StandardMatrixTag> Transpose(NekMatrix<NekDouble, StandardMatrixTag>& rhs);
815 
816  template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>;
817 
818  template<typename DataType>
820  {
821  for(unsigned int i = 0; i < m.GetRows(); ++i)
822  {
823  for(unsigned int j = 0; j < m.GetColumns(); ++j)
824  {
825  m(i,j) *= -1.0;
826  }
827  }
828  }
829 
830  template LIB_UTILITIES_EXPORT void NegateInPlace(NekMatrix<double, StandardMatrixTag>& v);
831 
832 // template<typename DataType>
833 // template<typename T, typename MatrixType>
834 // 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) :
835 // m_data(d),
836 // m_end(e),
837 // m_curRow(std::numeric_limits<unsigned int>::max()),
838 // m_curColumn(std::numeric_limits<unsigned int>::max()),
839 // m_matrix(NULL),
840 // m_curIndex(std::numeric_limits<unsigned int>::max()),
841 // m_transpose('N')
842 // {
843 // if( isEnd )
844 // {
845 // m_data = m_end;
846 // }
847 // }
848 
849 // template<typename DataType>
850 // template<typename T, typename MatrixType>
851 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::iterator_impl(MatrixType* m, char transpose, bool isEnd) :
852 // m_data(NULL),
853 // m_end(NULL),
854 // m_curRow(0),
855 // m_curColumn(0),
856 // m_matrix(m),
857 // m_curIndex(0),
858 // m_transpose(transpose)
859 // {
860 // if( isEnd )
861 // {
862 // m_curRow = std::numeric_limits<unsigned int>::max();
863 // m_curColumn = std::numeric_limits<unsigned int>::max();
864 // m_curIndex = std::numeric_limits<unsigned int>::max();
865 // }
866 // }
867 
868 // template<typename DataType>
869 // template<typename T, typename MatrixType>
870 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::iterator_impl(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs) :
871 // m_data(rhs.m_data),
872 // m_end(rhs.m_end),
873 // m_curRow(rhs.m_curRow),
874 // m_curColumn(rhs.m_curColumn),
875 // m_matrix(rhs.m_matrix),
876 // m_curIndex(rhs.m_curIndex),
877 // m_transpose(rhs.m_transpose)
878 // {
879 // }
880 
881 // template<typename DataType>
882 // template<typename T, typename MatrixType>
883 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>&
884 // NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator=(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
885 // {
886 // m_data = rhs.m_data;
887 // m_end = rhs.m_end;
888 // m_curRow = rhs.m_curRow;
889 // m_curColumn = rhs.m_curColumn;
890 // m_matrix = rhs.m_matrix;
891 // m_curIndex = rhs.m_curIndex;
892 // m_transpose = rhs.m_transpose;
893 // return *this;
894 // }
895 
896 // template<typename DataType>
897 // template<typename T, typename MatrixType>
898 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::reference NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator*()
899 // {
900 // if( m_data )
901 // {
902 // ASSERTL1(m_data < m_end, "Attempt to dereference matrix iterator after its end.");
903 // return *m_data;
904 // }
905 // else
906 // {
907 // return m_matrix->GetPtr()[m_curIndex];
908 // }
909 // }
910 
911 // template<typename DataType>
912 // template<typename T, typename MatrixType>
913 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>::const_reference NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator*() const
914 // {
915 // if( m_data )
916 // {
917 // ASSERTL1(m_data < m_end, "Attempt to dereference matrix iterator after its end.");
918 // return *m_data;
919 // }
920 // else
921 // {
922 // return m_matrix->GetPtr()[m_curIndex];
923 // }
924 // }
925 
926 // template<typename DataType>
927 // template<typename T, typename MatrixType>
928 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator++()
929 // {
930 // if( m_data )
931 // {
932 // ++m_data;
933 // }
934 // else
935 // {
936 // boost::tie(m_curRow, m_curColumn) =
937 // m_matrix->Advance(m_curRow, m_curColumn, m_transpose);
938 // if( m_curRow == std::numeric_limits<unsigned int>::max() )
939 // {
940 // m_curIndex = m_curRow;
941 // }
942 // else
943 // {
944 // m_curIndex = m_matrix->CalculateIndex(m_curRow, m_curColumn, m_transpose);
945 // }
946 // }
947 // return *this;
948 // }
949 
950 // template<typename DataType>
951 // template<typename T, typename MatrixType>
952 // typename NekMatrix<DataType>::template iterator_impl<T, MatrixType> NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator++(int)
953 // {
954 // iterator_impl<T, MatrixType> result = *this;
955 // ++(*this);
956 // return result;
957 // }
958 
959 // template<typename DataType>
960 // template<typename T, typename MatrixType>
961 // bool NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator==(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
962 // {
963 // return m_data == rhs.m_data &&
964 // m_end == rhs.m_end &&
965 // m_curRow == rhs.m_curRow &&
966 // m_curColumn == rhs.m_curColumn &&
967 // m_matrix == rhs.m_matrix &&
968 // m_curIndex == rhs.m_curIndex &&
969 // m_transpose == rhs.m_transpose;
970 // }
971 
972 // template<typename DataType>
973 // template<typename T, typename MatrixType>
974 // bool NekMatrix<DataType>::iterator_impl<T, MatrixType>::operator!=(const typename NekMatrix<DataType>::template iterator_impl<T, MatrixType>& rhs)
975 // {
976 // return !(*this == rhs);
977 // }
978 
979 // template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>::iterator_impl<const NekDouble, const NekMatrix<NekDouble, StandardMatrixTag> >;
980 // template LIB_UTILITIES_EXPORT class NekMatrix<NekDouble, StandardMatrixTag>::iterator_impl<NekDouble, NekMatrix<NekDouble, StandardMatrixTag> >;
981 
982 }
983 
984 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:191
iterator_impl< const DataType, const ThisType > const_iterator
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
Definition: MatrixFuncs.h:86
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:871
STL namespace.
char GetTransposeFlag() const
Definition: MatrixBase.cpp:114
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
unsigned int GetColumns() const
Definition: MatrixBase.cpp:78
size_type num_elements() const
Returns the array's size.
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
Matrix< DataType > & operator=(const Matrix< DataType > &rhs)
Definition: MatrixBase.cpp:306
static unsigned int CalculateIndex(MatrixStorage type, unsigned int row, unsigned int col, unsigned int numRows, unsigned int numColumns, const char transpose= 'N', unsigned int numSubDiags=0, unsigned int numSuperDiags=0)
Definition: MatrixBase.cpp:120
#define LIB_UTILITIES_EXPORT
void CopyArrayN(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest, unsigned int n)
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
int GetSize(const ConstArray< OneD, NekDouble > &x)
Definition: StdExpUtil.cpp:111
double NekDouble
NekMatrix< NekDouble > Invert(const NekMatrix< NekDouble > &matA)
Definition: StdExpUtil.cpp:79
NekPoint< DataType > operator-(const NekPoint< DataType > &lhs, const NekPoint< DataType > &rhs)
Definition: NekPoint.hpp:407
const Array< OneD, const DataType > & GetPtr() const
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
unsigned int GetRows() const
Definition: MatrixBase.cpp:59
MatrixStorage GetType() const
Definition: MatrixBase.hpp:67
static void EigenSolve(unsigned int n, const Array< OneD, const double > &A, Array< OneD, NekDouble > &EigValReal, Array< OneD, NekDouble > &EigValImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
Definition: MatrixFuncs.h:133
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn, char transpose= 'N')
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
static unsigned int GetRequiredStorageSize(MatrixStorage type, unsigned int rows, unsigned int columns, unsigned int subDiags=0, unsigned int superDiags=0)
Definition: MatrixBase.cpp:176
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:206
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:262
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
1D Array of constant elements with garbage collection and bounds checking.
Definition: SharedArray.hpp:60
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
Definition: VtkToFld.cpp:54
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
void NegateInPlace(NekVector< DataType > &v)
Definition: NekVector.cpp:962
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
Definition: MatrixFuncs.cpp:91