Loading...
Searching...
No Matches
Matrix.hpp
//===========================================================================
//
// File: Matrix.hpp
//
// Created: Tue Jun 30 11:25:46 2009
//
// Author(s): Bård Skaflestad <bard.skaflestad@sintef.no>
// Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenRS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_MATRIX_HEADER
#define OPENRS_MATRIX_HEADER
#include <algorithm>
#include <ostream>
#include <vector>
#include <dune/common/fvector.hh>
#include <opm/common/ErrorMacros.hpp>
#include <opm/porsol/common/fortran.hpp>
#include <opm/porsol/common/blas_lapack.hpp>
namespace Opm {
// ----------------------------------------------------------------------
// FullMatrix storage policies.
//
template<typename T>
class OwnData {
public:
T& operator[](int i) { return data_[i]; }
const T& operator[](int i) const { return data_[i]; }
int size() const { return data_.size(); }
T* data() { return &data_[0]; }
const T* data() const { return &data_[0]; }
protected:
OwnData(int sz, const T* data)
{
if (data) {
data_.assign(data, data + sz);
} else {
data_.resize(sz);
}
}
private:
std::vector<T> data_;
};
template<typename T>
class SharedData {
public:
T& operator[](int i) { return data_[i]; }
const T& operator[](int i) const { return data_[i]; }
int size() const { return sz_; }
T* data() { return data_; }
const T* data() const { return data_; }
protected:
SharedData(int sz, T* data)
: sz_(sz), data_(data)
{
assert ((sz == 0) == (data == 0));
}
private:
int sz_;
T* data_;
};
template<typename T>
class ImmutableSharedData {
public:
const T& operator[](int i) const { return data_[i]; }
int size() const { return sz_; }
const T* data() const { return data_; }
protected:
ImmutableSharedData(int sz, const T* data)
: sz_(sz), data_(data)
{
assert (data_ != 0);
}
private:
int sz_;
const T* data_;
};
// ----------------------------------------------------------------------
// FullMatrix ordering policies.
//
class OrderingBase {
public:
int numRows() const { return rows_; }
int numCols() const { return cols_; }
protected:
OrderingBase()
: rows_(0), cols_(0)
{}
OrderingBase(int rows, int cols)
: rows_(rows), cols_(cols)
{}
private:
int rows_, cols_;
};
class COrdering : public OrderingBase {
public:
int leadingDimension() const { return numCols(); }
int idx(int row, int col) const
{
assert ((0 <= row) && (row < numRows()));
assert ((0 <= col) && (col < numCols()));
return row*numCols() + col;
}
protected:
COrdering()
: OrderingBase()
{}
COrdering(int rows, int cols)
: OrderingBase(rows, cols)
{}
};
class FortranOrdering : public OrderingBase {
public:
int leadingDimension() const { return numRows(); }
int idx(int row, int col) const
{
assert ((0 <= row) && (row < numRows()));
assert ((0 <= col) && (col < numCols()));
return row + col*numRows();
}
protected:
FortranOrdering()
: OrderingBase()
{}
FortranOrdering(int rows, int cols)
: OrderingBase(rows, cols)
{}
};
// ----------------------------------------------------------------------
// Class FullMatrix.
//
template<typename T,
template<typename> class StoragePolicy,
class OrderingPolicy>
class FullMatrix : private StoragePolicy<T>,
private OrderingPolicy
{
public:
FullMatrix()
: StoragePolicy<T>(0, 0),
OrderingPolicy()
{}
template <typename DataPointer>
FullMatrix(int rows, int cols, DataPointer data_arg)
: StoragePolicy<T>(rows * cols, data_arg),
OrderingPolicy(rows, cols)
{}
template <template<typename> class OtherSP>
explicit FullMatrix(const FullMatrix<T, OtherSP, OrderingPolicy>& m)
: StoragePolicy<T>(m.numRows()*m.numCols(), m.data()),
OrderingPolicy(m.numRows(), m.numCols())
{
}
template <template<typename> class OtherSP, class OtherOP>
FullMatrix& operator=(const FullMatrix<T, OtherSP, OtherOP>& m)
{
assert(numRows() == m.numRows());
assert(numCols() == m.numCols());
for (int r = 0; r < numRows(); ++r) {
for (int c = 0; c < numCols(); ++c) {
this->operator()(r, c) = m(r,c);
}
}
return *this;
}
template <template<typename> class OtherSP, class OtherOP>
bool operator==(const FullMatrix<T, OtherSP, OtherOP>& m) const
{
if (numRows() != m.numRows() || numCols() != m.numCols()) {
return false;
}
for (int r = 0; r < numRows(); ++r) {
for (int c = 0; c < numCols(); ++c) {
if (this->operator()(r,c) != m(r,c)) {
return false;
}
}
}
return true;
}
template <template<typename> class OtherSP>
void operator+= (const FullMatrix<T, OtherSP, OrderingPolicy>& m)
{
assert(numRows() == m.numRows() && numCols() == m.numCols());
std::transform(data(), data() + this->size(),
m.data(), data(), std::plus<T>());
}
void operator*= (const T& scalar)
{
std::transform(data(), data() + this->size(),
data(), [scalar](const T& input) { return input*scalar; } );
}
typedef T value_type;
using StoragePolicy<T>::data;
using OrderingPolicy::numRows;
using OrderingPolicy::numCols;
using OrderingPolicy::leadingDimension;
value_type& operator()(int row, int col)
{
return this->operator[](this->idx(row, col));
}
const value_type& operator()(int row, int col) const
{
return this->operator[](this->idx(row, col));
}
};
// ----------------------------------------------------------------------
// FullMatrix operations.
//
// Convenience typedefs.
typedef FullMatrix<double, OwnData, COrdering> OwnCMatrix;
typedef FullMatrix<double, SharedData, COrdering> SharedCMatrix;
typedef const FullMatrix<double, ImmutableSharedData, COrdering> ImmutableCMatrix;
typedef FullMatrix<double, OwnData, FortranOrdering> OwnFortranMatrix;
typedef FullMatrix<double, SharedData, FortranOrdering> SharedFortranMatrix;
typedef const FullMatrix<double, ImmutableSharedData, FortranOrdering> ImmutableFortranMatrix;
template<class Matrix>
void zero(Matrix& A)
{
std::fill_n(A.data(), A.numRows() * A.numCols(),
typename Matrix::value_type(0.0));
}
template<class Matrix>
void eye(Matrix& A)
{
zero(A);
for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
A(i, i) = 1.0;
}
}
template<class Matrix>
typename Matrix::value_type
trace(const Matrix& A)
{
typename Matrix::value_type ret(0);
for (int i = 0; i < std::min(A.numRows(), A.numCols()); ++i) {
ret += A(i,i);
}
return ret;
}
template<class Matrix, int rows>
Dune::FieldVector<typename Matrix::value_type, rows>
prod(const Matrix& A, const Dune::FieldVector<typename Matrix::value_type,rows>& x)
{
const int cols = rows;
assert (A.numRows() == rows);
assert (A.numCols() == cols);
Dune::FieldVector<typename Matrix::value_type, rows> res(0.0);
for (int c = 0; c < cols; ++c) {
for (int r = 0; r < rows; ++r) {
res[r] += A(r, c)*x[c];
}
}
return res;
}
template<class Matrix1, class Matrix2, class MutableMatrix>
void prod(const Matrix1& A, const Matrix2& B, MutableMatrix& C)
{
int result_rows = A.numRows();
int result_cols = B.numCols();
int inner_dim = A.numCols();
assert (inner_dim == B.numRows());
assert(C.numRows() == result_rows);
assert(C.numCols() == result_cols);
for (int c = 0; c < result_cols; ++c) {
for (int r = 0; r < result_rows; ++r) {
C(r,c) = 0.0;
for (int i = 0; i < inner_dim; ++i) {
C(r,c) += A(r,i)*B(i,c);
}
}
}
}
template<typename T, template<typename> class StoragePolicy>
int orthogonalizeColumns(FullMatrix<T,StoragePolicy,FortranOrdering>& A)
{
typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
static std::vector<value_type> tau;
static std::vector<value_type> work;
if (int(tau .size()) < A.numCols()) tau .resize( A.numCols());
if (int(work.size()) < 64 * A.numRows()) work.resize(64 * A.numRows()); // 64 from ILAENV
int info = 0;
// Generate factorization. Matrix Q is implicitly represented.
Opm::BLAS_LAPACK::GEQRF(A.numRows(), A.numCols() ,
A.data() , A.leadingDimension() ,
&tau[0] , &work[0], int(work.size()),
info);
if (info == 0)
{
// QR factorization successfully generated--extract
// explict representation of orthogonal matrix Q.
Opm::BLAS_LAPACK::ORGQR(A.numRows(), A.numCols(), A.numCols() ,
A.data() , A.leadingDimension() ,
&tau[0] , &work[0], int(work.size()),
info);
}
return info;
}
template<typename T, template<typename> class StoragePolicy, class OrderingPolicy>
int invert(FullMatrix<T,StoragePolicy,OrderingPolicy>& A)
{
typedef typename FullMatrix<T,StoragePolicy,OrderingPolicy>::value_type value_type;
assert (A.numRows() == A.numCols());
std::vector<int> ipiv(A.numRows());
int info = 0;
// Correct both for COrdering and FortranOrdering (inv(A)' == inv(A')).
Opm::BLAS_LAPACK::GETRF(A.numRows(), A.numCols(), A.data(),
A.leadingDimension(), &ipiv[0], info);
if (info == 0) {
std::vector<value_type> work(A.numRows());
Opm::BLAS_LAPACK::GETRI(A.numRows(), A.data(), A.leadingDimension(),
&ipiv[0], &work[0], int(work.size()), info);
}
return info;
}
template<typename T, template<typename> class StoragePolicy>
void symmetricUpdate(const T& a1,
const FullMatrix<T,StoragePolicy,FortranOrdering>& A ,
const T& a2,
FullMatrix<T,StoragePolicy,FortranOrdering>& C )
{
Opm::BLAS_LAPACK::SYRK("Upper" , "No transpose" ,
C.numRows() , A.numCols() ,
a1, A.data(), A.leadingDimension(),
a2, C.data(), C.leadingDimension());
// Account for SYRK (in this case) only updating the upper
// leading n-by-n submatrix.
//
for (int j = 0; j < C.numCols(); ++j) {
for (int i = j+1; i < C.numRows(); ++i) {
C(i,j) = C(j,i);
}
}
}
// B <- A*B*A'
// Assumes T is an arithmetic (floating point) type, and that A==A'.
template<typename T, template<typename> class StoragePolicy>
void symmetricUpdate(const FullMatrix<T,StoragePolicy,FortranOrdering>& A,
FullMatrix<T,StoragePolicy,FortranOrdering>& B)
{
typedef typename FullMatrix<T,StoragePolicy,FortranOrdering>::value_type value_type;
// B <- A*B
Opm::BLAS_LAPACK::TRMM("Left" , "Upper", "No transpose", "Non-unit",
B.numRows(), B.numCols(), value_type(1.0),
A.data(), A.leadingDimension(),
B.data(), B.leadingDimension());
// B <- B*A (== A * B_orig * A)
Opm::BLAS_LAPACK::TRMM("Right", "Upper", "No transpose", "Non-unit",
B.numRows(), B.numCols(), value_type(1.0),
A.data(), A.leadingDimension(),
B.data(), B.leadingDimension());
// Account for TRMM (in this case) only updating the upper
// leading n-by-n submatrix.
//
for (int j = 0; j < B.numCols(); ++j) {
for (int i = j+1; i < B.numRows(); ++i) {
B(i,j) = B(j,i);
}
}
}
template<typename T ,
template<typename> class SP>
void vecMulAdd_N(const T& a1,
const FullMatrix<T,SP,FortranOrdering>& A ,
const std::vector<T>& x ,
const T& a2,
std::vector<T>& y)
{
assert(A.numRows() == y.size());
assert(A.numCols() == x.size());
Opm::BLAS_LAPACK::GEMV("No Transpose",
A.numRows(), A.numCols(),
a1, A.data(), A.leadingDimension(),
&x[0], 1, a2, &y[0], 1);
}
template<typename T ,
template<typename> class SP>
void vecMulAdd_N(const T& a1,
const FullMatrix<T,SP,FortranOrdering>& A ,
const T* x ,
const T& a2,
T* y)
{
Opm::BLAS_LAPACK::GEMV("No Transpose",
A.numRows(), A.numCols(),
a1, A.data(), A.leadingDimension(),
x, 1, a2, y, 1);
}
template<typename T ,
template<typename> class SP>
void vecMulAdd_T(const T& a1,
const FullMatrix<T,SP,FortranOrdering>& A ,
const std::vector<T>& x ,
const T& a2,
std::vector<T>& y)
{
assert (A.numCols() == y.size());
assert (A.numRows() == x.size());
Opm::BLAS_LAPACK::GEMV("Transpose",
A.numRows(), A.numCols(),
a1, A.data(), A.leadingDimension(),
&x[0], 1, a2, &y[0], 1);
}
template<typename T ,
template<typename> class SP>
void vecMulAdd_T(const T& a1,
const FullMatrix<T,SP,FortranOrdering>& A ,
const T* x ,
const T& a2,
T* y)
{
Opm::BLAS_LAPACK::GEMV("Transpose",
A.numRows(), A.numCols(),
a1, A.data(), A.leadingDimension(),
x, 1, a2, y, 1);
}
template<typename T ,
template<typename> class SP>
void vecMulAdd_N(const T& a1,
const FullMatrix<T,SP,COrdering>& A ,
const T* x ,
const T& a2,
T* y)
{
typedef FullMatrix<T, ImmutableSharedData, FortranOrdering> FMAT;
const FMAT At(A.numCols(), A.numRows(), A.data());
vecMulAdd_T(a1, At, x, a2, y);
}
template<typename T ,
template<typename> class SP1,
template<typename> class SP2,
template<typename> class SP3>
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)
{
assert(A.numRows() == C.numRows());
assert(A.numCols() == B.numRows());
assert(B.numCols() == C.numCols());
int m = A.numRows(); // Number of *rows* in A
int n = B.numCols(); // Number of *cols* in B
int k = A.numCols(); // Number of *cols* in A (== numer of *rows* in B)
Opm::BLAS_LAPACK::GEMM("No Transpose", "No Transpose", m, n, k,
a1, A.data(), A.leadingDimension(),
B.data(), B.leadingDimension(),
a2, C.data(), C.leadingDimension());
}
template<typename T ,
template<typename> class SP1,
template<typename> class SP2,
template<typename> class SP3>
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)
{
assert(A.numRows() == C.numRows());
assert(B.numRows() == C.numCols());
assert(A.numCols() == B.numCols());
int m = A.numRows(); // Number of *rows* in A
int n = B.numRows(); // Number of *cols* in B'
int k = A.numCols(); // Number of *cols* in A (== numer of *rows* in B')
Opm::BLAS_LAPACK::GEMM("No Transpose", "Transpose", m, n, k,
a1, A.data(), A.leadingDimension(),
B.data(), B.leadingDimension(),
a2, C.data(), C.leadingDimension());
}
template<typename T ,
template<typename> class SP1,
template<typename> class SP2,
template<typename> class SP3>
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)
{
assert (A.numCols() == C.numRows());
assert (A.numRows() == B.numRows());
assert (B.numCols() == C.numCols());
int m = A.numCols(); // Number of *rows* in A'
int n = B.numCols(); // Number of *cols* in B
int k = A.numRows(); // Number of *cols* in A' (== numer of *rows* in B)
Opm::BLAS_LAPACK::GEMM("Transpose", "No Transpose", m, n, k,
a1, A.data(), A.leadingDimension(),
B.data(), B.leadingDimension(),
a2, C.data(), C.leadingDimension());
}
template<typename T ,
template<typename> class SP1,
template<typename> class SP2,
template<typename> class SP3>
void matMulAdd_NN(const T& a1,
const FullMatrix<T,SP1,FortranOrdering>& A ,
const FullMatrix<T,SP2,COrdering>& B ,
const T& a2,
FullMatrix<T,SP3,FortranOrdering>& C)
{
typedef typename FullMatrix<T,SP2,COrdering>::value_type value_type;
typedef FullMatrix<value_type,ImmutableSharedData,FortranOrdering> FMat;
const FMat Bt(B.numCols(), B.numRows(), B.data());
matMulAdd_NT(a1, A, Bt, a2, C);
}
template<class charT, class traits,
typename T, template<typename> class SP, class OP>
std::basic_ostream<charT,traits>&
operator<<(std::basic_ostream<charT,traits>& os,
const FullMatrix<T,SP,OP>& A)
{
for (int i = 0; i < A.numRows(); ++i) {
for (int j = 0; j < A.numCols(); ++j)
os << A(i,j) << ' ';
os << '\n';
}
return os;
}
} // namespace Opm
#endif // OPENRS_MATRIX_HEADER
Opm::ImmutableSharedData::data
const T * data() const
Direct access to all data.
Definition Matrix.hpp:191
Opm::ImmutableSharedData::ImmutableSharedData
ImmutableSharedData(int sz, const T *data)
Constructor.
Definition Matrix.hpp:202
Opm::ImmutableSharedData::operator[]
const T & operator[](int i) const
Storage element access.
Definition Matrix.hpp:181
Opm::ImmutableSharedData::size
int size() const
Data size query.
Definition Matrix.hpp:186
Opm::OwnData::OwnData
OwnData(int sz, const T *data)
Constructor.
Definition Matrix.hpp:99
Opm::OwnData::operator[]
T & operator[](int i)
Storage element access.
Definition Matrix.hpp:72
Opm::OwnData::data
T * data()
Direct access to all data.
Definition Matrix.hpp:84
Opm::OwnData::size
int size() const
Data size query.
Definition Matrix.hpp:79
Opm::SharedData::operator[]
T & operator[](int i)
Storage element access.
Definition Matrix.hpp:129
Opm::SharedData::size
int size() const
Data size query.
Definition Matrix.hpp:135
Opm::SharedData::SharedData
SharedData(int sz, T *data)
Constructor.
Definition Matrix.hpp:153
Opm::SharedData::data
T * data()
Direct access to all data.
Definition Matrix.hpp:140
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::vecMulAdd_T
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
Opm::matMulAdd_NN
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
Opm::orthogonalizeColumns
int orthogonalizeColumns(FullMatrix< T, StoragePolicy, FortranOrdering > &A)
Construct orthonormal basis for matrix range (i.e., column space).
Definition Matrix.hpp:729
Opm::symmetricUpdate
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
Opm::trace
Matrix::value_type trace(const Matrix &A)
Compute matrix trace (i.e., sum(diag(A))).
Definition Matrix.hpp:638
Opm::matMulAdd_NT
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
Opm::zero
void zero(Matrix &A)
Zero-fill a.
Definition Matrix.hpp:602
Opm::OwnFortranMatrix
FullMatrix< double, OwnData, FortranOrdering > OwnFortranMatrix
Convenience typedefs for Fortran-ordered.
Definition Matrix.hpp:588
Opm::invert
int invert(FullMatrix< T, StoragePolicy, OrderingPolicy > &A)
Matrix inversion, .
Definition Matrix.hpp:780
Opm::operator<<
std::basic_ostream< charT, traits > & operator<<(std::basic_ostream< charT, traits > &os, const BCBase< T > &bc)
Stream insertion for BCBase.
Definition BoundaryConditions.hpp:110
Opm::vecMulAdd_N
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
Opm::matMulAdd_TN
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
Opm::eye
void eye(Matrix &A)
Make an identity.
Definition Matrix.hpp:617
Opm::OwnCMatrix
FullMatrix< double, OwnData, COrdering > OwnCMatrix
Convenience typedefs for C-ordered.
Definition Matrix.hpp:579
Opm::prod
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