575855f896
this class can do both, more and less than dune-istl's BTDMatrix. more because it supports entries on the lower left and upper right, less because it does not support a block structure. The primary motivation for this class were the spline classes which for which the former feature is required to implement periodic splines and the latter is not necessary...
850 lines
22 KiB
C++
850 lines
22 KiB
C++
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
|
|
// vi: set et ts=4 sw=4 sts=4:
|
|
/*****************************************************************************
|
|
* Copyright (C) 2010-2013 by Andreas Lauser *
|
|
* *
|
|
* This program 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. *
|
|
* *
|
|
* This program is distributed in the hope that it will be useful, *
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
* MERCHANTBILITY 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 this program. If not, see <http://www.gnu.org/licenses/>. *
|
|
*****************************************************************************/
|
|
/*!
|
|
* \file
|
|
* \copydetails Opm::TridiagonalMatrix
|
|
*/
|
|
#ifndef OPM_TRIDIAGONAL_MATRIX_HH
|
|
#define OPM_TRIDIAGONAL_MATRIX_HH
|
|
|
|
#include <iostream>
|
|
#include <vector>
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
|
|
#include <assert.h>
|
|
|
|
namespace Opm {
|
|
|
|
/*!
|
|
* \brief Provides a tridiagonal matrix that also supports non-zero
|
|
* entries in the upper right and lower left
|
|
*
|
|
* The entries in the lower left and upper right are supported to make
|
|
* implementing periodic systems easy.
|
|
*
|
|
* The API of this class is designed to be close to the one used by
|
|
* the DUNE matrix classes.
|
|
*/
|
|
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); }
|
|
|
|
/*!
|
|
* \brief Prefix increment operator
|
|
*/
|
|
TridiagRow_ &operator++()
|
|
{ ++ rowIdx_; return *this; }
|
|
|
|
/*!
|
|
* \brief Prefix decrement operator
|
|
*/
|
|
TridiagRow_ &operator--()
|
|
{ -- rowIdx_; return *this; }
|
|
|
|
/*!
|
|
* \brief Comparision operator
|
|
*/
|
|
bool operator==(const TridiagRow_ &other) const
|
|
{ return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
|
|
|
|
/*!
|
|
* \brief Comparision operator
|
|
*/
|
|
bool operator!=(const TridiagRow_ &other) const
|
|
{ return !operator==(other); }
|
|
|
|
/*!
|
|
* \brief Dereference operator
|
|
*/
|
|
TridiagRow_ &operator*()
|
|
{ return *this; }
|
|
|
|
/*!
|
|
* \brief Return the row index of the this row.
|
|
*
|
|
* 0 is the first row.
|
|
*/
|
|
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(int numRows = 0)
|
|
{
|
|
resize(numRows);
|
|
};
|
|
|
|
TridiagonalMatrix(int numRows, Scalar value)
|
|
{
|
|
resize(numRows);
|
|
this->operator=(value);
|
|
};
|
|
|
|
/*!
|
|
* \brief Copy constructor.
|
|
*/
|
|
TridiagonalMatrix(const TridiagonalMatrix &source)
|
|
{ *this = source; };
|
|
|
|
/*!
|
|
* \brief Return the number of rows/columns of the matrix.
|
|
*/
|
|
size_t size() const
|
|
{ return diag_[0].size(); }
|
|
|
|
/*!
|
|
* \brief Return the number of rows of the matrix.
|
|
*/
|
|
size_t rows() const
|
|
{ return size(); }
|
|
|
|
/*!
|
|
* \brief Return the number of columns of the matrix.
|
|
*/
|
|
size_t cols() const
|
|
{ return size(); }
|
|
|
|
/*!
|
|
* \brief Change the number of rows of the matrix.
|
|
*/
|
|
void resize(size_t n)
|
|
{
|
|
if (n == size())
|
|
return;
|
|
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
|
|
diag_[diagIdx].resize(n);
|
|
}
|
|
|
|
/*!
|
|
* \brief Access an entry.
|
|
*/
|
|
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];
|
|
}
|
|
|
|
int diagIdx = 1 + colIdx - rowIdx;
|
|
// make sure that the requested column is in range
|
|
assert(0 <= diagIdx && diagIdx < 3);
|
|
return diag_[diagIdx][colIdx];
|
|
}
|
|
|
|
/*!
|
|
* \brief Access an entry.
|
|
*/
|
|
Scalar at(size_t rowIdx, size_t colIdx) const
|
|
{
|
|
int 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];
|
|
|
|
int diagIdx = 1 + colIdx - rowIdx;
|
|
// make sure that the requested column is in range
|
|
assert(0 <= diagIdx && diagIdx < 3);
|
|
return diag_[diagIdx][colIdx];
|
|
}
|
|
|
|
/*!
|
|
* \brief Assignment operator from another tridiagonal matrix.
|
|
*/
|
|
TridiagonalMatrix &operator=(const TridiagonalMatrix &source)
|
|
{
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
|
|
diag_[diagIdx] = source.diag_[diagIdx];
|
|
|
|
return *this;
|
|
}
|
|
|
|
/*!
|
|
* \brief Assignment operator from a Scalar.
|
|
*/
|
|
TridiagonalMatrix &operator=(Scalar value)
|
|
{
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
|
|
diag_[diagIdx].assign(size(), value);
|
|
|
|
return *this;
|
|
}
|
|
|
|
/*!
|
|
* \begin Iterator for the first row
|
|
*/
|
|
iterator begin()
|
|
{ return TridiagRow_(*this, 0); }
|
|
|
|
/*!
|
|
* \begin Const iterator for the first row
|
|
*/
|
|
const_iterator begin() const
|
|
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), 0); }
|
|
|
|
/*!
|
|
* \begin Const iterator for the next-to-last row
|
|
*/
|
|
const_iterator end() const
|
|
{ return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), size()); }
|
|
|
|
/*!
|
|
* \brief Row access operator.
|
|
*/
|
|
TridiagRow_ operator[](size_t rowIdx)
|
|
{ return TridiagRow_(*this, rowIdx); }
|
|
|
|
/*!
|
|
* \brief Row access operator.
|
|
*/
|
|
const TridiagRow_ operator[](size_t rowIdx) const
|
|
{ return TridiagRow_(*this, rowIdx); }
|
|
|
|
/*!
|
|
* \brief Multiplication with a Scalar
|
|
*/
|
|
TridiagonalMatrix &operator*=(Scalar alpha)
|
|
{
|
|
int n = size();
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
|
|
for (int i = 0; i < n; ++i) {
|
|
diag_[diagIdx][i] *= alpha;
|
|
}
|
|
}
|
|
|
|
return *this;
|
|
}
|
|
|
|
/*!
|
|
* \brief Division by a Scalar
|
|
*/
|
|
TridiagonalMatrix &operator/=(Scalar alpha)
|
|
{
|
|
alpha = 1.0/alpha;
|
|
int n = size();
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) {
|
|
for (int i = 0; i < n; ++i) {
|
|
diag_[diagIdx][i] *= alpha;
|
|
}
|
|
}
|
|
|
|
return *this;
|
|
}
|
|
|
|
/*!
|
|
* \brief Subtraction operator
|
|
*/
|
|
TridiagonalMatrix &operator-=(const TridiagonalMatrix &other)
|
|
{ return axpy(-1.0, other); }
|
|
|
|
/*!
|
|
* \brief Addition operator
|
|
*/
|
|
TridiagonalMatrix &operator+=(const TridiagonalMatrix &other)
|
|
{ return axpy(1.0, other); }
|
|
|
|
|
|
/*!
|
|
* \brief Multiply and add the matrix entries of another
|
|
* tridiagonal matrix.
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.axpy(alpha, B)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* A += alpha*C
|
|
* \endcode
|
|
*/
|
|
TridiagonalMatrix &axpy(Scalar alpha, const TridiagonalMatrix &other)
|
|
{
|
|
assert(size() == other.size());
|
|
|
|
int n = size();
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
|
|
for (int i = 0; i < n; ++ i)
|
|
diag_[diagIdx][i] += alpha * other[diagIdx][i];
|
|
|
|
return *this;
|
|
}
|
|
|
|
/*!
|
|
* \brief Matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.mv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y = A*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Additive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.umv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y += A*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Subtractive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.mmv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y -= A*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Scaled additive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.usmv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y += alpha*(A*x)
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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]);
|
|
}
|
|
|
|
/*!
|
|
* \brief Transposed matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.mtv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y = A^T*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Transposed additive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.umtv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y += A^T*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Transposed subtractive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.mmtv(x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y -= A^T*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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];
|
|
}
|
|
|
|
/*!
|
|
* \brief Transposed scaled additive matrix-vector product
|
|
*
|
|
* This means that
|
|
* \code
|
|
* A.umtv(alpha, x, y)
|
|
* \endcode
|
|
* is equivalent to
|
|
* \code
|
|
* y += alpha*A^T*x
|
|
* \endcode
|
|
*/
|
|
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
|
|
int n = size();
|
|
for (int 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]);
|
|
}
|
|
|
|
/*!
|
|
* \brief Calculate the frobenius norm
|
|
*
|
|
* i.e., the square root of the sum of all squared entries. This
|
|
* corresponds to the euclidean norm for vectors.
|
|
*/
|
|
Scalar frobeniusNorm() const
|
|
{ return std::sqrt(frobeniusNormSquared()); }
|
|
|
|
/*!
|
|
* \brief Calculate the squared frobenius norm
|
|
*
|
|
* i.e., the sum of all squared entries.
|
|
*/
|
|
Scalar frobeniusNormSquared() const
|
|
{
|
|
Scalar result = 0;
|
|
|
|
int n = size();
|
|
for (int i = 0; i < n; ++ i)
|
|
for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
|
|
result += diag_[diagIdx][i];
|
|
|
|
return result;
|
|
}
|
|
|
|
/*!
|
|
* \brief Calculate the infinity norm
|
|
*
|
|
* i.e., the maximum of the sum of the absolute values of all rows.
|
|
*/
|
|
Scalar infinityNorm() const
|
|
{
|
|
Scalar result=0;
|
|
|
|
// deal with rows 1 .. n-2
|
|
int n = size();
|
|
for (int 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;
|
|
}
|
|
|
|
/*!
|
|
* \brief Calculate the solution for a linear system of equations
|
|
*
|
|
* i.e., calculate x, so that it solves Ax = b, where A is a
|
|
* tridiagonal matrix.
|
|
*/
|
|
template <class XVector, class BVector>
|
|
void solve(XVector &x, const BVector &b) const
|
|
{
|
|
if (size() > 2 && diag_[2][0] != 0)
|
|
solveWithUpperRight_(x, b);
|
|
else
|
|
solveWithoutUpperRight_(x, b);
|
|
}
|
|
|
|
/*!
|
|
* \brief Print the matrix to a given output stream.
|
|
*/
|
|
void print(std::ostream &os = std::cout) const
|
|
{
|
|
int n = size();
|
|
|
|
// row 0
|
|
os << at(0, 0) << "\t"
|
|
<< at(0, 1) << "\t";
|
|
|
|
if (n > 3)
|
|
os << "\t";
|
|
if (n > 2)
|
|
os << at(0, n-1);
|
|
os << "\n";
|
|
|
|
// row 1 .. n - 2
|
|
for (int rowIdx = 1; rowIdx < n-1; ++rowIdx) {
|
|
if (rowIdx > 1)
|
|
os << "\t";
|
|
if (rowIdx == n - 2)
|
|
os << "\t";
|
|
|
|
os << at(rowIdx, rowIdx - 1) << "\t"
|
|
<< at(rowIdx, rowIdx) << "\t"
|
|
<< at(rowIdx, rowIdx + 1) << "\n";
|
|
}
|
|
|
|
// row n - 1
|
|
if (n > 2)
|
|
os << at(n-1, 0) << "\t";
|
|
if (n > 3)
|
|
os << "\t";
|
|
if (n > 4)
|
|
os << "\t";
|
|
os << at(n-1, n-2) << "\t"
|
|
<< at(n-1, n-1) << "\n";
|
|
}
|
|
|
|
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 = n - 2; i >= 0; --i) {
|
|
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1] - x[n-1]*lastColumn[i])/mainDiag[i];
|
|
}
|
|
}
|
|
|
|
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 = n - 2; i >= 0; --i) {
|
|
x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1])/mainDiag[i];
|
|
}
|
|
}
|
|
|
|
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
|