Go to the documentation of this file.
36#ifndef OPENRS_MATRIX_HEADER
37#define OPENRS_MATRIX_HEADER
43#include <dune/common/fvector.hh>
44#include <opm/common/ErrorMacros.hpp>
79 int size() const { return data_.size(); }
84 T* data() { return &data_[0]; }
85 const T* data() const { return &data_[0]; }
109 std::vector<T> data_;
135 int size() const { return sz_; }
141 const T* data() const { return data_; }
154 : sz_(sz), data_( data)
156 assert ((sz == 0) == ( data == 0));
186 int size() const { return sz_; }
191 const T* data() const { return data_; }
203 : sz_(sz), data_( data)
230 int numRows() const { return rows_; }
237 int numCols() const { return cols_; }
254 OrderingBase( int rows, int cols)
255 : rows_(rows), cols_(cols)
266 class COrdering : public OrderingBase {
275 int leadingDimension() const { return numCols(); }
290 int idx( int row, int col) const
292 assert ((0 <= row) && (row < numRows()));
293 assert ((0 <= col) && (col < numCols()));
295 return row*numCols() + col;
312 COrdering( int rows, int cols)
313 : OrderingBase(rows, cols)
321 class FortranOrdering : public OrderingBase {
330 int leadingDimension() const { return numRows(); }
345 int idx( int row, int col) const
347 assert ((0 <= row) && (row < numRows()));
348 assert ((0 <= col) && (col < numCols()));
350 return row + col*numRows();
367 FortranOrdering( int rows, int cols)
368 : OrderingBase(rows, cols)
395 template< typename> class StoragePolicy,
396 class OrderingPolicy>
397 class FullMatrix : private StoragePolicy<T>,
398 private OrderingPolicy
404 : StoragePolicy<T>(0, 0),
425 template < typename DataPo inter>
426 FullMatrix( int rows, int cols, DataPointer data_arg)
427 : StoragePolicy<T>(rows * cols, data_arg),
428 OrderingPolicy(rows, cols)
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())
457 template < template< typename> class OtherSP, class OtherOP>
458 FullMatrix& operator=( const FullMatrix<T, OtherSP, OtherOP>& m)
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);
481 template < template< typename> class OtherSP, class OtherOP>
482 bool operator==( const FullMatrix<T, OtherSP, OtherOP>& m)
484 if (numRows() != m.numRows() || numCols() != m.numCols()) {
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)) {
505 template < template< typename> class OtherSP>
506 void operator+= ( const FullMatrix<T, OtherSP, OrderingPolicy>& m)
508 assert(numRows() == m.numRows() && numCols() == m.numCols());
509 std::transform(data(), data() + this->size(),
510 m.data(), data(), std::plus<T>());
518 void operator*= ( const T& scalar)
520 std::transform(data(), data() + this->size(),
521 data(), [scalar]( const T& input) { return input*scalar; } );
526 typedef T value_type;
528 using StoragePolicy<T>::data;
529 using OrderingPolicy::numRows;
530 using OrderingPolicy::numCols;
531 using OrderingPolicy::leadingDimension;
544 value_type& operator()( int row, int col)
546 return this->operator[](this->idx(row, col));
560 const value_type& operator()( int row, int col) const
562 return this->operator[](this->idx(row, col));
601 template< class Matrix>
604 std::fill_n(A.data(), A.numRows() * A.numCols(),
605 typename Matrix::value_type(0.0));
616 template< class Matrix>
620 for ( int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
636 template< class Matrix>
637 typename Matrix::value_type
640 typename Matrix::value_type ret(0);
642 for ( int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
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)
670 const int cols = rows;
671 assert (A.numRows() == rows);
672 assert (A.numCols() == cols);
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];
690 template< class Matrix1, class Matrix2, class MutableMatrix>
691 void prod( const Matrix1& A, const Matrix2& B, MutableMatrix& C)
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);
700 for ( int c = 0; c < result_cols; ++c) {
701 for ( int r = 0; r < result_rows; ++r) {
703 for ( int i = 0; i < inner_dim; ++i) {
704 C(r,c) += A(r,i)*B(i,c);
728 template< typename T, template< typename> class StoragePolicy>
731 typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
733 static std::vector<value_type> tau;
734 static std::vector<value_type> work;
736 if ( int(tau .size()) < A.numCols()) tau .resize( A.numCols());
737 if ( int(work.size()) < 64 * A.numRows()) work.resize(64 * A.numRows());
743 A.data() , A.leadingDimension() ,
744 &tau[0] , &work[0], int(work.size()),
752 A.data() , A.leadingDimension() ,
753 &tau[0] , &work[0], int(work.size()),
779 template< typename T, template< typename> class StoragePolicy, class OrderingPolicy>
780 int invert(FullMatrix<T,StoragePolicy,OrderingPolicy>& A)
782 typedef typename FullMatrix<T,StoragePolicy,OrderingPolicy>::value_type value_type;
784 assert (A.numRows() == A.numCols());
786 std::vector<int> ipiv(A.numRows());
791 A.leadingDimension(), &ipiv[0], info);
794 std::vector<value_type> work(A.numRows());
797 &ipiv[0], &work[0], int(work.size()), info);
828 template< typename T, template< typename> class StoragePolicy>
830 const FullMatrix<T,StoragePolicy,FortranOrdering>& A ,
832 FullMatrix<T,StoragePolicy,FortranOrdering>& C )
835 C.numRows() , A.numCols() ,
836 a1, A.data(), A.leadingDimension(),
837 a2, C.data(), C.leadingDimension());
842 for ( int j = 0; j < C.numCols(); ++ j) {
843 for ( int i = j+1; i < C.numRows(); ++i) {
856 template< typename T, template< typename> class StoragePolicy>
858 FullMatrix<T,StoragePolicy,FortranOrdering>& B)
860 typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
864 B.numRows(), B.numCols(), value_type(1.0),
865 A.data(), A.leadingDimension(),
866 B.data(), B.leadingDimension());
870 B.numRows(), B.numCols(), value_type(1.0),
871 A.data(), A.leadingDimension(),
872 B.data(), B.leadingDimension());
877 for ( int j = 0; j < B.numCols(); ++ j) {
878 for ( int i = j+1; i < B.numRows(); ++i) {
911 template< typename T ,
912 template< typename> class SP>
914 const FullMatrix<T,SP,FortranOrdering>& A ,
915 const std::vector<T>& x ,
919 assert(A.numRows() == y.size());
920 assert(A.numCols() == x.size());
923 A.numRows(), A.numCols(),
924 a1, A.data(), A.leadingDimension(),
925 &x[0], 1, a2, &y[0], 1);
955 template< typename T ,
956 template< typename> class SP>
958 const FullMatrix<T,SP,FortranOrdering>& A ,
964 A.numRows(), A.numCols(),
965 a1, A.data(), A.leadingDimension(),
997 template< typename T ,
998 template< typename> class SP>
1000 const FullMatrix<T,SP,FortranOrdering>& A ,
1001 const std::vector<T>& x ,
1005 assert (A.numCols() == y.size());
1006 assert (A.numRows() == x.size());
1009 A.numRows(), A.numCols(),
1010 a1, A.data(), A.leadingDimension(),
1011 &x[0], 1, a2, &y[0], 1);
1042 template< typename T ,
1043 template< typename> class SP>
1045 const FullMatrix<T,SP,FortranOrdering>& A ,
1051 A.numRows(), A.numCols(),
1052 a1, A.data(), A.leadingDimension(),
1084 template< typename T ,
1085 template< typename> class SP>
1087 const FullMatrix<T,SP,COrdering>& A ,
1092 typedef FullMatrix<T, ImmutableSharedData, FortranOrdering> FMAT;
1094 const FMAT At(A.numCols(), A.numRows(), A.data());
1132 template< typename T ,
1133 template< typename> class SP1,
1134 template< typename> class SP2,
1135 template< typename> class SP3>
1137 const FullMatrix<T,SP1,FortranOrdering>& A ,
1138 const FullMatrix<T,SP2,FortranOrdering>& B ,
1140 FullMatrix<T,SP3,FortranOrdering>& C)
1142 assert(A.numRows() == C.numRows());
1143 assert(A.numCols() == B.numRows());
1144 assert(B.numCols() == C.numCols());
1146 int m = A.numRows();
1147 int n = B.numCols();
1148 int k = A.numCols();
1151 a1, A.data(), A.leadingDimension(),
1152 B.data(), B.leadingDimension(),
1153 a2, C.data(), C.leadingDimension());
1190 template< typename T ,
1191 template< typename> class SP1,
1192 template< typename> class SP2,
1193 template< typename> class SP3>
1195 const FullMatrix<T,SP1,FortranOrdering>& A ,
1196 const FullMatrix<T,SP2,FortranOrdering>& B ,
1198 FullMatrix<T,SP3,FortranOrdering>& C)
1200 assert(A.numRows() == C.numRows());
1201 assert(B.numRows() == C.numCols());
1202 assert(A.numCols() == B.numCols());
1204 int m = A.numRows();
1205 int n = B.numRows();
1206 int k = A.numCols();
1209 a1, A.data(), A.leadingDimension(),
1210 B.data(), B.leadingDimension(),
1211 a2, C.data(), C.leadingDimension());
1248 template< typename T ,
1249 template< typename> class SP1,
1250 template< typename> class SP2,
1251 template< typename> class SP3>
1253 const FullMatrix<T,SP1,FortranOrdering>& A ,
1254 const FullMatrix<T,SP2,FortranOrdering>& B ,
1256 FullMatrix<T,SP3,FortranOrdering>& C)
1258 assert (A.numCols() == C.numRows());
1259 assert (A.numRows() == B.numRows());
1260 assert (B.numCols() == C.numCols());
1262 int m = A.numCols();
1263 int n = B.numCols();
1264 int k = A.numRows();
1267 a1, A.data(), A.leadingDimension(),
1268 B.data(), B.leadingDimension(),
1269 a2, C.data(), C.leadingDimension());
1306 template< typename T ,
1307 template< typename> class SP1,
1308 template< typename> class SP2,
1309 template< typename> class SP3>
1311 const FullMatrix<T,SP1,FortranOrdering>& A ,
1312 const FullMatrix<T,SP2,COrdering>& B ,
1314 FullMatrix<T,SP3,FortranOrdering>& C)
1316 typedef typename FullMatrix<T,SP2,COrdering>::value_type value_type;
1317 typedef FullMatrix<value_type,ImmutableSharedData,FortranOrdering> FMat;
1319 const FMat Bt(B.numCols(), B.numRows(), B.data());
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)
1357 for ( int i = 0; i < A.numRows(); ++i) {
1358 for ( int j = 0; j < A.numCols(); ++ j)
1359 os << A(i, j) << ' ';
FullMatrix StoragePolicy which provides immutable object sharing semantics. Definition: Matrix.hpp:172
const T * data() const Direct access to all data. Definition: Matrix.hpp:191
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
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
const T * data() const Definition: Matrix.hpp:85
const T & operator[](int i) const Definition: Matrix.hpp:73
FullMatrix StoragePolicy which provides object sharing semantics. Definition: Matrix.hpp:120
const T & operator[](int i) const Definition: Matrix.hpp:130
T & operator[](int i) Storage element access. Definition: Matrix.hpp:129
int size() const Data size query. Definition: Matrix.hpp:135
SharedData(int sz, T *data) Constructor. Definition: Matrix.hpp:153
T * data() Direct access to all data. Definition: Matrix.hpp:140
const T * data() const Definition: Matrix.hpp:141
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).
min[0] Definition: elasticity_upscale_impl.hpp:146
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix A sparse matrix holding our operator. Definition: matrixops.hpp:27
int j Definition: elasticity_upscale_impl.hpp:658
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). Specifically, . 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. Specificlly, . Definition: Matrix.hpp:1136
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: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. Specifically, . Definition: Matrix.hpp:829
FullMatrix< double, SharedData, COrdering > SharedCMatrix Definition: Matrix.hpp:580
const FullMatrix< double, ImmutableSharedData, FortranOrdering > ImmutableFortranMatrix Definition: Matrix.hpp:590
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. Specificlly, . Definition: Matrix.hpp:1194
const FullMatrix< double, ImmutableSharedData, COrdering > ImmutableCMatrix Definition: Matrix.hpp:581
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
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:913
FullMatrix< double, SharedData, FortranOrdering > SharedFortranMatrix Definition: Matrix.hpp:589
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:1252
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
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
|