Go to the documentation of this file.
36#ifndef OPENRS_MATRIX_HEADER
37#define OPENRS_MATRIX_HEADER
42#include <boost/bind.hpp>
44#include <dune/common/fvector.hh>
45#include <opm/common/ErrorMacros.hpp>
80 int size() const { return data_.size(); }
85 T* data() { return &data_[0]; }
86 const T* data() const { return &data_[0]; }
110 std::vector<T> data_;
136 int size() const { return sz_; }
142 const T* data() const { return data_; }
155 : sz_(sz), data_( data)
157 assert ((sz == 0) == ( data == 0));
187 int size() const { return sz_; }
192 const T* data() const { return data_; }
204 : sz_(sz), data_( data)
231 int numRows() const { return rows_; }
238 int numCols() const { return cols_; }
255 OrderingBase( int rows, int cols)
256 : rows_(rows), cols_(cols)
267 class COrdering : public OrderingBase {
276 int leadingDimension() const { return numCols(); }
291 int idx( int row, int col) const
293 assert ((0 <= row) && (row < numRows()));
294 assert ((0 <= col) && (col < numCols()));
296 return row*numCols() + col;
313 COrdering( int rows, int cols)
314 : OrderingBase(rows, cols)
322 class FortranOrdering : public OrderingBase {
331 int leadingDimension() const { return numRows(); }
346 int idx( int row, int col) const
348 assert ((0 <= row) && (row < numRows()));
349 assert ((0 <= col) && (col < numCols()));
351 return row + col*numRows();
368 FortranOrdering( int rows, int cols)
369 : OrderingBase(rows, cols)
396 template< typename> class StoragePolicy,
397 class OrderingPolicy>
398 class FullMatrix : private StoragePolicy<T>,
399 private OrderingPolicy
405 : StoragePolicy<T>(0, 0),
426 template < typename DataPo inter>
427 FullMatrix( int rows, int cols, DataPointer data)
428 : StoragePolicy<T>(rows * cols, data),
429 OrderingPolicy(rows, cols)
440 template < template< typename> class OtherSP>
441 explicit FullMatrix( const FullMatrix<T, OtherSP, OrderingPolicy>& m)
442 : StoragePolicy<T>(m.numRows()*m.numCols(), m.data()),
443 OrderingPolicy(m.numRows(), m.numCols())
458 template < template< typename> class OtherSP, class OtherOP>
459 FullMatrix& operator=( const FullMatrix<T, OtherSP, OtherOP>& m)
461 assert(numRows() == m.numRows());
462 assert(numCols() == m.numCols());
463 for ( int r = 0; r < numRows(); ++r) {
464 for ( int c = 0; c < numCols(); ++c) {
465 this->operator()(r, c) = m(r,c);
482 template < template< typename> class OtherSP, class OtherOP>
483 bool operator==( const FullMatrix<T, OtherSP, OtherOP>& m)
485 if (numRows() != m.numRows() || numCols() != m.numCols()) {
488 for ( int r = 0; r < numRows(); ++r) {
489 for ( int c = 0; c < numCols(); ++c) {
490 if (this-> operator()(r,c) != m(r,c)) {
506 template < template< typename> class OtherSP>
507 void operator+= ( const FullMatrix<T, OtherSP, OrderingPolicy>& m)
509 assert(numRows() == m.numRows() && numCols() == m.numCols());
510 std::transform(data(), data() + this->size(),
511 m.data(), data(), std::plus<T>());
519 void operator*= ( const T& scalar)
521 std::transform(data(), data() + this->size(),
522 data(), boost::bind(std::multiplies<T>(), _1, scalar));
527 typedef T value_type;
529 using StoragePolicy<T>::data;
530 using OrderingPolicy::numRows;
531 using OrderingPolicy::numCols;
532 using OrderingPolicy::leadingDimension;
545 value_type& operator()( int row, int col)
547 return this->operator[](this->idx(row, col));
561 const value_type& operator()( int row, int col) const
563 return this->operator[](this->idx(row, col));
602 template< class Matrix>
605 std::fill_n(A.data(), A.numRows() * A.numCols(),
606 typename Matrix::value_type(0.0));
617 template< class Matrix>
621 for ( int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
637 template< class Matrix>
638 typename Matrix::value_type
641 typename Matrix::value_type ret(0);
643 for ( int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
667 template< class Matrix, int rows>
668 Dune::FieldVector<typename Matrix::value_type, rows>
669 prod( const Matrix& A, const Dune::FieldVector<typename Matrix::value_type,rows>& x)
671 const int cols = rows;
672 assert (A.numRows() == rows);
673 assert (A.numCols() == cols);
675 Dune::FieldVector<typename Matrix::value_type, rows> res(0.0);
676 for ( int c = 0; c < cols; ++c) {
677 for ( int r = 0; r < rows; ++r) {
678 res[r] += A(r, c)*x[c];
691 template< class Matrix1, class Matrix2, class MutableMatrix>
692 void prod( const Matrix1& A, const Matrix2& B, MutableMatrix& C)
694 int result_rows = A.numRows();
695 int result_cols = B.numCols();
696 int inner_dim = A.numCols();
697 assert (inner_dim == B.numRows());
698 assert(C.numRows() == result_rows);
699 assert(C.numCols() == result_cols);
701 for ( int c = 0; c < result_cols; ++c) {
702 for ( int r = 0; r < result_rows; ++r) {
704 for ( int i = 0; i < inner_dim; ++i) {
705 C(r,c) += A(r,i)*B(i,c);
729 template< typename T, template< typename> class StoragePolicy>
732 typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
734 static std::vector<value_type> tau;
735 static std::vector<value_type> work;
737 if ( int(tau .size()) < A.numCols()) tau .resize( A.numCols());
738 if ( int(work.size()) < 64 * A.numRows()) work.resize(64 * A.numRows());
744 A.data() , A.leadingDimension() ,
745 &tau[0] , &work[0], int(work.size()),
753 A.data() , A.leadingDimension() ,
754 &tau[0] , &work[0], int(work.size()),
780 template< typename T, template< typename> class StoragePolicy, class OrderingPolicy>
781 int invert(FullMatrix<T,StoragePolicy,OrderingPolicy>& A)
783 typedef typename FullMatrix<T,StoragePolicy,OrderingPolicy>::value_type value_type;
785 assert (A.numRows() == A.numCols());
787 std::vector<int> ipiv(A.numRows());
792 A.leadingDimension(), &ipiv[0], info);
795 std::vector<value_type> work(A.numRows());
798 &ipiv[0], &work[0], int(work.size()), info);
829 template< typename T, template< typename> class StoragePolicy>
831 const FullMatrix<T,StoragePolicy,FortranOrdering>& A ,
833 FullMatrix<T,StoragePolicy,FortranOrdering>& C )
836 C.numRows() , A.numCols() ,
837 a1, A.data(), A.leadingDimension(),
838 a2, C.data(), C.leadingDimension());
843 for ( int j = 0; j < C.numCols(); ++j) {
844 for ( int i = j+1; i < C.numRows(); ++i) {
857 template< typename T, template< typename> class StoragePolicy>
859 FullMatrix<T,StoragePolicy,FortranOrdering>& B)
861 typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
865 B.numRows(), B.numCols(), value_type(1.0),
866 A.data(), A.leadingDimension(),
867 B.data(), B.leadingDimension());
871 B.numRows(), B.numCols(), value_type(1.0),
872 A.data(), A.leadingDimension(),
873 B.data(), B.leadingDimension());
878 for ( int j = 0; j < B.numCols(); ++j) {
879 for ( int i = j+1; i < B.numRows(); ++i) {
912 template< typename T ,
913 template< typename> class SP>
915 const FullMatrix<T,SP,FortranOrdering>& A ,
916 const std::vector<T>& x ,
920 assert(A.numRows() == y.size());
921 assert(A.numCols() == x.size());
924 A.numRows(), A.numCols(),
925 a1, A.data(), A.leadingDimension(),
926 &x[0], 1, a2, &y[0], 1);
956 template< typename T ,
957 template< typename> class SP>
959 const FullMatrix<T,SP,FortranOrdering>& A ,
965 A.numRows(), A.numCols(),
966 a1, A.data(), A.leadingDimension(),
998 template< typename T ,
999 template< typename> class SP>
1001 const FullMatrix<T,SP,FortranOrdering>& A ,
1002 const std::vector<T>& x ,
1006 assert (A.numCols() == y.size());
1007 assert (A.numRows() == x.size());
1010 A.numRows(), A.numCols(),
1011 a1, A.data(), A.leadingDimension(),
1012 &x[0], 1, a2, &y[0], 1);
1043 template< typename T ,
1044 template< typename> class SP>
1046 const FullMatrix<T,SP,FortranOrdering>& A ,
1052 A.numRows(), A.numCols(),
1053 a1, A.data(), A.leadingDimension(),
1085 template< typename T ,
1086 template< typename> class SP>
1088 const FullMatrix<T,SP,COrdering>& A ,
1093 typedef FullMatrix<T, ImmutableSharedData, FortranOrdering> FMAT;
1095 const FMAT At(A.numCols(), A.numRows(), A.data());
1133 template< typename T ,
1134 template< typename> class SP1,
1135 template< typename> class SP2,
1136 template< typename> class SP3>
1138 const FullMatrix<T,SP1,FortranOrdering>& A ,
1139 const FullMatrix<T,SP2,FortranOrdering>& B ,
1141 FullMatrix<T,SP3,FortranOrdering>& C)
1143 assert(A.numRows() == C.numRows());
1144 assert(A.numCols() == B.numRows());
1145 assert(B.numCols() == C.numCols());
1147 int m = A.numRows();
1148 int n = B.numCols();
1149 int k = A.numCols();
1152 a1, A.data(), A.leadingDimension(),
1153 B.data(), B.leadingDimension(),
1154 a2, C.data(), C.leadingDimension());
1191 template< typename T ,
1192 template< typename> class SP1,
1193 template< typename> class SP2,
1194 template< typename> class SP3>
1196 const FullMatrix<T,SP1,FortranOrdering>& A ,
1197 const FullMatrix<T,SP2,FortranOrdering>& B ,
1199 FullMatrix<T,SP3,FortranOrdering>& C)
1201 assert(A.numRows() == C.numRows());
1202 assert(B.numRows() == C.numCols());
1203 assert(A.numCols() == B.numCols());
1205 int m = A.numRows();
1206 int n = B.numRows();
1207 int k = A.numCols();
1210 a1, A.data(), A.leadingDimension(),
1211 B.data(), B.leadingDimension(),
1212 a2, C.data(), C.leadingDimension());
1249 template< typename T ,
1250 template< typename> class SP1,
1251 template< typename> class SP2,
1252 template< typename> class SP3>
1254 const FullMatrix<T,SP1,FortranOrdering>& A ,
1255 const FullMatrix<T,SP2,FortranOrdering>& B ,
1257 FullMatrix<T,SP3,FortranOrdering>& C)
1259 assert (A.numCols() == C.numRows());
1260 assert (A.numRows() == B.numRows());
1261 assert (B.numCols() == C.numCols());
1263 int m = A.numCols();
1264 int n = B.numCols();
1265 int k = A.numRows();
1268 a1, A.data(), A.leadingDimension(),
1269 B.data(), B.leadingDimension(),
1270 a2, C.data(), C.leadingDimension());
1307 template< typename T ,
1308 template< typename> class SP1,
1309 template< typename> class SP2,
1310 template< typename> class SP3>
1312 const FullMatrix<T,SP1,FortranOrdering>& A ,
1313 const FullMatrix<T,SP2,COrdering>& B ,
1315 FullMatrix<T,SP3,FortranOrdering>& C)
1317 typedef typename FullMatrix<T,SP2,COrdering>::value_type value_type;
1318 typedef FullMatrix<value_type,ImmutableSharedData,FortranOrdering> FMat;
1320 const FMat Bt(B.numCols(), B.numRows(), B.data());
1352 template< class charT, class traits,
1353 typename T, template< typename> class SP, class OP>
1354 std::basic_ostream<charT,traits>&
1355 operator<<(std::basic_ostream<charT,traits>& os,
1356 const FullMatrix<T,SP,OP>& A)
1358 for ( int i = 0; i < A.numRows(); ++i) {
1359 for ( int j = 0; j < A.numCols(); ++j)
1360 os << A(i,j) << ' ';
FullMatrix StoragePolicy which provides immutable object sharing semantics. Definition: Matrix.hpp:173
const T * data() const Direct access to all data. Definition: Matrix.hpp:192
ImmutableSharedData(int sz, const T *data) Constructor. Definition: Matrix.hpp:203
const T & operator[](int i) const Storage element access. Definition: Matrix.hpp:182
int size() const Data size query. Definition: Matrix.hpp:187
FullMatrix StoragePolicy which provides object owning semantics. Definition: Matrix.hpp:64
OwnData(int sz, const T *data) Constructor. Definition: Matrix.hpp:100
T & operator[](int i) Storage element access. Definition: Matrix.hpp:73
T * data() Direct access to all data. Definition: Matrix.hpp:85
int size() const Data size query. Definition: Matrix.hpp:80
const T * data() const Definition: Matrix.hpp:86
const T & operator[](int i) const Definition: Matrix.hpp:74
FullMatrix StoragePolicy which provides object sharing semantics. Definition: Matrix.hpp:121
const T & operator[](int i) const Definition: Matrix.hpp:131
T & operator[](int i) Storage element access. Definition: Matrix.hpp:130
int size() const Data size query. Definition: Matrix.hpp:136
SharedData(int sz, T *data) Constructor. Definition: Matrix.hpp:154
T * data() Direct access to all data. Definition: Matrix.hpp:141
const T * data() const Definition: Matrix.hpp:142
void GEMM(const char *transA, const char *transB, const int m, const int n, const int k, const T &a1, const T *A, const int ldA, const T *B, const int ldB, const T &a2, T *C, const int ldC) GEneral Matrix Matrix product (Level 3 BLAS).
void GEMV(const char *transA, const int m, const int n, const T &a1, const T *A, const int ldA, const T *x, const int incX, const T &a2, T *y, const int incY) GEneral Matrix Vector product (Level 2 BLAS).
void GEQRF(const int m, const int n, T *A, const int ld, T *tau, T *work, const int lwork, int &info) GEneral matrix QR Factorization (LAPACK)
void SYRK(const char *uplo, const char *trans, const int n, const int k, const T &a1, const T *A, const int ldA, const T &a2, T *C, const int ldC) SYmmetric Rank K update of symmetric matrix (Level 3 BLAS)
void GETRI(const int n, T *A, const int ld, const int *ipiv, T *work, int lwork, int &info) GEneral matrix TRiangular Inversion (LAPACK).
void TRMM(const char *side, const char *uplo, const char *transA, const char *diag, const int m, const int n, const T &a, const T *A, const int ldA, T *B, const int ldB) TRiangular Matrix Matrix product (Level 2 BLAS)
void GETRF(const int m, const int n, T *A, const int ld, int *ipiv, int &info) GEneral matrix TRiangular Factorization (LAPACK).
void ORGQR(const int m, const int n, const int k, T *A, const int ld, const T *tau, T *work, const int lwork, int &info) ORthogonal matrix Generator from QR factorization (LAPACK).
Definition: BlackoilFluid.hpp:32
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). Specifically, . Definition: Matrix.hpp:1000
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. Specificlly, . Definition: Matrix.hpp:1137
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A) Construct orthonormal basis for matrix range (i.e., column space). Based on a QR factorization of the... Definition: Matrix.hpp:730
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. Specifically, . Definition: Matrix.hpp:830
FullMatrix< double, SharedData, COrdering > SharedCMatrix Definition: Matrix.hpp:581
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix Definition: Matrix.hpp:591
Matrix::value_type trace(const Matrix &A) Compute matrix trace (i.e., sum(diag(A))). Definition: Matrix.hpp:639
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. Specificlly, . Definition: Matrix.hpp:1195
const FullMatrix< double, ImmutableSharedData, COrdering > ImmutableCMatrix Definition: Matrix.hpp:582
void zero(Matrix &A) Zero-fill a. Definition: Matrix.hpp:603
FullMatrix< double, OwnData, FortranOrdering > OwnFortranMatrix Convenience typedefs for Fortran-ordered. Definition: Matrix.hpp:589
int invert(FullMatrix< T, StoragePolicy, OrderingPolicy > &A) Matrix inversion, . Definition: Matrix.hpp:781
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 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). Specifically, . Definition: Matrix.hpp:914
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix Definition: Matrix.hpp:590
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. Specificlly, . Definition: Matrix.hpp:1253
void eye(Matrix &A) Make an identity. Definition: Matrix.hpp:618
FullMatrix< double, OwnData, COrdering > OwnCMatrix Convenience typedefs for C-ordered. Definition: Matrix.hpp:580
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:669
|