My Project
Matrix.hpp
1 //===========================================================================
2 //
3 // File: Matrix.hpp
4 //
5 // Created: Tue Jun 30 11:25:46 2009
6 //
7 // Author(s): Bård Skaflestad <bard.skaflestad@sintef.no>
8 // Atgeirr F Rasmussen <atgeirr@sintef.no>
9 //
10 // $Date$
11 //
12 // $Revision$
13 //
14 //===========================================================================
15 
16 /*
17  Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18  Copyright 2009, 2010 Statoil ASA.
19 
20  This file is part of The Open Reservoir Simulator Project (OpenRS).
21 
22  OpenRS is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OpenRS is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPENRS_MATRIX_HEADER
37 #define OPENRS_MATRIX_HEADER
38 
39 #include <algorithm>
40 #include <ostream>
41 #include <vector>
42 
43 #include <dune/common/fvector.hh>
44 #include <opm/common/ErrorMacros.hpp>
45 
46 #include <opm/porsol/common/fortran.hpp>
47 #include <opm/porsol/common/blas_lapack.hpp>
48 
49 namespace Opm {
50 
51  // ----------------------------------------------------------------------
52  // FullMatrix storage policies.
53  //
54 
62  template<typename T>
63  class OwnData {
64  public:
72  T& operator[](int i) { return data_[i]; }
73  const T& operator[](int i) const { return data_[i]; }
74 
75 
79  int size() const { return data_.size(); }
80 
84  T* data() { return &data_[0]; }
85  const T* data() const { return &data_[0]; }
86 
87  protected:
99  OwnData(int sz, const T* data)
100  {
101  if (data) {
102  data_.assign(data, data + sz);
103  } else {
104  data_.resize(sz);
105  }
106  }
107 
108  private:
109  std::vector<T> data_;
110  };
111 
119  template<typename T>
120  class SharedData {
121  public:
129  T& operator[](int i) { return data_[i]; }
130  const T& operator[](int i) const { return data_[i]; }
131 
135  int size() const { return sz_; }
136 
140  T* data() { return data_; }
141  const T* data() const { return data_; }
142 
143  protected:
153  SharedData(int sz, T* data)
154  : sz_(sz), data_(data)
155  {
156  assert ((sz == 0) == (data == 0));
157  }
158 
159  private:
160  int sz_;
161  T* data_;
162  };
163 
171  template<typename T>
173  public:
181  const T& operator[](int i) const { return data_[i]; }
182 
186  int size() const { return sz_; }
187 
191  const T* data() const { return data_; }
192 
193  protected:
202  ImmutableSharedData(int sz, const T* data)
203  : sz_(sz), data_(data)
204  {
205  assert (data_ != 0);
206  }
207 
208  private:
209  int sz_;
210  const T* data_;
211  };
212 
213 
214 
215 
216 
217  // ----------------------------------------------------------------------
218  // FullMatrix ordering policies.
219  //
220 
223  class OrderingBase {
224  public:
230  int numRows() const { return rows_; }
231 
237  int numCols() const { return cols_; }
238 
239  protected:
242  OrderingBase()
243  : rows_(0), cols_(0)
244  {}
245 
254  OrderingBase(int rows, int cols)
255  : rows_(rows), cols_(cols)
256  {}
257 
258  private:
259  int rows_, cols_;
260  };
261 
262 
266  class COrdering : public OrderingBase {
267  public:
275  int leadingDimension() const { return numCols(); }
276 
290  int idx(int row, int col) const
291  {
292  assert ((0 <= row) && (row < numRows()));
293  assert ((0 <= col) && (col < numCols()));
294 
295  return row*numCols() + col;
296  }
297 
298  protected:
300  COrdering()
301  : OrderingBase()
302  {}
303 
312  COrdering(int rows, int cols)
313  : OrderingBase(rows, cols)
314  {}
315  };
316 
317 
321  class FortranOrdering : public OrderingBase {
322  public:
330  int leadingDimension() const { return numRows(); }
331 
345  int idx(int row, int col) const
346  {
347  assert ((0 <= row) && (row < numRows()));
348  assert ((0 <= col) && (col < numCols()));
349 
350  return row + col*numRows();
351  }
352 
353  protected:
355  FortranOrdering()
356  : OrderingBase()
357  {}
358 
367  FortranOrdering(int rows, int cols)
368  : OrderingBase(rows, cols)
369  {}
370  };
371 
372 
373 
374 
375  // ----------------------------------------------------------------------
376  // Class FullMatrix.
377  //
378 
379 
394  template<typename T,
395  template<typename> class StoragePolicy,
396  class OrderingPolicy>
397  class FullMatrix : private StoragePolicy<T>,
398  private OrderingPolicy
399  {
400  public:
403  FullMatrix()
404  : StoragePolicy<T>(0, 0),
405  OrderingPolicy()
406  {}
407 
425  template <typename DataPointer>
426  FullMatrix(int rows, int cols, DataPointer data_arg)
427  : StoragePolicy<T>(rows * cols, data_arg),
428  OrderingPolicy(rows, cols)
429  {}
430 
439  template <template<typename> class OtherSP>
440  explicit FullMatrix(const FullMatrix<T, OtherSP, OrderingPolicy>& m)
441  : StoragePolicy<T>(m.numRows()*m.numCols(), m.data()),
442  OrderingPolicy(m.numRows(), m.numCols())
443  {
444  }
445 
457  template <template<typename> class OtherSP, class OtherOP>
458  FullMatrix& operator=(const FullMatrix<T, OtherSP, OtherOP>& m)
459  {
460  assert(numRows() == m.numRows());
461  assert(numCols() == m.numCols());
462  for (int r = 0; r < numRows(); ++r) {
463  for (int c = 0; c < numCols(); ++c) {
464  this->operator()(r, c) = m(r,c);
465  }
466  }
467  return *this;
468  }
469 
481  template <template<typename> class OtherSP, class OtherOP>
482  bool operator==(const FullMatrix<T, OtherSP, OtherOP>& m)
483  {
484  if (numRows() != m.numRows() || numCols() != m.numCols()) {
485  return false;
486  }
487  for (int r = 0; r < numRows(); ++r) {
488  for (int c = 0; c < numCols(); ++c) {
489  if (this->operator()(r,c) != m(r,c)) {
490  return false;
491  }
492  }
493  }
494  return true;
495  }
496 
505  template <template<typename> class OtherSP>
506  void operator+= (const FullMatrix<T, OtherSP, OrderingPolicy>& m)
507  {
508  assert(numRows() == m.numRows() && numCols() == m.numCols());
509  std::transform(data(), data() + this->size(),
510  m.data(), data(), std::plus<T>());
511  }
512 
518  void operator*= (const T& scalar)
519  {
520  std::transform(data(), data() + this->size(),
521  data(), [scalar](const T& input) { return input*scalar; } );
522  }
523 
526  typedef T value_type;
527 
528  using StoragePolicy<T>::data;
529  using OrderingPolicy::numRows;
530  using OrderingPolicy::numCols;
531  using OrderingPolicy::leadingDimension;
532 
544  value_type& operator()(int row, int col)
545  {
546  return this->operator[](this->idx(row, col));
547  }
548 
560  const value_type& operator()(int row, int col) const
561  {
562  return this->operator[](this->idx(row, col));
563  }
564  };
565 
566 
567 
568  // ----------------------------------------------------------------------
569  // FullMatrix operations.
570  //
571 
572 
573  // Convenience typedefs.
574 
579  typedef FullMatrix<double, OwnData, COrdering> OwnCMatrix;
580  typedef FullMatrix<double, SharedData, COrdering> SharedCMatrix;
581  typedef const FullMatrix<double, ImmutableSharedData, COrdering> ImmutableCMatrix;
582 
583 
588  typedef FullMatrix<double, OwnData, FortranOrdering> OwnFortranMatrix;
589  typedef FullMatrix<double, SharedData, FortranOrdering> SharedFortranMatrix;
590  typedef const FullMatrix<double, ImmutableSharedData, FortranOrdering> ImmutableFortranMatrix;
591 
592 
601  template<class Matrix>
602  void zero(Matrix& A)
603  {
604  std::fill_n(A.data(), A.numRows() * A.numCols(),
605  typename Matrix::value_type(0.0));
606  }
607 
616  template<class Matrix>
617  void eye(Matrix& A)
618  {
619  zero(A);
620  for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
621  A(i, i) = 1.0;
622  }
623  }
624 
636  template<class Matrix>
637  typename Matrix::value_type
638  trace(const Matrix& A)
639  {
640  typename Matrix::value_type ret(0);
641 
642  for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
643  ret += A(i,i);
644  }
645  return ret;
646  }
647 
648 
666  template<class Matrix, int rows>
667  Dune::FieldVector<typename Matrix::value_type, rows>
668  prod(const Matrix& A, const Dune::FieldVector<typename Matrix::value_type,rows>& x)
669  {
670  const int cols = rows;
671  assert (A.numRows() == rows);
672  assert (A.numCols() == cols);
673 
674  Dune::FieldVector<typename Matrix::value_type, rows> res(0.0);
675  for (int c = 0; c < cols; ++c) {
676  for (int r = 0; r < rows; ++r) {
677  res[r] += A(r, c)*x[c];
678  }
679  }
680  return res;
681  }
682 
690  template<class Matrix1, class Matrix2, class MutableMatrix>
691  void prod(const Matrix1& A, const Matrix2& B, MutableMatrix& C)
692  {
693  int result_rows = A.numRows();
694  int result_cols = B.numCols();
695  int inner_dim = A.numCols();
696  assert (inner_dim == B.numRows());
697  assert(C.numRows() == result_rows);
698  assert(C.numCols() == result_cols);
699 
700  for (int c = 0; c < result_cols; ++c) {
701  for (int r = 0; r < result_rows; ++r) {
702  C(r,c) = 0.0;
703  for (int i = 0; i < inner_dim; ++i) {
704  C(r,c) += A(r,i)*B(i,c);
705  }
706  }
707  }
708  }
709 
710 
728  template<typename T, template<typename> class StoragePolicy>
729  int orthogonalizeColumns(FullMatrix<T,StoragePolicy,FortranOrdering>& A)
730  {
731  typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
732 
733  static std::vector<value_type> tau;
734  static std::vector<value_type> work;
735 
736  if (int(tau .size()) < A.numCols()) tau .resize( A.numCols());
737  if (int(work.size()) < 64 * A.numRows()) work.resize(64 * A.numRows()); // 64 from ILAENV
738 
739  int info = 0;
740 
741  // Generate factorization. Matrix Q is implicitly represented.
742  Opm::BLAS_LAPACK::GEQRF(A.numRows(), A.numCols() ,
743  A.data() , A.leadingDimension() ,
744  &tau[0] , &work[0], int(work.size()),
745  info);
746 
747  if (info == 0)
748  {
749  // QR factorization successfully generated--extract
750  // explict representation of orthogonal matrix Q.
751  Opm::BLAS_LAPACK::ORGQR(A.numRows(), A.numCols(), A.numCols() ,
752  A.data() , A.leadingDimension() ,
753  &tau[0] , &work[0], int(work.size()),
754  info);
755  }
756 
757  return info;
758  }
759 
760 
779  template<typename T, template<typename> class StoragePolicy, class OrderingPolicy>
780  int invert(FullMatrix<T,StoragePolicy,OrderingPolicy>& A)
781  {
782  typedef typename FullMatrix<T,StoragePolicy,OrderingPolicy>::value_type value_type;
783 
784  assert (A.numRows() == A.numCols());
785 
786  std::vector<int> ipiv(A.numRows());
787  int info = 0;
788 
789  // Correct both for COrdering and FortranOrdering (inv(A)' == inv(A')).
790  Opm::BLAS_LAPACK::GETRF(A.numRows(), A.numCols(), A.data(),
791  A.leadingDimension(), &ipiv[0], info);
792 
793  if (info == 0) {
794  std::vector<value_type> work(A.numRows());
795 
796  Opm::BLAS_LAPACK::GETRI(A.numRows(), A.data(), A.leadingDimension(),
797  &ipiv[0], &work[0], int(work.size()), info);
798  }
799 
800  return info;
801  }
802 
803 
828  template<typename T, template<typename> class StoragePolicy>
829  void symmetricUpdate(const T& a1,
830  const FullMatrix<T,StoragePolicy,FortranOrdering>& A ,
831  const T& a2,
832  FullMatrix<T,StoragePolicy,FortranOrdering>& C )
833  {
834  Opm::BLAS_LAPACK::SYRK("Upper" , "No transpose" ,
835  C.numRows() , A.numCols() ,
836  a1, A.data(), A.leadingDimension(),
837  a2, C.data(), C.leadingDimension());
838 
839  // Account for SYRK (in this case) only updating the upper
840  // leading n-by-n submatrix.
841  //
842  for (int j = 0; j < C.numCols(); ++j) {
843  for (int i = j+1; i < C.numRows(); ++i) {
844  C(i,j) = C(j,i);
845  }
846  }
847  }
848 
849 
850  // B <- A*B*A'
851  // Assumes T is an arithmetic (floating point) type, and that A==A'.
856  template<typename T, template<typename> class StoragePolicy>
857  void symmetricUpdate(const FullMatrix<T,StoragePolicy,FortranOrdering>& A,
858  FullMatrix<T,StoragePolicy,FortranOrdering>& B)
859  {
860  typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
861 
862  // B <- A*B
863  Opm::BLAS_LAPACK::TRMM("Left" , "Upper", "No transpose", "Non-unit",
864  B.numRows(), B.numCols(), value_type(1.0),
865  A.data(), A.leadingDimension(),
866  B.data(), B.leadingDimension());
867 
868  // B <- B*A (== A * B_orig * A)
869  Opm::BLAS_LAPACK::TRMM("Right", "Upper", "No transpose", "Non-unit",
870  B.numRows(), B.numCols(), value_type(1.0),
871  A.data(), A.leadingDimension(),
872  B.data(), B.leadingDimension());
873 
874  // Account for TRMM (in this case) only updating the upper
875  // leading n-by-n submatrix.
876  //
877  for (int j = 0; j < B.numCols(); ++j) {
878  for (int i = j+1; i < B.numRows(); ++i) {
879  B(i,j) = B(j,i);
880  }
881  }
882  }
883 
884 
911  template<typename T ,
912  template<typename> class SP>
913  void vecMulAdd_N(const T& a1,
914  const FullMatrix<T,SP,FortranOrdering>& A ,
915  const std::vector<T>& x ,
916  const T& a2,
917  std::vector<T>& y)
918  {
919  assert(A.numRows() == y.size());
920  assert(A.numCols() == x.size());
921 
922  Opm::BLAS_LAPACK::GEMV("No Transpose",
923  A.numRows(), A.numCols(),
924  a1, A.data(), A.leadingDimension(),
925  &x[0], 1, a2, &y[0], 1);
926  }
927 
928 
955  template<typename T ,
956  template<typename> class SP>
957  void vecMulAdd_N(const T& a1,
958  const FullMatrix<T,SP,FortranOrdering>& A ,
959  const T* x ,
960  const T& a2,
961  T* y)
962  {
963  Opm::BLAS_LAPACK::GEMV("No Transpose",
964  A.numRows(), A.numCols(),
965  a1, A.data(), A.leadingDimension(),
966  x, 1, a2, y, 1);
967  }
968 
969 
997  template<typename T ,
998  template<typename> class SP>
999  void vecMulAdd_T(const T& a1,
1000  const FullMatrix<T,SP,FortranOrdering>& A ,
1001  const std::vector<T>& x ,
1002  const T& a2,
1003  std::vector<T>& y)
1004  {
1005  assert (A.numCols() == y.size());
1006  assert (A.numRows() == x.size());
1007 
1008  Opm::BLAS_LAPACK::GEMV("Transpose",
1009  A.numRows(), A.numCols(),
1010  a1, A.data(), A.leadingDimension(),
1011  &x[0], 1, a2, &y[0], 1);
1012  }
1013 
1014 
1042  template<typename T ,
1043  template<typename> class SP>
1044  void vecMulAdd_T(const T& a1,
1045  const FullMatrix<T,SP,FortranOrdering>& A ,
1046  const T* x ,
1047  const T& a2,
1048  T* y)
1049  {
1050  Opm::BLAS_LAPACK::GEMV("Transpose",
1051  A.numRows(), A.numCols(),
1052  a1, A.data(), A.leadingDimension(),
1053  x, 1, a2, y, 1);
1054  }
1055 
1056 
1084  template<typename T ,
1085  template<typename> class SP>
1086  void vecMulAdd_N(const T& a1,
1087  const FullMatrix<T,SP,COrdering>& A ,
1088  const T* x ,
1089  const T& a2,
1090  T* y)
1091  {
1092  typedef FullMatrix<T, ImmutableSharedData, FortranOrdering> FMAT;
1093 
1094  const FMAT At(A.numCols(), A.numRows(), A.data());
1095 
1096  vecMulAdd_T(a1, At, x, a2, y);
1097  }
1098 
1099 
1132  template<typename T ,
1133  template<typename> class SP1,
1134  template<typename> class SP2,
1135  template<typename> class SP3>
1136  void matMulAdd_NN(const T& a1,
1137  const FullMatrix<T,SP1,FortranOrdering>& A ,
1138  const FullMatrix<T,SP2,FortranOrdering>& B ,
1139  const T& a2,
1140  FullMatrix<T,SP3,FortranOrdering>& C)
1141  {
1142  assert(A.numRows() == C.numRows());
1143  assert(A.numCols() == B.numRows());
1144  assert(B.numCols() == C.numCols());
1145 
1146  int m = A.numRows(); // Number of *rows* in A
1147  int n = B.numCols(); // Number of *cols* in B
1148  int k = A.numCols(); // Number of *cols* in A (== numer of *rows* in B)
1149 
1150  Opm::BLAS_LAPACK::GEMM("No Transpose", "No Transpose", m, n, k,
1151  a1, A.data(), A.leadingDimension(),
1152  B.data(), B.leadingDimension(),
1153  a2, C.data(), C.leadingDimension());
1154  }
1155 
1156 
1190  template<typename T ,
1191  template<typename> class SP1,
1192  template<typename> class SP2,
1193  template<typename> class SP3>
1194  void matMulAdd_NT(const T& a1,
1195  const FullMatrix<T,SP1,FortranOrdering>& A ,
1196  const FullMatrix<T,SP2,FortranOrdering>& B ,
1197  const T& a2,
1198  FullMatrix<T,SP3,FortranOrdering>& C)
1199  {
1200  assert(A.numRows() == C.numRows());
1201  assert(B.numRows() == C.numCols());
1202  assert(A.numCols() == B.numCols());
1203 
1204  int m = A.numRows(); // Number of *rows* in A
1205  int n = B.numRows(); // Number of *cols* in B'
1206  int k = A.numCols(); // Number of *cols* in A (== numer of *rows* in B')
1207 
1208  Opm::BLAS_LAPACK::GEMM("No Transpose", "Transpose", m, n, k,
1209  a1, A.data(), A.leadingDimension(),
1210  B.data(), B.leadingDimension(),
1211  a2, C.data(), C.leadingDimension());
1212  }
1213 
1214 
1248  template<typename T ,
1249  template<typename> class SP1,
1250  template<typename> class SP2,
1251  template<typename> class SP3>
1252  void matMulAdd_TN(const T& a1,
1253  const FullMatrix<T,SP1,FortranOrdering>& A ,
1254  const FullMatrix<T,SP2,FortranOrdering>& B ,
1255  const T& a2,
1256  FullMatrix<T,SP3,FortranOrdering>& C)
1257  {
1258  assert (A.numCols() == C.numRows());
1259  assert (A.numRows() == B.numRows());
1260  assert (B.numCols() == C.numCols());
1261 
1262  int m = A.numCols(); // Number of *rows* in A'
1263  int n = B.numCols(); // Number of *cols* in B
1264  int k = A.numRows(); // Number of *cols* in A' (== numer of *rows* in B)
1265 
1266  Opm::BLAS_LAPACK::GEMM("Transpose", "No Transpose", m, n, k,
1267  a1, A.data(), A.leadingDimension(),
1268  B.data(), B.leadingDimension(),
1269  a2, C.data(), C.leadingDimension());
1270  }
1271 
1272 
1306  template<typename T ,
1307  template<typename> class SP1,
1308  template<typename> class SP2,
1309  template<typename> class SP3>
1310  void matMulAdd_NN(const T& a1,
1311  const FullMatrix<T,SP1,FortranOrdering>& A ,
1312  const FullMatrix<T,SP2,COrdering>& B ,
1313  const T& a2,
1314  FullMatrix<T,SP3,FortranOrdering>& C)
1315  {
1316  typedef typename FullMatrix<T,SP2,COrdering>::value_type value_type;
1317  typedef FullMatrix<value_type,ImmutableSharedData,FortranOrdering> FMat;
1318 
1319  const FMat Bt(B.numCols(), B.numRows(), B.data());
1320 
1321  matMulAdd_NT(a1, A, Bt, a2, C);
1322  }
1323 
1324 
1351  template<class charT, class traits,
1352  typename T, template<typename> class SP, class OP>
1353  std::basic_ostream<charT,traits>&
1354  operator<<(std::basic_ostream<charT,traits>& os,
1355  const FullMatrix<T,SP,OP>& A)
1356  {
1357  for (int i = 0; i < A.numRows(); ++i) {
1358  for (int j = 0; j < A.numCols(); ++j)
1359  os << A(i,j) << ' ';
1360  os << '\n';
1361  }
1362 
1363  return os;
1364  }
1365 } // namespace Opm
1366 #endif // OPENRS_MATRIX_HEADER
FullMatrix StoragePolicy which provides immutable object sharing semantics.
Definition: Matrix.hpp:172
ImmutableSharedData(int sz, const T *data)
Constructor.
Definition: Matrix.hpp:202
const T & operator[](int i) const
Storage element access.
Definition: Matrix.hpp:181
int size() const
Data size query.
Definition: Matrix.hpp:186
const T * data() const
Direct access to all data.
Definition: Matrix.hpp:191
FullMatrix StoragePolicy which provides object owning semantics.
Definition: Matrix.hpp:63
OwnData(int sz, const T *data)
Constructor.
Definition: Matrix.hpp:99
T & operator[](int i)
Storage element access.
Definition: Matrix.hpp:72
T * data()
Direct access to all data.
Definition: Matrix.hpp:84
int size() const
Data size query.
Definition: Matrix.hpp:79
FullMatrix StoragePolicy which provides object sharing semantics.
Definition: Matrix.hpp:120
T & operator[](int i)
Storage element access.
Definition: Matrix.hpp:129
int size() const
Data size query.
Definition: Matrix.hpp:135
T * data()
Direct access to all data.
Definition: Matrix.hpp:140
SharedData(int sz, T *data)
Constructor.
Definition: Matrix.hpp:153
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
void vecMulAdd_T(const T &a1, const FullMatrix< T, SP, FortranOrdering > &A, const std::vector< T > &x, const T &a2, std::vector< T > &y)
GEneral Matrix-Vector product (GAXPY operation).
Definition: Matrix.hpp:999
void matMulAdd_NN(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition: Matrix.hpp:1136
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition: Matrix.hpp:729
void symmetricUpdate(const T &a1, const FullMatrix< T, StoragePolicy, FortranOrdering > &A, const T &a2, FullMatrix< T, StoragePolicy, FortranOrdering > &C)
Symmetric, rank update of symmetric matrix.
Definition: Matrix.hpp:829
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:638
void matMulAdd_NT(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition: Matrix.hpp:1194
std::basic_ostream< charT, traits > & operator<<(std::basic_ostream< charT, traits > &os, const BCBase< T > &bc)
Stream insertion for BCBase.
Definition: BoundaryConditions.hpp:110
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:602
FullMatrix< double, OwnData, FortranOrdering > OwnFortranMatrix
Convenience typedefs for Fortran-ordered.
Definition: Matrix.hpp:588
int invert(FullMatrix< T, StoragePolicy, OrderingPolicy > &A)
Matrix inversion, .
Definition: Matrix.hpp:780
void vecMulAdd_N(const T &a1, const FullMatrix< T, SP, FortranOrdering > &A, const std::vector< T > &x, const T &a2, std::vector< T > &y)
GEneral Matrix-Vector product (GAXPY operation).
Definition: Matrix.hpp:913
void matMulAdd_TN(const T &a1, const FullMatrix< T, SP1, FortranOrdering > &A, const FullMatrix< T, SP2, FortranOrdering > &B, const T &a2, FullMatrix< T, SP3, FortranOrdering > &C)
GEneral Matrix-Matrix product update of other matrix.
Definition: Matrix.hpp:1252
Dune::FieldVector< typename Matrix::value_type, rows > prod(const Matrix &A, const Dune::FieldVector< typename Matrix::value_type, rows > &x)
Matrix applied to a vector.
Definition: Matrix.hpp:668
void eye(Matrix &A)
Make an identity.
Definition: Matrix.hpp:617
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition: Matrix.hpp:579