27#ifndef OPM_TRIDIAGONAL_MATRIX_HH
28#define OPM_TRIDIAGONAL_MATRIX_HH
48template <
class Scalar>
57 Scalar& operator[](
size_t colIdx)
58 {
return matrix_.at(rowIdx_, colIdx); }
60 Scalar operator[](
size_t colIdx)
const
61 {
return matrix_.at(rowIdx_, colIdx); }
66 TridiagRow_& operator++()
67 { ++ rowIdx_;
return *
this; }
72 TridiagRow_& operator--()
73 { -- rowIdx_;
return *
this; }
78 bool operator==(
const TridiagRow_& other)
const
79 {
return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
84 bool operator!=(
const TridiagRow_& other)
const
85 {
return !operator==(other); }
90 TridiagRow_& operator*()
103 mutable size_t rowIdx_;
107 typedef Scalar FieldType;
108 typedef TridiagRow_ RowType;
109 typedef size_t SizeType;
110 typedef TridiagRow_ iterator;
111 typedef TridiagRow_ const_iterator;
134 {
return diag_[0].size(); }
156 for (
int diagIdx = 0; diagIdx < 3; ++ diagIdx)
163 Scalar&
at(
size_t rowIdx,
size_t colIdx)
169 if (rowIdx == 0 && colIdx == n - 1)
171 if (rowIdx == n - 1 && colIdx == 0)
172 return diag_[0][n - 1];
175 size_t diagIdx = 1 + colIdx - rowIdx;
178 return diag_[diagIdx][colIdx];
184 Scalar
at(
size_t rowIdx,
size_t colIdx)
const
189 if (rowIdx == 0 && colIdx == n - 1)
191 if (rowIdx == n - 1 && colIdx == 0)
192 return diag_[0][n - 1];
194 size_t diagIdx = 1 + colIdx - rowIdx;
197 return diag_[diagIdx][colIdx];
205 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
206 diag_[diagIdx] = source.diag_[diagIdx];
216 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
217 diag_[diagIdx].assign(
size(), value);
226 {
return TridiagRow_(*
this, 0); }
237 const_iterator
end()
const
244 {
return TridiagRow_(*
this, rowIdx); }
250 {
return TridiagRow_(*
this, rowIdx); }
258 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
259 for (
unsigned i = 0; i < n; ++i) {
260 diag_[diagIdx][i] *= alpha;
274 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
275 for (
unsigned i = 0; i < n; ++i) {
276 diag_[diagIdx][i] *= alpha;
287 {
return axpy(-1.0, other); }
293 {
return axpy(1.0, other); }
314 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
315 for (
unsigned i = 0; i < n; ++ i)
316 diag_[diagIdx][i] += alpha * other[diagIdx][i];
333 template<
class Vector>
334 void mv(
const Vector& source, Vector& dest)
const
336 assert(source.size() ==
size());
337 assert(dest.size() ==
size());
342 for (
unsigned i = 1; i < n - 1; ++ i) {
344 diag_[0][i - 1]*source[i-1] +
345 diag_[1][i]*source[i] +
346 diag_[2][i + 1]*source[i + 1];
351 diag_[1][0]*source[0] +
352 diag_[2][1]*source[1] +
353 diag_[2][0]*source[n - 1];
356 diag_[0][n-1]*source[0] +
357 diag_[0][n-2]*source[n-2] +
358 diag_[1][n-1]*source[n-1];
373 template<
class Vector>
374 void umv(
const Vector& source, Vector& dest)
const
376 assert(source.size() ==
size());
377 assert(dest.size() ==
size());
382 for (
unsigned i = 1; i < n - 1; ++ i) {
384 diag_[0][i - 1]*source[i-1] +
385 diag_[1][i]*source[i] +
386 diag_[2][i + 1]*source[i + 1];
391 diag_[1][0]*source[0] +
392 diag_[2][1]*source[1] +
393 diag_[2][0]*source[n - 1];
396 diag_[0][n-1]*source[0] +
397 diag_[0][n-2]*source[n-2] +
398 diag_[1][n-1]*source[n-1];
413 template<
class Vector>
414 void mmv(
const Vector& source, Vector& dest)
const
416 assert(source.size() ==
size());
417 assert(dest.size() ==
size());
422 for (
unsigned i = 1; i < n - 1; ++ i) {
424 diag_[0][i - 1]*source[i-1] +
425 diag_[1][i]*source[i] +
426 diag_[2][i + 1]*source[i + 1];
431 diag_[1][0]*source[0] +
432 diag_[2][1]*source[1] +
433 diag_[2][0]*source[n - 1];
436 diag_[0][n-1]*source[0] +
437 diag_[0][n-2]*source[n-2] +
438 diag_[1][n-1]*source[n-1];
453 template<
class Vector>
454 void usmv(Scalar alpha,
const Vector& source, Vector& dest)
const
456 assert(source.size() ==
size());
457 assert(dest.size() ==
size());
462 for (
unsigned i = 1; i < n - 1; ++ i) {
465 diag_[0][i - 1]*source[i-1] +
466 diag_[1][i]*source[i] +
467 diag_[2][i + 1]*source[i + 1]);
473 diag_[1][0]*source[0] +
474 diag_[2][1]*source[1] +
475 diag_[2][0]*source[n - 1]);
479 diag_[0][n-1]*source[0] +
480 diag_[0][n-2]*source[n-2] +
481 diag_[1][n-1]*source[n-1]);
496 template<
class Vector>
497 void mtv(
const Vector& source, Vector& dest)
const
499 assert(source.size() ==
size());
500 assert(dest.size() ==
size());
505 for (
unsigned i = 1; i < n - 1; ++ i) {
507 diag_[2][i + 1]*source[i-1] +
508 diag_[1][i]*source[i] +
509 diag_[0][i - 1]*source[i + 1];
514 diag_[1][0]*source[0] +
515 diag_[0][1]*source[1] +
516 diag_[0][n-1]*source[n - 1];
519 diag_[2][0]*source[0] +
520 diag_[2][n-1]*source[n-2] +
521 diag_[1][n-1]*source[n-1];
536 template<
class Vector>
537 void umtv(
const Vector& source, Vector& dest)
const
539 assert(source.size() ==
size());
540 assert(dest.size() ==
size());
545 for (
unsigned i = 1; i < n - 1; ++ i) {
547 diag_[2][i + 1]*source[i-1] +
548 diag_[1][i]*source[i] +
549 diag_[0][i - 1]*source[i + 1];
554 diag_[1][0]*source[0] +
555 diag_[0][1]*source[1] +
556 diag_[0][n-1]*source[n - 1];
559 diag_[2][0]*source[0] +
560 diag_[2][n-1]*source[n-2] +
561 diag_[1][n-1]*source[n-1];
576 template<
class Vector>
577 void mmtv (
const Vector& source, Vector& dest)
const
579 assert(source.size() ==
size());
580 assert(dest.size() ==
size());
585 for (
unsigned i = 1; i < n - 1; ++ i) {
587 diag_[2][i + 1]*source[i-1] +
588 diag_[1][i]*source[i] +
589 diag_[0][i - 1]*source[i + 1];
594 diag_[1][0]*source[0] +
595 diag_[0][1]*source[1] +
596 diag_[0][n-1]*source[n - 1];
599 diag_[2][0]*source[0] +
600 diag_[2][n-1]*source[n-2] +
601 diag_[1][n-1]*source[n-1];
616 template<
class Vector>
617 void usmtv(Scalar alpha,
const Vector& source, Vector& dest)
const
619 assert(source.size() ==
size());
620 assert(dest.size() ==
size());
625 for (
unsigned i = 1; i < n - 1; ++ i) {
628 diag_[2][i + 1]*source[i-1] +
629 diag_[1][i]*source[i] +
630 diag_[0][i - 1]*source[i + 1]);
636 diag_[1][0]*source[0] +
637 diag_[0][1]*source[1] +
638 diag_[0][n-1]*source[n - 1]);
642 diag_[2][0]*source[0] +
643 diag_[2][n-1]*source[n-2] +
644 diag_[1][n-1]*source[n-1]);
666 for (
unsigned i = 0; i < n; ++ i)
667 for (
unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
668 result += diag_[diagIdx][i];
684 for (
unsigned i = 1; i < n - 1; ++ i) {
685 result = std::max(result,
686 std::abs(diag_[0][i - 1]) +
687 std::abs(diag_[1][i]) +
688 std::abs(diag_[2][i + 1]));
692 result = std::max(result,
693 std::abs(diag_[1][0]) +
694 std::abs(diag_[2][1]) +
695 std::abs(diag_[2][0]));
698 result = std::max(result,
699 std::abs(diag_[0][n-1]) +
700 std::abs(diag_[0][n-2]) +
701 std::abs(diag_[1][n-2]));
712 template <
class XVector,
class BVector>
713 void solve(XVector& x,
const BVector& b)
const
715 if (
size() > 2 && std::abs(diag_[2][0]) < 1e-30)
716 solveWithUpperRight_(x, b);
718 solveWithoutUpperRight_(x, b);
724 void print(std::ostream& os)
const;
727 template <
class XVector,
class BVector>
728 void solveWithUpperRight_(XVector& x,
const BVector& b)
const
732 std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n);
733 std::vector<Scalar> bStar(n);
734 std::copy(b.begin(), b.end(), bStar.begin());
736 lastColumn[0] = upperDiag[0];
739 for (
size_t i = 1; i < n; ++i) {
740 Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
742 lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
743 mainDiag[i] -= alpha * upperDiag[i];
745 bStar[i] -= alpha * bStar[i - 1];
749 if (lowerDiag[n - 1] != 0.0 &&
size() > 2) {
750 Scalar lastRow = lowerDiag[n - 1];
751 for (
size_t i = 0; i < n - 1; ++i) {
752 Scalar alpha = lastRow/mainDiag[i];
753 lastRow = - alpha*upperDiag[i + 1];
754 bStar[n - 1] -= alpha * bStar[i];
757 mainDiag[n-1] += lastRow;
761 x[n - 1] = bStar[n - 1]/mainDiag[n-1];
762 for (
int i =
static_cast<int>(n) - 2; i >= 0; --i) {
763 unsigned iu =
static_cast<unsigned>(i);
764 x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1] - x[n-1]*lastColumn[iu])/mainDiag[iu];
768 template <
class XVector,
class BVector>
769 void solveWithoutUpperRight_(XVector& x,
const BVector& b)
const
773 std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]);
774 std::vector<Scalar> bStar(n);
775 std::copy(b.begin(), b.end(), bStar.begin());
778 for (
size_t i = 1; i < n; ++i) {
779 Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
781 lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
782 mainDiag[i] -= alpha * upperDiag[i];
784 bStar[i] -= alpha * bStar[i - 1];
788 if (lowerDiag[n - 1] != 0.0 &&
size() > 2) {
789 Scalar lastRow = lowerDiag[n - 1];
790 for (
size_t i = 0; i < n - 1; ++i) {
791 Scalar alpha = lastRow/mainDiag[i];
792 lastRow = - alpha*upperDiag[i + 1];
793 bStar[n - 1] -= alpha * bStar[i];
796 mainDiag[n-1] += lastRow;
800 x[n - 1] = bStar[n - 1]/mainDiag[n-1];
801 for (
int i =
static_cast<int>(n) - 2; i >= 0; --i) {
802 unsigned iu =
static_cast<unsigned>(i);
803 x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1])/mainDiag[iu];
807 mutable std::vector<Scalar> diag_[3];
812template <
class Scalar>
Provides a tridiagonal matrix that also supports non-zero entries in the upper right and lower left.
Definition TridiagonalMatrix.hpp:50
const_iterator end() const
\begin Const iterator for the next-to-last row
Definition TridiagonalMatrix.hpp:237
TridiagonalMatrix(const TridiagonalMatrix &source)
Copy constructor.
Definition TridiagonalMatrix.hpp:127
TridiagRow_ operator[](size_t rowIdx)
Row access operator.
Definition TridiagonalMatrix.hpp:243
size_t cols() const
Return the number of columns of the matrix.
Definition TridiagonalMatrix.hpp:145
size_t rows() const
Return the number of rows of the matrix.
Definition TridiagonalMatrix.hpp:139
Scalar frobeniusNormSquared() const
Calculate the squared frobenius norm.
Definition TridiagonalMatrix.hpp:661
TridiagonalMatrix & operator/=(Scalar alpha)
Division by a Scalar.
Definition TridiagonalMatrix.hpp:270
TridiagonalMatrix & operator=(Scalar value)
Assignment operator from a Scalar.
Definition TridiagonalMatrix.hpp:214
void umtv(const Vector &source, Vector &dest) const
Transposed additive matrix-vector product.
Definition TridiagonalMatrix.hpp:537
const_iterator begin() const
\begin Const iterator for the first row
Definition TridiagonalMatrix.hpp:231
iterator begin()
\begin Iterator for the first row
Definition TridiagonalMatrix.hpp:225
void usmtv(Scalar alpha, const Vector &source, Vector &dest) const
Transposed scaled additive matrix-vector product.
Definition TridiagonalMatrix.hpp:617
TridiagonalMatrix & operator+=(const TridiagonalMatrix &other)
Addition operator.
Definition TridiagonalMatrix.hpp:292
void mv(const Vector &source, Vector &dest) const
Matrix-vector product.
Definition TridiagonalMatrix.hpp:334
Scalar infinityNorm() const
Calculate the infinity norm.
Definition TridiagonalMatrix.hpp:678
Scalar at(size_t rowIdx, size_t colIdx) const
Access an entry.
Definition TridiagonalMatrix.hpp:184
const TridiagRow_ operator[](size_t rowIdx) const
Row access operator.
Definition TridiagonalMatrix.hpp:249
void mmv(const Vector &source, Vector &dest) const
Subtractive matrix-vector product.
Definition TridiagonalMatrix.hpp:414
void usmv(Scalar alpha, const Vector &source, Vector &dest) const
Scaled additive matrix-vector product.
Definition TridiagonalMatrix.hpp:454
void mmtv(const Vector &source, Vector &dest) const
Transposed subtractive matrix-vector product.
Definition TridiagonalMatrix.hpp:577
void print(std::ostream &os) const
Print the matrix to a given output stream.
Definition TridiagonalMatrix.cpp:32
TridiagonalMatrix & operator-=(const TridiagonalMatrix &other)
Subtraction operator.
Definition TridiagonalMatrix.hpp:286
size_t size() const
Return the number of rows/columns of the matrix.
Definition TridiagonalMatrix.hpp:133
TridiagonalMatrix & axpy(Scalar alpha, const TridiagonalMatrix &other)
Multiply and add the matrix entries of another tridiagonal matrix.
Definition TridiagonalMatrix.hpp:309
TridiagonalMatrix & operator=(const TridiagonalMatrix &source)
Assignment operator from another tridiagonal matrix.
Definition TridiagonalMatrix.hpp:203
void solve(XVector &x, const BVector &b) const
Calculate the solution for a linear system of equations.
Definition TridiagonalMatrix.hpp:713
void resize(size_t n)
Change the number of rows of the matrix.
Definition TridiagonalMatrix.hpp:151
void mtv(const Vector &source, Vector &dest) const
Transposed matrix-vector product.
Definition TridiagonalMatrix.hpp:497
Scalar frobeniusNorm() const
Calculate the frobenius norm.
Definition TridiagonalMatrix.hpp:653
void umv(const Vector &source, Vector &dest) const
Additive matrix-vector product.
Definition TridiagonalMatrix.hpp:374
Scalar & at(size_t rowIdx, size_t colIdx)
Access an entry.
Definition TridiagonalMatrix.hpp:163
TridiagonalMatrix & operator*=(Scalar alpha)
Multiplication with a Scalar.
Definition TridiagonalMatrix.hpp:255
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30