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>
47 #include <opm/porsol/common/fortran.hpp>
48 #include <opm/porsol/common/blas_lapack.hpp>
74 const T&
operator[](
int i)
const {
return data_[i]; }
80 int size()
const {
return data_.size(); }
85 T*
data() {
return &data_[0]; }
86 const T*
data()
const {
return &data_[0]; }
103 data_.assign(data, data + sz);
110 std::vector<T> data_;
131 const T&
operator[](
int i)
const {
return data_[i]; }
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));
581 typedef FullMatrix<double, SharedData, COrdering> SharedCMatrix;
582 typedef const FullMatrix<double, ImmutableSharedData, COrdering> ImmutableCMatrix;
590 typedef FullMatrix<double, SharedData, FortranOrdering> SharedFortranMatrix;
591 typedef const FullMatrix<double, ImmutableSharedData, FortranOrdering> ImmutableFortranMatrix;
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());
743 Opm::BLAS_LAPACK::GEQRF(A.numRows(), A.numCols() ,
744 A.data() , A.leadingDimension() ,
745 &tau[0] , &work[0], int(work.size()),
752 Opm::BLAS_LAPACK::ORGQR(A.numRows(), A.numCols(), A.numCols() ,
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());
791 Opm::BLAS_LAPACK::GETRF(A.numRows(), A.numCols(), A.data(),
792 A.leadingDimension(), &ipiv[0], info);
795 std::vector<value_type> work(A.numRows());
797 Opm::BLAS_LAPACK::GETRI(A.numRows(), A.data(), A.leadingDimension(),
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 )
835 Opm::BLAS_LAPACK::SYRK(
"Upper" ,
"No transpose" ,
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;
864 Opm::BLAS_LAPACK::TRMM(
"Left" ,
"Upper",
"No transpose",
"Non-unit",
865 B.numRows(), B.numCols(), value_type(1.0),
866 A.data(), A.leadingDimension(),
867 B.data(), B.leadingDimension());
870 Opm::BLAS_LAPACK::TRMM(
"Right",
"Upper",
"No transpose",
"Non-unit",
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());
923 Opm::BLAS_LAPACK::GEMV(
"No Transpose",
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 ,
964 Opm::BLAS_LAPACK::GEMV(
"No Transpose",
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());
1009 Opm::BLAS_LAPACK::GEMV(
"Transpose",
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 ,
1051 Opm::BLAS_LAPACK::GEMV(
"Transpose",
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();
1151 Opm::BLAS_LAPACK::GEMM(
"No Transpose",
"No Transpose", m, n, k,
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();
1209 Opm::BLAS_LAPACK::GEMM(
"No Transpose",
"Transpose", m, n, k,
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();
1267 Opm::BLAS_LAPACK::GEMM(
"Transpose",
"No Transpose", m, n, k,
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) <<
' ';
1367 #endif // OPENRS_MATRIX_HEADER
FullMatrix StoragePolicy which provides immutable object sharing semantics.
Definition: Matrix.hpp:173
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition: Matrix.hpp:580
T * data()
Direct access to all data.
Definition: Matrix.hpp:141
FullMatrix StoragePolicy which provides object owning semantics.
Definition: Matrix.hpp:64
int size() const
Data size query.
Definition: Matrix.hpp:136
int invert(FullMatrix< T, StoragePolicy, OrderingPolicy > &A)
Matrix inversion, .
Definition: Matrix.hpp:781
FullMatrix StoragePolicy which provides object sharing semantics.
Definition: Matrix.hpp:121
T & operator[](int i)
Storage element access.
Definition: Matrix.hpp:130
ImmutableSharedData(int sz, const T *data)
Constructor.
Definition: Matrix.hpp:203
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:1253
T & operator[](int i)
Storage element access.
Definition: Matrix.hpp:73
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:1000
int size() const
Data size query.
Definition: Matrix.hpp:187
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:1137
void eye(Matrix &A)
Make an identity.
Definition: Matrix.hpp:618
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition: Matrix.hpp:639
T * data()
Direct access to all data.
Definition: Matrix.hpp:85
OwnData(int sz, const T *data)
Constructor.
Definition: Matrix.hpp:100
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
const T * data() const
Direct access to all data.
Definition: Matrix.hpp:192
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:830
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:914
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition: Matrix.hpp:730
const T & operator[](int i) const
Storage element access.
Definition: Matrix.hpp:182
int size() const
Data size query.
Definition: Matrix.hpp:80
FullMatrix< double, OwnData, FortranOrdering > OwnFortranMatrix
Convenience typedefs for Fortran-ordered.
Definition: Matrix.hpp:589
SharedData(int sz, T *data)
Constructor.
Definition: Matrix.hpp:154
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
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:1195