Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
12 // copy of this software and associated documentation files (the "Software"),
13 // to deal in the Software without restriction, including without limitation
14 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 // and/or sell copies of the Software, and to permit persons to whom the
16 // Software is furnished to do so, subject to the following conditions:
17 //
18 // The above copyright notice and this permission notice shall be included
19 // in all copies or substantial portions of the Software.
20 //
21 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 // DEALINGS IN THE SOFTWARE.
28 //
29 ///////////////////////////////////////////////////////////////////////////////
30 
32 
33 namespace Nektar
34 {
35 
36 template <typename DataType>
38  : Matrix<DataType>(0, 0, eFULL), m_data(), m_wrapperType(eCopy),
39  m_numberOfSuperDiagonals(std::numeric_limits<unsigned int>::max()),
40  m_numberOfSubDiagonals(std::numeric_limits<unsigned int>::max()),
41  m_tempSpace()
42 {
43  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
44 }
45 
46 template <typename DataType>
48  unsigned int columns,
49  MatrixStorage policy,
50  unsigned int subDiagonals,
51  unsigned int superDiagonals)
52  : Matrix<DataType>(rows, columns, policy), m_data(), m_wrapperType(eCopy),
53  m_numberOfSuperDiagonals(superDiagonals),
54  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
55 {
56  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
57 }
58 
59 template <typename DataType>
61  unsigned int columns,
62  MatrixStorage policy,
63  unsigned int subDiagonals,
64  unsigned int superDiagonals,
65  unsigned int capacity)
66  : Matrix<DataType>(rows, columns, policy), m_data(), m_wrapperType(eCopy),
67  m_numberOfSuperDiagonals(superDiagonals),
68  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
69 {
70  unsigned int requiredStorage = this->GetRequiredStorageSize();
71  unsigned int actualSize = std::max(requiredStorage, capacity);
72  m_data = Array<OneD, DataType>(actualSize);
73 }
74 
75 template <typename DataType>
77  unsigned int rows, unsigned int columns,
78  typename boost::call_traits<DataType>::const_reference initValue,
79  MatrixStorage policy, unsigned int subDiagonals,
80  unsigned int superDiagonals)
81  : Matrix<DataType>(rows, columns, policy), m_data(), m_wrapperType(eCopy),
82  m_numberOfSuperDiagonals(superDiagonals),
83  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
84 {
85  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
86  std::fill(m_data.begin(), m_data.end(), initValue);
87 }
88 
89 template <typename DataType>
91  unsigned int columns,
92  const DataType *data,
93  MatrixStorage policy,
94  unsigned int subDiagonals,
95  unsigned int superDiagonals)
96  : Matrix<DataType>(rows, columns, policy), m_data(), m_wrapperType(eCopy),
97  m_numberOfSuperDiagonals(superDiagonals),
98  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
99 {
100  unsigned int size = GetRequiredStorageSize();
101  m_data = Array<OneD, DataType>(size);
102  std::copy(data, data + size, m_data.begin());
103 }
104 
105 template <typename DataType>
107  unsigned int rows, unsigned int columns,
109  unsigned int subDiagonals, unsigned int superDiagonals)
110  : Matrix<DataType>(rows, columns, policy), m_data(), m_wrapperType(eCopy),
111  m_numberOfSuperDiagonals(superDiagonals),
112  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
113 {
114  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
115  CopyArrayN(d, m_data, m_data.size());
116 }
117 
118 template <typename DataType>
120  unsigned int rows, unsigned int columns, const Array<OneD, DataType> &d,
121  PointerWrapper wrapperType, MatrixStorage policy, unsigned int subDiagonals,
122  unsigned int superDiagonals)
123  : Matrix<DataType>(rows, columns, policy), m_data(),
124  m_wrapperType(wrapperType), m_numberOfSuperDiagonals(superDiagonals),
125  m_numberOfSubDiagonals(subDiagonals), m_tempSpace()
126 {
127  if (wrapperType == eWrapper)
128  {
130  }
131  else
132  {
133  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
134  CopyArrayN(d, m_data, m_data.size());
135  }
136 }
137 
138 template <typename DataType>
141  : Matrix<DataType>(rhs), m_data(), m_wrapperType(rhs.m_wrapperType),
142  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
143  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals), m_tempSpace()
144 {
145  PerformCopyConstruction(rhs);
146 }
147 
148 template <typename DataType>
152 {
153  if (this == &rhs)
154  {
155  return *this;
156  }
157 
159  m_numberOfSubDiagonals = rhs.m_numberOfSubDiagonals;
160  m_numberOfSuperDiagonals = rhs.m_numberOfSuperDiagonals;
161 
162  ResizeDataArrayIfNeeded();
163 
164  unsigned int requiredStorageSize = GetRequiredStorageSize();
165  std::copy(rhs.m_data.data(), rhs.m_data.data() + requiredStorageSize,
166  m_data.data());
167 
168  return *this;
169 }
170 
171 /// Fill matrix with scalar
172 template <typename DataType>
175 {
176  unsigned int requiredStorageSize = GetRequiredStorageSize();
177 
178  DataType *lhs_array = m_data.data();
179 
180  Vmath::Fill(requiredStorageSize, rhs, lhs_array, 1);
181 
182  return *this;
183 }
184 
185 template <typename DataType>
187  DataType, StandardMatrixTag>::operator()(unsigned int row,
188  unsigned int column) const
189 {
190  ASSERTL2(row < this->GetRows(),
191  std::string("Row ") + std::to_string(row) +
192  std::string(" requested in a matrix with a maximum of ") +
193  std::to_string(this->GetRows()) + std::string(" rows"));
194  ASSERTL2(column < this->GetColumns(),
195  std::string("Column ") + std::to_string(column) +
196  std::string(" requested in a matrix with a maximum of ") +
197  std::to_string(this->GetColumns()) + std::string(" columns"));
198 
199  return this->GetValue(row, column, this->GetTransposeFlag());
200 }
201 
202 template <typename DataType>
204  DataType, StandardMatrixTag>::operator()(unsigned int row,
205  unsigned int column,
206  char transpose) const
207 {
208  return this->GetValue(row, column, transpose);
209 }
210 
211 template <typename DataType>
213  const
214 {
216  this->GetStorageType(), this->GetRows(), this->GetColumns(),
217  this->GetNumberOfSubDiagonals(), this->GetNumberOfSuperDiagonals());
218 }
219 
220 template <typename DataType>
222  unsigned int row, unsigned int col, const char transpose) const
223 {
224  unsigned int numRows = this->GetSize()[0];
225  unsigned int numColumns = this->GetSize()[1];
227  this->GetStorageType(), row, col, numRows, numColumns, transpose,
228  m_numberOfSubDiagonals, m_numberOfSuperDiagonals);
229 }
230 
231 template <typename DataType>
232 typename boost::call_traits<DataType>::const_reference NekMatrix<
233  DataType, StandardMatrixTag>::GetValue(unsigned int row,
234  unsigned int column) const
235 {
236  ASSERTL2(row < this->GetRows(),
237  std::string("Row ") + std::to_string(row) +
238  std::string(" requested in a matrix with a maximum of ") +
239  std::to_string(this->GetRows()) + std::string(" rows"));
240  ASSERTL2(column < this->GetColumns(),
241  std::string("Column ") + std::to_string(column) +
242  std::string(" requested in a matrix with a maximum of ") +
243  std::to_string(this->GetColumns()) + std::string(" columns"));
244 
245  return GetValue(row, column, this->GetTransposeFlag());
246 }
247 
248 template <typename DataType>
250  DataType, StandardMatrixTag>::GetValue(unsigned int row,
251  unsigned int column,
252  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>
269 {
270  return m_data;
271 }
272 
273 template <typename DataType>
275 {
276  return DataType(1);
277 }
278 
279 template <typename DataType>
281 {
282  return m_data.data();
283 }
284 
285 template <typename DataType>
287  DataType, StandardMatrixTag>::begin() const
288 {
289  return begin(this->GetTransposeFlag());
290 }
291 
292 template <typename DataType>
294  DataType, StandardMatrixTag>::begin(char transpose) const
295 {
296  if (transpose == 'N')
297  {
298  return const_iterator(m_data.data(), m_data.data() + m_data.size());
299  }
300  else
301  {
302  return const_iterator(this, transpose);
303  }
304 }
305 
306 template <typename DataType>
308  DataType, StandardMatrixTag>::end() const
309 {
310  return end(this->GetTransposeFlag());
311 }
312 
313 template <typename DataType>
315  DataType, StandardMatrixTag>::end(char transpose) const
316 {
317  if (transpose == 'N')
318  {
319  return const_iterator(m_data.data(), m_data.data() + m_data.size(),
320  true);
321  }
322  else
323  {
324  return const_iterator(this, transpose, true);
325  }
326 }
327 
328 template <typename DataType>
330 {
331  return m_data.size();
332 }
333 
334 template <typename DataType>
336  const
337 {
338  if (m_numberOfSubDiagonals != std::numeric_limits<unsigned int>::max())
339  {
340  return m_numberOfSubDiagonals;
341  }
342  else if (this->GetRows() > 0)
343  {
344  return this->GetRows() - 1;
345  }
346  else
347  {
348  return 0;
349  }
350 }
351 
352 template <typename DataType>
354  const
355 {
356  if (m_numberOfSuperDiagonals != std::numeric_limits<unsigned int>::max())
357  {
358  return m_numberOfSuperDiagonals;
359  }
360  else if (this->GetRows() > 0)
361  {
362  return this->GetRows() - 1;
363  }
364  else
365  {
366  return 0;
367  }
368 }
369 
370 template <typename DataType>
372  const
373 {
374  return GetNumberOfSubDiagonals() + GetNumberOfSuperDiagonals() + 1;
375 }
376 
377 template <typename DataType>
380 {
381  if (this->GetRows() != rhs.GetRows() ||
382  this->GetColumns() != rhs.GetColumns())
383  {
384  return false;
385  }
386 
387  if (this->GetTransposeFlag() == rhs.GetTransposeFlag())
388  {
389  return std::equal(begin(), end(), rhs.begin());
390  }
391  else
392  {
393  for (unsigned int i = 0; i < this->GetRows(); ++i)
394  {
395  for (unsigned int j = 0; j < this->GetColumns(); ++j)
396  {
397  if ((*this)(i, j) != rhs(i, j))
398  {
399  return false;
400  }
401  }
402  }
403  }
404 
405  return true;
406 }
407 
408 template <typename DataType>
410 {
411  return m_wrapperType;
412 }
413 
414 template <typename DataType>
415 std::tuple<unsigned int, unsigned int> NekMatrix<
416  DataType, StandardMatrixTag>::Advance(unsigned int curRow,
417  unsigned int curColumn) const
418 {
419  return Advance(curRow, curColumn, this->GetTransposeFlag());
420 }
421 
422 template <typename DataType>
423 std::tuple<unsigned int, unsigned int> NekMatrix<
424  DataType, StandardMatrixTag>::Advance(unsigned int curRow,
425  unsigned int curColumn,
426  char transpose) const
427 {
428  unsigned int numRows = this->GetTransposedRows(transpose);
429  unsigned int numColumns = this->GetTransposedColumns(transpose);
430 
431  switch (this->GetStorageType())
432  {
433  case eFULL:
434  return FullMatrixFuncs::Advance(numRows, numColumns, curRow,
435  curColumn);
436  break;
437  case eDIAGONAL:
438  return DiagonalMatrixFuncs::Advance(numRows, numColumns, curRow,
439  curColumn);
440  break;
441  case eUPPER_TRIANGULAR:
442  return UpperTriangularMatrixFuncs::Advance(numRows, numColumns,
443  curRow, curColumn);
444  break;
445 
446  case eLOWER_TRIANGULAR:
447  return LowerTriangularMatrixFuncs::Advance(numRows, numColumns,
448  curRow, curColumn);
449  break;
450 
451  case eSYMMETRIC:
453  return SymmetricMatrixFuncs::Advance(numRows, numColumns, curRow,
454  curColumn);
455  break;
456  case eBANDED:
457  return BandedMatrixFuncs::Advance(numRows, numColumns, curRow,
458  curColumn);
459  break;
460  case eSYMMETRIC_BANDED:
462  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
463  break;
465  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
466  break;
468  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
469  break;
470 
471  default:
472  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
473  }
474  return std::tuple<unsigned int, unsigned int>(curRow, curColumn);
475 }
476 
477 template <typename DataType>
480 {
482 }
483 
484 template <typename DataType>
485 std::shared_ptr<NekMatrix<DataType, StandardMatrixTag>> NekMatrix<
486  DataType, StandardMatrixTag>::
487  CreateWrapper(
488  const std::shared_ptr<NekMatrix<DataType, StandardMatrixTag>> &rhs)
489 {
490  return std::shared_ptr<NekMatrix<DataType, StandardMatrixTag>>(
492 }
493 
494 template <typename DataType>
497  PointerWrapper wrapperType)
498  : BaseType(rhs), m_data(), m_wrapperType(wrapperType),
499  m_numberOfSuperDiagonals(rhs.m_numberOfSuperDiagonals),
500  m_numberOfSubDiagonals(rhs.m_numberOfSubDiagonals)
501 {
502  PerformCopyConstruction(rhs);
503 }
504 
505 template <typename DataType>
507 {
508  return m_data;
509 }
510 
511 template <typename DataType>
513 {
514  if (m_wrapperType == eCopy)
515  {
516  unsigned int requiredStorageSize = GetRequiredStorageSize();
517  if (m_data.size() > requiredStorageSize)
518  {
519  Array<OneD, DataType> newArray(requiredStorageSize);
520  CopyArrayN(m_data, newArray, requiredStorageSize);
521  m_data = newArray;
522  }
523  }
524  else if (m_wrapperType == eWrapper)
525  {
526  ASSERTL0(true, "Can't call RemoveExcessCapacity on a wrapped matrix.");
527  }
528 }
529 
530 template <typename DataType>
532  unsigned int requiredStorageSize)
533 {
534  if (m_wrapperType == eCopy)
535  {
536  if (m_data.size() < requiredStorageSize)
537  {
538  Array<OneD, DataType> newData(requiredStorageSize);
539  std::copy(m_data.data(), m_data.data() + m_data.size(),
540  newData.data());
541  m_data = newData;
542  }
543  }
544  else if (m_wrapperType == eWrapper)
545  {
546  // If the current matrix is wrapped, then just copy over the top,
547  // but the sizes of the two matrices must be the same.
548  ASSERTL0(
549  m_data.size() >= requiredStorageSize,
550  "Wrapped NekMatrices must have the same dimension in operator=");
551  }
552 }
553 
554 template <typename DataType>
556 {
557  unsigned int requiredStorageSize = GetRequiredStorageSize();
558  ResizeDataArrayIfNeeded(requiredStorageSize);
559 }
560 
561 template <typename DataType>
563  const ThisType &rhs)
564 {
565  if (m_wrapperType == eWrapper)
566  {
567  m_data = rhs.m_data;
568  }
569  else
570  {
571  m_data = Array<OneD, DataType>(GetRequiredStorageSize());
572  CopyArrayN(rhs.m_data, m_data, m_data.size());
573  }
574 }
575 
576 template <typename DataType>
577 typename boost::call_traits<DataType>::value_type NekMatrix<
578  DataType, StandardMatrixTag>::v_GetValue(unsigned int row,
579  unsigned int column) const
580 {
582 }
583 
584 template <typename DataType>
586 {
588 }
589 
590 template <typename DataType>
592  unsigned int row, unsigned int column,
593  typename boost::call_traits<DataType>::const_reference d)
594 {
596 }
597 
598 template <typename DataType>
600  unsigned int cols)
601 {
602  this->Resize(rows, cols);
603 
604  // Some places in Nektar++ access the matrix data array and
605  // use size() to see how big it is. When using
606  // expression templates, the data array's capacity is often larger
607  // than the actual number of elements, so this statement is
608  // required to report the correct number of elements.
609  this->GetData().ChangeSize(this->GetRequiredStorageSize());
610  ASSERTL0(this->GetRequiredStorageSize() <= this->GetData().size(),
611  "Can't resize matrices if there is not enough capacity.");
612 }
613 
614 template <typename DataType>
616  DataType, StandardMatrixTag>::operator()(unsigned int row,
617  unsigned int column)
618 {
619  ASSERTL2(row < this->GetRows(),
620  std::string("Row ") + std::to_string(row) +
621  std::string(" requested in a matrix with a maximum of ") +
622  std::to_string(this->GetRows()) + std::string(" rows"));
623  ASSERTL2(column < this->GetColumns(),
624  std::string("Column ") + std::to_string(column) +
625  std::string(" requested in a matrix with a maximum of ") +
626  std::to_string(this->GetColumns()) + std::string(" columns"));
627 
628  return (*this)(row, column, this->GetTransposeFlag());
629 }
630 
631 template <typename DataType>
633  DataType, StandardMatrixTag>::operator()(unsigned int row,
634  unsigned int column,
635  char transpose)
636 {
637  unsigned int index = this->CalculateIndex(row, column, transpose);
638  if (index != std::numeric_limits<unsigned int>::max())
639  {
640  return Proxy(this->GetData()[index]);
641  }
642  else
643  {
644  return Proxy();
645  }
646 }
647 
648 template <typename DataType>
650  unsigned int row, unsigned int column,
651  typename boost::call_traits<DataType>::const_reference d)
652 {
653  ASSERTL2(row < this->GetRows(),
654  std::string("Row ") + std::to_string(row) +
655  std::string(" requested in a matrix with a maximum of ") +
656  std::to_string(this->GetRows()) + std::string(" rows"));
657  ASSERTL2(column < this->GetColumns(),
658  std::string("Column ") + std::to_string(column) +
659  std::string(" requested in a matrix with a maximum of ") +
660  std::to_string(this->GetColumns()) + std::string(" columns"));
661  SetValue(row, column, d, this->GetTransposeFlag());
662 }
663 
664 template <typename DataType>
666  unsigned int row, unsigned int column,
667  typename boost::call_traits<DataType>::const_reference d, char transpose)
668 {
669  unsigned int index = this->CalculateIndex(row, column, transpose);
670  if (index != std::numeric_limits<unsigned int>::max())
671  {
672  this->GetData()[index] = d;
673  }
674  else
675  {
676  NEKERROR(
678  "Can't assign values into zeroed elements of a special array.");
679  }
680 }
681 
682 template <typename DataType>
684 {
685  return this->GetData();
686 }
687 
688 template <typename DataType>
690 {
691  return this->GetData().data();
692 }
693 
694 template <typename DataType>
696  DataType, StandardMatrixTag>::begin()
697 {
698  return begin(this->GetTransposeFlag());
699 }
700 
701 template <typename DataType>
703  DataType, StandardMatrixTag>::begin(char transpose)
704 {
705  if (transpose == 'N')
706  {
707  return iterator(this->GetData().data(),
708  this->GetData().data() + this->GetData().size());
709  }
710  else
711  {
712  return iterator(this, transpose);
713  }
714 }
715 
716 template <typename DataType>
718  DataType, StandardMatrixTag>::end()
719 {
720  return end(this->GetTransposeFlag());
721 }
722 
723 template <typename DataType>
725  DataType, StandardMatrixTag>::end(char transpose)
726 {
727  if (transpose == 'N')
728  {
729  return iterator(this->GetData().data(),
730  this->GetData().data() + this->GetData().size(), true);
731  }
732  else
733  {
734  return iterator(this, transpose, true);
735  }
736 }
737 
738 template <typename DataType>
740  void)
741 {
742  DataType returnval;
743  int nvals = this->GetColumns();
744  Array<OneD, DataType> EigValReal(nvals);
745  Array<OneD, DataType> EigValImag(nvals);
746  Array<OneD, DataType> Evecs;
747 
748  EigenSolve(EigValReal, EigValImag, Evecs);
749 
750  Vmath::Vmul(nvals, EigValReal, 1, EigValReal, 1, EigValReal, 1);
751  Vmath::Vmul(nvals, EigValImag, 1, EigValImag, 1, EigValImag, 1);
752  Vmath::Vadd(nvals, EigValReal, 1, EigValImag, 1, EigValReal, 1);
753 
754  returnval = sqrt(Vmath::Vmax(nvals, EigValReal, 1) /
755  Vmath::Vmin(nvals, EigValReal, 1));
756 
757  return returnval;
758 }
759 
760 template <typename DataType>
762  Array<OneD, DataType> &EigValReal, Array<OneD, DataType> &EigValImag,
763  Array<OneD, DataType> &EigVecs)
764 {
765  ASSERTL0(this->GetRows() == this->GetColumns(),
766  "Only square matrices can be called");
767 
768  switch (this->GetType())
769  {
770  case eFULL:
771  FullMatrixFuncs::EigenSolve(this->GetRows(), this->GetData(),
772  EigValReal, EigValImag, EigVecs);
773  break;
774  case eDIAGONAL:
775  Vmath::Vcopy(this->GetRows(), &(this->GetData())[0], 1,
776  &EigValReal[0], 1);
777  Vmath::Zero(this->GetRows(), &EigValImag[0], 1);
778  break;
779  case eUPPER_TRIANGULAR:
780  case eLOWER_TRIANGULAR:
781  case eSYMMETRIC:
782  case eBANDED:
783  case eSYMMETRIC_BANDED:
786  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
787  break;
788 
789  default:
790  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
791  }
792 }
793 
794 template <typename DataType>
796 {
797  ASSERTL0(this->GetRows() == this->GetColumns(),
798  "Only square matrices can be inverted.");
799  ASSERTL0(this->GetTransposeFlag() == 'N',
800  "Only untransposed matrices may be inverted.");
801 
802  switch (this->GetType())
803  {
804  case eFULL:
805  FullMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
806  this->GetData(), this->GetTransposeFlag());
807  break;
808  case eDIAGONAL:
809  DiagonalMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
810  this->GetData());
811  break;
812  case eSYMMETRIC:
813  SymmetricMatrixFuncs::Invert(this->GetRows(), this->GetColumns(),
814  this->GetData());
815  break;
816  case eUPPER_TRIANGULAR:
817  case eLOWER_TRIANGULAR:
818  case eBANDED:
819  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
820  break;
821  case eSYMMETRIC_BANDED:
822  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
823  break;
825  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
826  break;
828  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
829  break;
830 
831  default:
832  NEKERROR(ErrorUtil::efatal, "Unhandled matrix type");
833  }
834 }
835 
836 template <typename DataType>
838 {
839  if (m_tempSpace.capacity() == 0)
840  {
841  m_tempSpace = Array<OneD, DataType>(this->GetData().capacity());
842  }
843  return m_tempSpace;
844 }
845 
846 template <typename DataType>
848 {
849  std::swap(m_tempSpace, this->GetData());
850 }
851 
852 template <typename DataType>
854  DataType, StandardMatrixTag>::operator-() const
855 {
857  NegateInPlace(result);
858  return result;
859 }
860 
861 template <typename DataType>
864 {
865  for (unsigned int i = 0; i < this->GetPtr().size(); ++i)
866  {
867  this->GetPtr()[i] *= s;
868  }
869  return *this;
870 }
871 
872 template <typename DataType>
875 {
877  rhs.GetRows(), rhs.GetColumns(), rhs.GetPtr(), eWrapper, rhs.GetType(),
879  result.Transpose();
880  return result;
881 }
882 
885 
887 
890 
892 
893 template <typename DataType>
895 {
896  for (unsigned int i = 0; i < m.GetRows(); ++i)
897  {
898  for (unsigned int j = 0; j < m.GetColumns(); ++j)
899  {
900  m(i, j) *= -1.0;
901  }
902  }
903 }
904 
909 
910 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define LIB_UTILITIES_EXPORT
bool operator==(const VertexSharedPtr &v1, const VertexSharedPtr &v2)
Define comparison operator for the vertex struct.
Definition: VtkToFld.cpp:61
1D Array of constant elements with garbage collection and bounds checking.
Definition: SharedArray.hpp:58
size_type size() const
Returns the array's size.
MatrixStorage GetType() const
Definition: MatrixBase.hpp:64
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:124
unsigned int GetRows() const
Definition: MatrixBase.cpp:61
static unsigned int GetRequiredStorageSize(MatrixStorage type, unsigned int rows, unsigned int columns, unsigned int subDiags=0, unsigned int superDiags=0)
Definition: MatrixBase.cpp:182
char GetTransposeFlag() const
Definition: MatrixBase.cpp:118
unsigned int GetColumns() const
Definition: MatrixBase.cpp:80
Matrix< DataType > & operator=(const Matrix< DataType > &rhs)
Definition: MatrixBase.cpp:319
iterator_impl< const DataType, const ThisType > const_iterator
iterator_impl< DataType, ThisType > iterator
const Array< OneD, const DataType > & GetPtr() const
def copy(self)
Definition: pycml.py:2663
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
@ eVECTOR_WRAPPER
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
void NegateInPlace(NekVector< DataType > &v)
Definition: NekVector.cpp:1160
@ eLOWER_TRIANGULAR_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
@ eUPPER_TRIANGULAR_BANDED
void CopyArrayN(const Array< OneD, ConstDataType > &source, Array< OneD, DataType > &dest, typename Array< OneD, DataType >::size_type n)
PointerWrapper
Specifies if the pointer passed to a NekMatrix or NekVector is copied into an internal representation...
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:209
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:1050
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:359
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
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:945
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
static std::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static std::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:285
static void EigenSolve(unsigned int n, const Array< OneD, const DataType > &A, Array< OneD, DataType > &EigValReal, Array< OneD, DataType > &EigValImag, Array< OneD, DataType > &EigVecs)
Definition: MatrixFuncs.h:159
static std::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, const char transpose)
Definition: MatrixFuncs.h:91
static std::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')
static std::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:230
static std::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)