Loading...
Searching...
No Matches
TridiagonalMatrix.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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 2 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_TRIDIAGONAL_MATRIX_HH
#define OPM_TRIDIAGONAL_MATRIX_HH
#include <algorithm>
#include <cassert>
#include <cmath>
#include <iosfwd>
#include <vector>
namespace Opm {
template <class Scalar>
class TridiagonalMatrix
{
struct TridiagRow_ {
TridiagRow_(TridiagonalMatrix& m, size_t rowIdx)
: matrix_(m)
, rowIdx_(rowIdx)
{}
Scalar& operator[](size_t colIdx)
{ return matrix_.at(rowIdx_, colIdx); }
Scalar operator[](size_t colIdx) const
{ return matrix_.at(rowIdx_, colIdx); }
TridiagRow_& operator++()
{ ++ rowIdx_; return *this; }
TridiagRow_& operator--()
{ -- rowIdx_; return *this; }
bool operator==(const TridiagRow_& other) const
{ return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
bool operator!=(const TridiagRow_& other) const
{ return !operator==(other); }
TridiagRow_& operator*()
{ return *this; }
size_t index() const
{ return rowIdx_; }
private:
TridiagonalMatrix& matrix_;
mutable size_t rowIdx_;
};
public:
typedef Scalar FieldType;
typedef TridiagRow_ RowType;
typedef size_t SizeType;
typedef TridiagRow_ iterator;
typedef TridiagRow_ const_iterator;
explicit TridiagonalMatrix(size_t numRows = 0)
{
resize(numRows);
}
TridiagonalMatrix(size_t numRows, Scalar value)
{
resize(numRows);
this->operator=(value);
}
TridiagonalMatrix(const TridiagonalMatrix& source)
{ *this = source; }
size_t size() const
{ return diag_[0].size(); }
size_t rows() const
{ return size(); }
size_t cols() const
{ return size(); }
void resize(size_t n)
{
if (n == size())
return;
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].resize(n);
}
Scalar& at(size_t rowIdx, size_t colIdx)
{
size_t n = size();
// special cases
if (n > 2) {
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
}
size_t diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(diagIdx < 3);
return diag_[diagIdx][colIdx];
}
Scalar at(size_t rowIdx, size_t colIdx) const
{
size_t n = size();
// special cases
if (rowIdx == 0 && colIdx == n - 1)
return diag_[2][0];
if (rowIdx == n - 1 && colIdx == 0)
return diag_[0][n - 1];
size_t diagIdx = 1 + colIdx - rowIdx;
// make sure that the requested column is in range
assert(diagIdx < 3);
return diag_[diagIdx][colIdx];
}
TridiagonalMatrix& operator=(const TridiagonalMatrix& source)
{
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx] = source.diag_[diagIdx];
return *this;
}
TridiagonalMatrix& operator=(Scalar value)
{
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
diag_[diagIdx].assign(size(), value);
return *this;
}
iterator begin()
{ return TridiagRow_(*this, 0); }
const_iterator begin() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), 0); }
const_iterator end() const
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), size()); }
TridiagRow_ operator[](size_t rowIdx)
{ return TridiagRow_(*this, rowIdx); }
const TridiagRow_ operator[](size_t rowIdx) const
{ return TridiagRow_(*this, rowIdx); }
TridiagonalMatrix& operator*=(Scalar alpha)
{
unsigned n = size();
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (unsigned i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
TridiagonalMatrix& operator/=(Scalar alpha)
{
alpha = 1.0/alpha;
unsigned n = size();
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
for (unsigned i = 0; i < n; ++i) {
diag_[diagIdx][i] *= alpha;
}
}
return *this;
}
TridiagonalMatrix& operator-=(const TridiagonalMatrix& other)
{ return axpy(-1.0, other); }
TridiagonalMatrix& operator+=(const TridiagonalMatrix& other)
{ return axpy(1.0, other); }
TridiagonalMatrix& axpy(Scalar alpha, const TridiagonalMatrix& other)
{
assert(size() == other.size());
unsigned n = size();
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
for (unsigned i = 0; i < n; ++ i)
diag_[diagIdx][i] += alpha * other[diagIdx][i];
return *this;
}
template<class Vector>
void mv(const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] =
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void umv(const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] +=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void mmv(const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1];
dest[n - 1] -=
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void usmv(Scalar alpha, const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[0][i - 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[2][i + 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[2][1]*source[1] +
diag_[2][0]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[0][n-1]*source[0] +
diag_[0][n-2]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
template<class Vector>
void mtv(const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] =
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] =
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] =
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void umtv(const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] +=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] +=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] +=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void mmtv (const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] -=
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1];
}
// rows 0 and n-1
dest[0] -=
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1];
dest[n - 1] -=
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1];
}
template<class Vector>
void usmtv(Scalar alpha, const Vector& source, Vector& dest) const
{
assert(source.size() == size());
assert(dest.size() == size());
assert(size() > 1);
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
dest[i] +=
alpha*(
diag_[2][i + 1]*source[i-1] +
diag_[1][i]*source[i] +
diag_[0][i - 1]*source[i + 1]);
}
// rows 0 and n-1
dest[0] +=
alpha*(
diag_[1][0]*source[0] +
diag_[0][1]*source[1] +
diag_[0][n-1]*source[n - 1]);
dest[n - 1] +=
alpha*(
diag_[2][0]*source[0] +
diag_[2][n-1]*source[n-2] +
diag_[1][n-1]*source[n-1]);
}
Scalar frobeniusNorm() const
{ return std::sqrt(frobeniusNormSquared()); }
Scalar frobeniusNormSquared() const
{
Scalar result = 0;
unsigned n = size();
for (unsigned i = 0; i < n; ++ i)
for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
result += diag_[diagIdx][i];
return result;
}
Scalar infinityNorm() const
{
Scalar result=0;
// deal with rows 1 .. n-2
unsigned n = size();
for (unsigned i = 1; i < n - 1; ++ i) {
result = std::max(result,
std::abs(diag_[0][i - 1]) +
std::abs(diag_[1][i]) +
std::abs(diag_[2][i + 1]));
}
// rows 0 and n-1
result = std::max(result,
std::abs(diag_[1][0]) +
std::abs(diag_[2][1]) +
std::abs(diag_[2][0]));
result = std::max(result,
std::abs(diag_[0][n-1]) +
std::abs(diag_[0][n-2]) +
std::abs(diag_[1][n-2]));
return result;
}
template <class XVector, class BVector>
void solve(XVector& x, const BVector& b) const
{
if (size() > 2 && std::abs(diag_[2][0]) < 1e-30)
solveWithUpperRight_(x, b);
else
solveWithoutUpperRight_(x, b);
}
void print(std::ostream& os) const;
private:
template <class XVector, class BVector>
void solveWithUpperRight_(XVector& x, const BVector& b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
lastColumn[0] = upperDiag[0];
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = static_cast<int>(n) - 2; i >= 0; --i) {
unsigned iu = static_cast<unsigned>(i);
x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1] - x[n-1]*lastColumn[iu])/mainDiag[iu];
}
}
template <class XVector, class BVector>
void solveWithoutUpperRight_(XVector& x, const BVector& b) const
{
size_t n = size();
std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]);
std::vector<Scalar> bStar(n);
std::copy(b.begin(), b.end(), bStar.begin());
// forward elimination
for (size_t i = 1; i < n; ++i) {
Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
mainDiag[i] -= alpha * upperDiag[i];
bStar[i] -= alpha * bStar[i - 1];
};
// deal with the last row if the entry on the lower left is not zero
if (lowerDiag[n - 1] != 0.0 && size() > 2) {
Scalar lastRow = lowerDiag[n - 1];
for (size_t i = 0; i < n - 1; ++i) {
Scalar alpha = lastRow/mainDiag[i];
lastRow = - alpha*upperDiag[i + 1];
bStar[n - 1] -= alpha * bStar[i];
}
mainDiag[n-1] += lastRow;
}
// backward elimination
x[n - 1] = bStar[n - 1]/mainDiag[n-1];
for (int i = static_cast<int>(n) - 2; i >= 0; --i) {
unsigned iu = static_cast<unsigned>(i);
x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1])/mainDiag[iu];
}
}
mutable std::vector<Scalar> diag_[3];
};
} // namespace Opm
template <class Scalar>
std::ostream& operator<<(std::ostream& os, const Opm::TridiagonalMatrix<Scalar>& mat)
{
mat.print(os);
return os;
}
#endif
Opm::TridiagonalMatrix
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition TridiagonalMatrix.hpp:50
Opm::TridiagonalMatrix::end
const_iterator end() const
\begin Const iterator for the next-to-last row
Definition TridiagonalMatrix.hpp:237
Opm::TridiagonalMatrix::TridiagonalMatrix
TridiagonalMatrix(const TridiagonalMatrix &source)
Copy constructor.
Definition TridiagonalMatrix.hpp:127
Opm::TridiagonalMatrix::operator[]
TridiagRow_ operator[](size_t rowIdx)
Row access operator.
Definition TridiagonalMatrix.hpp:243
Opm::TridiagonalMatrix::cols
size_t cols() const
Return the number of columns of the matrix.
Definition TridiagonalMatrix.hpp:145
Opm::TridiagonalMatrix::rows
size_t rows() const
Return the number of rows of the matrix.
Definition TridiagonalMatrix.hpp:139
Opm::TridiagonalMatrix::frobeniusNormSquared
Scalar frobeniusNormSquared() const
Calculate the squared frobenius norm.
Definition TridiagonalMatrix.hpp:661
Opm::TridiagonalMatrix::operator/=
TridiagonalMatrix & operator/=(Scalar alpha)
Division by a Scalar.
Definition TridiagonalMatrix.hpp:270
Opm::TridiagonalMatrix::operator=
TridiagonalMatrix & operator=(Scalar value)
Assignment operator from a Scalar.
Definition TridiagonalMatrix.hpp:214
Opm::TridiagonalMatrix::umtv
void umtv(const Vector &source, Vector &dest) const
Transposed additive matrix-vector product.
Definition TridiagonalMatrix.hpp:537
Opm::TridiagonalMatrix::begin
const_iterator begin() const
\begin Const iterator for the first row
Definition TridiagonalMatrix.hpp:231
Opm::TridiagonalMatrix::begin
iterator begin()
\begin Iterator for the first row
Definition TridiagonalMatrix.hpp:225
Opm::TridiagonalMatrix::usmtv
void usmtv(Scalar alpha, const Vector &source, Vector &dest) const
Transposed scaled additive matrix-vector product.
Definition TridiagonalMatrix.hpp:617
Opm::TridiagonalMatrix::operator+=
TridiagonalMatrix & operator+=(const TridiagonalMatrix &other)
Addition operator.
Definition TridiagonalMatrix.hpp:292
Opm::TridiagonalMatrix::mv
void mv(const Vector &source, Vector &dest) const
Matrix-vector product.
Definition TridiagonalMatrix.hpp:334
Opm::TridiagonalMatrix::infinityNorm
Scalar infinityNorm() const
Calculate the infinity norm.
Definition TridiagonalMatrix.hpp:678
Opm::TridiagonalMatrix::at
Scalar at(size_t rowIdx, size_t colIdx) const
Access an entry.
Definition TridiagonalMatrix.hpp:184
Opm::TridiagonalMatrix::operator[]
const TridiagRow_ operator[](size_t rowIdx) const
Row access operator.
Definition TridiagonalMatrix.hpp:249
Opm::TridiagonalMatrix::mmv
void mmv(const Vector &source, Vector &dest) const
Subtractive matrix-vector product.
Definition TridiagonalMatrix.hpp:414
Opm::TridiagonalMatrix::usmv
void usmv(Scalar alpha, const Vector &source, Vector &dest) const
Scaled additive matrix-vector product.
Definition TridiagonalMatrix.hpp:454
Opm::TridiagonalMatrix::mmtv
void mmtv(const Vector &source, Vector &dest) const
Transposed subtractive matrix-vector product.
Definition TridiagonalMatrix.hpp:577
Opm::TridiagonalMatrix::print
void print(std::ostream &os) const
Print the matrix to a given output stream.
Definition TridiagonalMatrix.cpp:32
Opm::TridiagonalMatrix::operator-=
TridiagonalMatrix & operator-=(const TridiagonalMatrix &other)
Subtraction operator.
Definition TridiagonalMatrix.hpp:286
Opm::TridiagonalMatrix::size
size_t size() const
Return the number of rows/columns of the matrix.
Definition TridiagonalMatrix.hpp:133
Opm::TridiagonalMatrix::axpy
TridiagonalMatrix & axpy(Scalar alpha, const TridiagonalMatrix &other)
Multiply and add the matrix entries of another tridiagonal matrix.
Definition TridiagonalMatrix.hpp:309
Opm::TridiagonalMatrix::operator=
TridiagonalMatrix & operator=(const TridiagonalMatrix &source)
Assignment operator from another tridiagonal matrix.
Definition TridiagonalMatrix.hpp:203
Opm::TridiagonalMatrix::solve
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition TridiagonalMatrix.hpp:713
Opm::TridiagonalMatrix::resize
void resize(size_t n)
Change the number of rows of the matrix.
Definition TridiagonalMatrix.hpp:151
Opm::TridiagonalMatrix::mtv
void mtv(const Vector &source, Vector &dest) const
Transposed matrix-vector product.
Definition TridiagonalMatrix.hpp:497
Opm::TridiagonalMatrix::frobeniusNorm
Scalar frobeniusNorm() const
Calculate the frobenius norm.
Definition TridiagonalMatrix.hpp:653
Opm::TridiagonalMatrix::umv
void umv(const Vector &source, Vector &dest) const
Additive matrix-vector product.
Definition TridiagonalMatrix.hpp:374
Opm::TridiagonalMatrix::at
Scalar & at(size_t rowIdx, size_t colIdx)
Access an entry.
Definition TridiagonalMatrix.hpp:163
Opm::TridiagonalMatrix::operator*=
TridiagonalMatrix & operator*=(Scalar alpha)
Multiplication with a Scalar.
Definition TridiagonalMatrix.hpp:255
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30