00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00027 #ifndef OPM_TRIDIAGONAL_MATRIX_HH
00028 #define OPM_TRIDIAGONAL_MATRIX_HH
00029
00030 #include <iostream>
00031 #include <vector>
00032 #include <algorithm>
00033 #include <cmath>
00034
00035 #include <assert.h>
00036
00037 namespace Opm {
00038
00049 template <class Scalar>
00050 class TridiagonalMatrix
00051 {
00052 struct TridiagRow_ {
00053 TridiagRow_(TridiagonalMatrix& m, size_t rowIdx)
00054 : matrix_(m)
00055 , rowIdx_(rowIdx)
00056 {}
00057
00058 Scalar& operator[](size_t colIdx)
00059 { return matrix_.at(rowIdx_, colIdx); }
00060
00061 Scalar operator[](size_t colIdx) const
00062 { return matrix_.at(rowIdx_, colIdx); }
00063
00067 TridiagRow_& operator++()
00068 { ++ rowIdx_; return *this; }
00069
00073 TridiagRow_& operator--()
00074 { -- rowIdx_; return *this; }
00075
00079 bool operator==(const TridiagRow_& other) const
00080 { return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; }
00081
00085 bool operator!=(const TridiagRow_& other) const
00086 { return !operator==(other); }
00087
00091 TridiagRow_& operator*()
00092 { return *this; }
00093
00099 size_t index() const
00100 { return rowIdx_; }
00101
00102 private:
00103 TridiagonalMatrix& matrix_;
00104 mutable size_t rowIdx_;
00105 };
00106
00107 public:
00108 typedef Scalar FieldType;
00109 typedef TridiagRow_ RowType;
00110 typedef size_t SizeType;
00111 typedef TridiagRow_ iterator;
00112 typedef TridiagRow_ const_iterator;
00113
00114 explicit TridiagonalMatrix(size_t numRows = 0)
00115 {
00116 resize(numRows);
00117 }
00118
00119 TridiagonalMatrix(size_t numRows, Scalar value)
00120 {
00121 resize(numRows);
00122 this->operator=(value);
00123 }
00124
00128 TridiagonalMatrix(const TridiagonalMatrix& source)
00129 { *this = source; }
00130
00134 size_t size() const
00135 { return diag_[0].size(); }
00136
00140 size_t rows() const
00141 { return size(); }
00142
00146 size_t cols() const
00147 { return size(); }
00148
00152 void resize(size_t n)
00153 {
00154 if (n == size())
00155 return;
00156
00157 for (int diagIdx = 0; diagIdx < 3; ++ diagIdx)
00158 diag_[diagIdx].resize(n);
00159 }
00160
00164 Scalar& at(size_t rowIdx, size_t colIdx)
00165 {
00166 size_t n = size();
00167
00168
00169 if (n > 2) {
00170 if (rowIdx == 0 && colIdx == n - 1)
00171 return diag_[2][0];
00172 if (rowIdx == n - 1 && colIdx == 0)
00173 return diag_[0][n - 1];
00174 }
00175
00176 size_t diagIdx = 1 + colIdx - rowIdx;
00177
00178 assert(0 <= diagIdx && diagIdx < 3);
00179 return diag_[diagIdx][colIdx];
00180 }
00181
00185 Scalar at(size_t rowIdx, size_t colIdx) const
00186 {
00187 size_t n = size();
00188
00189
00190 if (rowIdx == 0 && colIdx == n - 1)
00191 return diag_[2][0];
00192 if (rowIdx == n - 1 && colIdx == 0)
00193 return diag_[0][n - 1];
00194
00195 size_t diagIdx = 1 + colIdx - rowIdx;
00196
00197 assert(0 <= diagIdx && diagIdx < 3);
00198 return diag_[diagIdx][colIdx];
00199 }
00200
00204 TridiagonalMatrix& operator=(const TridiagonalMatrix& source)
00205 {
00206 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
00207 diag_[diagIdx] = source.diag_[diagIdx];
00208
00209 return *this;
00210 }
00211
00215 TridiagonalMatrix& operator=(Scalar value)
00216 {
00217 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
00218 diag_[diagIdx].assign(size(), value);
00219
00220 return *this;
00221 }
00222
00226 iterator begin()
00227 { return TridiagRow_(*this, 0); }
00228
00232 const_iterator begin() const
00233 { return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), 0); }
00234
00238 const_iterator end() const
00239 { return TridiagRow_(const_cast<TridiagonalMatrix&>(*this), size()); }
00240
00244 TridiagRow_ operator[](size_t rowIdx)
00245 { return TridiagRow_(*this, rowIdx); }
00246
00250 const TridiagRow_ operator[](size_t rowIdx) const
00251 { return TridiagRow_(*this, rowIdx); }
00252
00256 TridiagonalMatrix& operator*=(Scalar alpha)
00257 {
00258 unsigned n = size();
00259 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
00260 for (unsigned i = 0; i < n; ++i) {
00261 diag_[diagIdx][i] *= alpha;
00262 }
00263 }
00264
00265 return *this;
00266 }
00267
00271 TridiagonalMatrix& operator/=(Scalar alpha)
00272 {
00273 alpha = 1.0/alpha;
00274 unsigned n = size();
00275 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx) {
00276 for (unsigned i = 0; i < n; ++i) {
00277 diag_[diagIdx][i] *= alpha;
00278 }
00279 }
00280
00281 return *this;
00282 }
00283
00287 TridiagonalMatrix& operator-=(const TridiagonalMatrix& other)
00288 { return axpy(-1.0, other); }
00289
00293 TridiagonalMatrix& operator+=(const TridiagonalMatrix& other)
00294 { return axpy(1.0, other); }
00295
00296
00310 TridiagonalMatrix& axpy(Scalar alpha, const TridiagonalMatrix& other)
00311 {
00312 assert(size() == other.size());
00313
00314 unsigned n = size();
00315 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
00316 for (unsigned i = 0; i < n; ++ i)
00317 diag_[diagIdx][i] += alpha * other[diagIdx][i];
00318
00319 return *this;
00320 }
00321
00334 template<class Vector>
00335 void mv(const Vector& source, Vector& dest) const
00336 {
00337 assert(source.size() == size());
00338 assert(dest.size() == size());
00339 assert(size() > 1);
00340
00341
00342 unsigned n = size();
00343 for (unsigned i = 1; i < n - 1; ++ i) {
00344 dest[i] =
00345 diag_[0][i - 1]*source[i-1] +
00346 diag_[1][i]*source[i] +
00347 diag_[2][i + 1]*source[i + 1];
00348 }
00349
00350
00351 dest[0] =
00352 diag_[1][0]*source[0] +
00353 diag_[2][1]*source[1] +
00354 diag_[2][0]*source[n - 1];
00355
00356 dest[n - 1] =
00357 diag_[0][n-1]*source[0] +
00358 diag_[0][n-2]*source[n-2] +
00359 diag_[1][n-1]*source[n-1];
00360 }
00361
00374 template<class Vector>
00375 void umv(const Vector& source, Vector& dest) const
00376 {
00377 assert(source.size() == size());
00378 assert(dest.size() == size());
00379 assert(size() > 1);
00380
00381
00382 unsigned n = size();
00383 for (unsigned i = 1; i < n - 1; ++ i) {
00384 dest[i] +=
00385 diag_[0][i - 1]*source[i-1] +
00386 diag_[1][i]*source[i] +
00387 diag_[2][i + 1]*source[i + 1];
00388 }
00389
00390
00391 dest[0] +=
00392 diag_[1][0]*source[0] +
00393 diag_[2][1]*source[1] +
00394 diag_[2][0]*source[n - 1];
00395
00396 dest[n - 1] +=
00397 diag_[0][n-1]*source[0] +
00398 diag_[0][n-2]*source[n-2] +
00399 diag_[1][n-1]*source[n-1];
00400 }
00401
00414 template<class Vector>
00415 void mmv(const Vector& source, Vector& dest) const
00416 {
00417 assert(source.size() == size());
00418 assert(dest.size() == size());
00419 assert(size() > 1);
00420
00421
00422 unsigned n = size();
00423 for (unsigned i = 1; i < n - 1; ++ i) {
00424 dest[i] -=
00425 diag_[0][i - 1]*source[i-1] +
00426 diag_[1][i]*source[i] +
00427 diag_[2][i + 1]*source[i + 1];
00428 }
00429
00430
00431 dest[0] -=
00432 diag_[1][0]*source[0] +
00433 diag_[2][1]*source[1] +
00434 diag_[2][0]*source[n - 1];
00435
00436 dest[n - 1] -=
00437 diag_[0][n-1]*source[0] +
00438 diag_[0][n-2]*source[n-2] +
00439 diag_[1][n-1]*source[n-1];
00440 }
00441
00454 template<class Vector>
00455 void usmv(Scalar alpha, const Vector& source, Vector& dest) const
00456 {
00457 assert(source.size() == size());
00458 assert(dest.size() == size());
00459 assert(size() > 1);
00460
00461
00462 unsigned n = size();
00463 for (unsigned i = 1; i < n - 1; ++ i) {
00464 dest[i] +=
00465 alpha*(
00466 diag_[0][i - 1]*source[i-1] +
00467 diag_[1][i]*source[i] +
00468 diag_[2][i + 1]*source[i + 1]);
00469 }
00470
00471
00472 dest[0] +=
00473 alpha*(
00474 diag_[1][0]*source[0] +
00475 diag_[2][1]*source[1] +
00476 diag_[2][0]*source[n - 1]);
00477
00478 dest[n - 1] +=
00479 alpha*(
00480 diag_[0][n-1]*source[0] +
00481 diag_[0][n-2]*source[n-2] +
00482 diag_[1][n-1]*source[n-1]);
00483 }
00484
00497 template<class Vector>
00498 void mtv(const Vector& source, Vector& dest) const
00499 {
00500 assert(source.size() == size());
00501 assert(dest.size() == size());
00502 assert(size() > 1);
00503
00504
00505 unsigned n = size();
00506 for (unsigned i = 1; i < n - 1; ++ i) {
00507 dest[i] =
00508 diag_[2][i + 1]*source[i-1] +
00509 diag_[1][i]*source[i] +
00510 diag_[0][i - 1]*source[i + 1];
00511 }
00512
00513
00514 dest[0] =
00515 diag_[1][0]*source[0] +
00516 diag_[0][1]*source[1] +
00517 diag_[0][n-1]*source[n - 1];
00518
00519 dest[n - 1] =
00520 diag_[2][0]*source[0] +
00521 diag_[2][n-1]*source[n-2] +
00522 diag_[1][n-1]*source[n-1];
00523 }
00524
00537 template<class Vector>
00538 void umtv(const Vector& source, Vector& dest) const
00539 {
00540 assert(source.size() == size());
00541 assert(dest.size() == size());
00542 assert(size() > 1);
00543
00544
00545 unsigned n = size();
00546 for (unsigned i = 1; i < n - 1; ++ i) {
00547 dest[i] +=
00548 diag_[2][i + 1]*source[i-1] +
00549 diag_[1][i]*source[i] +
00550 diag_[0][i - 1]*source[i + 1];
00551 }
00552
00553
00554 dest[0] +=
00555 diag_[1][0]*source[0] +
00556 diag_[0][1]*source[1] +
00557 diag_[0][n-1]*source[n - 1];
00558
00559 dest[n - 1] +=
00560 diag_[2][0]*source[0] +
00561 diag_[2][n-1]*source[n-2] +
00562 diag_[1][n-1]*source[n-1];
00563 }
00564
00577 template<class Vector>
00578 void mmtv (const Vector& source, Vector& dest) const
00579 {
00580 assert(source.size() == size());
00581 assert(dest.size() == size());
00582 assert(size() > 1);
00583
00584
00585 unsigned n = size();
00586 for (unsigned i = 1; i < n - 1; ++ i) {
00587 dest[i] -=
00588 diag_[2][i + 1]*source[i-1] +
00589 diag_[1][i]*source[i] +
00590 diag_[0][i - 1]*source[i + 1];
00591 }
00592
00593
00594 dest[0] -=
00595 diag_[1][0]*source[0] +
00596 diag_[0][1]*source[1] +
00597 diag_[0][n-1]*source[n - 1];
00598
00599 dest[n - 1] -=
00600 diag_[2][0]*source[0] +
00601 diag_[2][n-1]*source[n-2] +
00602 diag_[1][n-1]*source[n-1];
00603 }
00604
00617 template<class Vector>
00618 void usmtv(Scalar alpha, const Vector& source, Vector& dest) const
00619 {
00620 assert(source.size() == size());
00621 assert(dest.size() == size());
00622 assert(size() > 1);
00623
00624
00625 unsigned n = size();
00626 for (unsigned i = 1; i < n - 1; ++ i) {
00627 dest[i] +=
00628 alpha*(
00629 diag_[2][i + 1]*source[i-1] +
00630 diag_[1][i]*source[i] +
00631 diag_[0][i - 1]*source[i + 1]);
00632 }
00633
00634
00635 dest[0] +=
00636 alpha*(
00637 diag_[1][0]*source[0] +
00638 diag_[0][1]*source[1] +
00639 diag_[0][n-1]*source[n - 1]);
00640
00641 dest[n - 1] +=
00642 alpha*(
00643 diag_[2][0]*source[0] +
00644 diag_[2][n-1]*source[n-2] +
00645 diag_[1][n-1]*source[n-1]);
00646 }
00647
00654 Scalar frobeniusNorm() const
00655 { return std::sqrt(frobeniusNormSquared()); }
00656
00662 Scalar frobeniusNormSquared() const
00663 {
00664 Scalar result = 0;
00665
00666 unsigned n = size();
00667 for (unsigned i = 0; i < n; ++ i)
00668 for (unsigned diagIdx = 0; diagIdx < 3; ++ diagIdx)
00669 result += diag_[diagIdx][i];
00670
00671 return result;
00672 }
00673
00679 Scalar infinityNorm() const
00680 {
00681 Scalar result=0;
00682
00683
00684 unsigned n = size();
00685 for (unsigned i = 1; i < n - 1; ++ i) {
00686 result = std::max(result,
00687 std::abs(diag_[0][i - 1]) +
00688 std::abs(diag_[1][i]) +
00689 std::abs(diag_[2][i + 1]));
00690 }
00691
00692
00693 result = std::max(result,
00694 std::abs(diag_[1][0]) +
00695 std::abs(diag_[2][1]) +
00696 std::abs(diag_[2][0]));
00697
00698
00699 result = std::max(result,
00700 std::abs(diag_[0][n-1]) +
00701 std::abs(diag_[0][n-2]) +
00702 std::abs(diag_[1][n-2]));
00703
00704 return result;
00705 }
00706
00713 template <class XVector, class BVector>
00714 void solve(XVector& x, const BVector& b) const
00715 {
00716 if (size() > 2 && std::abs(diag_[2][0]) < 1e-30)
00717 solveWithUpperRight_(x, b);
00718 else
00719 solveWithoutUpperRight_(x, b);
00720 }
00721
00725 void print(std::ostream& os = std::cout) const
00726 {
00727 size_t n = size();
00728
00729
00730 os << at(0, 0) << "\t"
00731 << at(0, 1) << "\t";
00732
00733 if (n > 3)
00734 os << "\t";
00735 if (n > 2)
00736 os << at(0, n-1);
00737 os << "\n";
00738
00739
00740 for (unsigned rowIdx = 1; rowIdx < n-1; ++rowIdx) {
00741 if (rowIdx > 1)
00742 os << "\t";
00743 if (rowIdx == n - 2)
00744 os << "\t";
00745
00746 os << at(rowIdx, rowIdx - 1) << "\t"
00747 << at(rowIdx, rowIdx) << "\t"
00748 << at(rowIdx, rowIdx + 1) << "\n";
00749 }
00750
00751
00752 if (n > 2)
00753 os << at(n-1, 0) << "\t";
00754 if (n > 3)
00755 os << "\t";
00756 if (n > 4)
00757 os << "\t";
00758 os << at(n-1, n-2) << "\t"
00759 << at(n-1, n-1) << "\n";
00760 }
00761
00762 private:
00763 template <class XVector, class BVector>
00764 void solveWithUpperRight_(XVector& x, const BVector& b) const
00765 {
00766 size_t n = size();
00767
00768 std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n);
00769 std::vector<Scalar> bStar(n);
00770 std::copy(b.begin(), b.end(), bStar.begin());
00771
00772 lastColumn[0] = upperDiag[0];
00773
00774
00775 for (size_t i = 1; i < n; ++i) {
00776 Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
00777
00778 lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
00779 mainDiag[i] -= alpha * upperDiag[i];
00780
00781 bStar[i] -= alpha * bStar[i - 1];
00782 };
00783
00784
00785 if (lowerDiag[n - 1] != 0.0 && size() > 2) {
00786 Scalar lastRow = lowerDiag[n - 1];
00787 for (size_t i = 0; i < n - 1; ++i) {
00788 Scalar alpha = lastRow/mainDiag[i];
00789 lastRow = - alpha*upperDiag[i + 1];
00790 bStar[n - 1] -= alpha * bStar[i];
00791 }
00792
00793 mainDiag[n-1] += lastRow;
00794 }
00795
00796
00797 x[n - 1] = bStar[n - 1]/mainDiag[n-1];
00798 for (int i = static_cast<int>(n) - 2; i >= 0; --i) {
00799 unsigned iu = static_cast<unsigned>(i);
00800 x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1] - x[n-1]*lastColumn[iu])/mainDiag[iu];
00801 }
00802 }
00803
00804 template <class XVector, class BVector>
00805 void solveWithoutUpperRight_(XVector& x, const BVector& b) const
00806 {
00807 size_t n = size();
00808
00809 std::vector<Scalar> lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]);
00810 std::vector<Scalar> bStar(n);
00811 std::copy(b.begin(), b.end(), bStar.begin());
00812
00813
00814 for (size_t i = 1; i < n; ++i) {
00815 Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1];
00816
00817 lowerDiag[i - 1] -= alpha * mainDiag[i - 1];
00818 mainDiag[i] -= alpha * upperDiag[i];
00819
00820 bStar[i] -= alpha * bStar[i - 1];
00821 };
00822
00823
00824 if (lowerDiag[n - 1] != 0.0 && size() > 2) {
00825 Scalar lastRow = lowerDiag[n - 1];
00826 for (size_t i = 0; i < n - 1; ++i) {
00827 Scalar alpha = lastRow/mainDiag[i];
00828 lastRow = - alpha*upperDiag[i + 1];
00829 bStar[n - 1] -= alpha * bStar[i];
00830 }
00831
00832 mainDiag[n-1] += lastRow;
00833 }
00834
00835
00836 x[n - 1] = bStar[n - 1]/mainDiag[n-1];
00837 for (int i = static_cast<int>(n) - 2; i >= 0; --i) {
00838 unsigned iu = static_cast<unsigned>(i);
00839 x[iu] = (bStar[iu] - x[iu + 1]*upperDiag[iu+1])/mainDiag[iu];
00840 }
00841 }
00842
00843 mutable std::vector<Scalar> diag_[3];
00844 };
00845
00846 }
00847
00848 template <class Scalar>
00849 std::ostream& operator<<(std::ostream& os, const Opm::TridiagonalMatrix<Scalar>& mat)
00850 {
00851 mat.print(os);
00852 return os;
00853 }
00854
00855 #endif