36#ifndef OPENRS_MATRIX_HEADER
37#define OPENRS_MATRIX_HEADER
43#include <dune/common/fvector.hh>
44#include <opm/common/ErrorMacros.hpp>
46#include <opm/porsol/common/fortran.hpp>
47#include <opm/porsol/common/blas_lapack.hpp>
73 const T&
operator[](
int i)
const {
return data_[i]; }
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_;
130 const T&
operator[](
int i)
const {
return data_[i]; }
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));
580 typedef FullMatrix<double, SharedData, COrdering> SharedCMatrix;
581 typedef const FullMatrix<double, ImmutableSharedData, COrdering> ImmutableCMatrix;
589 typedef FullMatrix<double, SharedData, FortranOrdering> SharedFortranMatrix;
590 typedef const FullMatrix<double, ImmutableSharedData, FortranOrdering> ImmutableFortranMatrix;
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());
742 Opm::BLAS_LAPACK::GEQRF(A.numRows(), A.numCols() ,
743 A.data() , A.leadingDimension() ,
744 &tau[0] , &work[0],
int(work.size()),
751 Opm::BLAS_LAPACK::ORGQR(A.numRows(), A.numCols(), A.numCols() ,
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());
790 Opm::BLAS_LAPACK::GETRF(A.numRows(), A.numCols(), A.data(),
791 A.leadingDimension(), &ipiv[0], info);
794 std::vector<value_type> work(A.numRows());
796 Opm::BLAS_LAPACK::GETRI(A.numRows(), A.data(), A.leadingDimension(),
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 )
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());
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;
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());
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());
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());
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);
955 template<
typename T ,
956 template<
typename>
class SP>
958 const FullMatrix<T,SP,FortranOrdering>& A ,
963 Opm::BLAS_LAPACK::GEMV(
"No Transpose",
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());
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);
1042 template<
typename T ,
1043 template<
typename>
class SP>
1045 const FullMatrix<T,SP,FortranOrdering>& A ,
1050 Opm::BLAS_LAPACK::GEMV(
"Transpose",
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();
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());
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();
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());
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();
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());
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
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
SharedData(int sz, T *data)
Constructor.
Definition Matrix.hpp:153
T * data()
Direct access to all data.
Definition Matrix.hpp:140
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
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).
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
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