123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173 |
- #include <Eigen/Eigen>
- // The banded system class is used for solving
- // banded linear system Ax=b efficiently.
- // A is an N*N band matrix with lower band width lowerBw
- // and upper band width upperBw.
- // Banded LU factorization has O(N) time complexity.
- // 带状系统类用于解决
- // 有效的带状线性系统 Ax=b。
- // A 是一个 N*N 带矩阵,具有较低带宽 lowerBw
- // 和上带宽 upperBw。
- // 带状 LU 分解的时间复杂度为 O(N)。
- #ifndef BANDEDSYSTEM_H
- #define BANDEDSYSTEM_H
- class BandedSystem // 这里没太看懂
- {
- public:
- // The size of A, as well as the lower/upper
- // banded width p/q are needed
- inline void create(const int &n, const int &p, const int &q)
- {
- // In case of re-creating before destroying
- destroy();
- N = n;
- lowerBw = p;
- upperBw = q;
- int actualSize = N * (lowerBw + upperBw + 1);
- ptrData = new double[actualSize];
- std::fill_n(ptrData, actualSize, 0.0);
- return;
- }
- inline void destroy()
- {
- if (ptrData != nullptr)
- {
- delete[] ptrData;
- ptrData = nullptr;
- }
- return;
- }
- private:
- int N;
- int lowerBw;
- int upperBw;
- // Compulsory nullptr initialization here
- double *ptrData = nullptr;
- public:
- // Reset the matrix to zero
- inline void reset(void) // 这个reset有问题,只有主对角线和设置的upbound lowbound对应的位置是0,其他都是不确定的
- {
- std::fill_n(ptrData, N * (lowerBw + upperBw + 1), 0.0);
- return;
- }
- // The band matrix is stored as suggested in "Matrix Computation"
- inline const double &operator()(const int &i, const int &j) const
- {
- return ptrData[(i - j + upperBw) * N + j];
- }
- inline double &operator()(const int &i, const int &j)
- {
- return ptrData[(i - j + upperBw) * N + j];
- }
- // This function conducts banded LU factorization in place
- // Note that NO PIVOT is applied on the matrix "A" for efficiency!!!
- inline void factorizeLU()
- {
- int iM, jM;
- double cVl;
- for (int k = 0; k <= N - 2; ++k)
- {
- iM = std::min(k + lowerBw, N - 1);
- cVl = operator()(k, k);
- for (int i = k + 1; i <= iM; ++i)
- {
- if (operator()(i, k) != 0.0)
- {
- operator()(i, k) /= cVl;
- }
- }
- jM = std::min(k + upperBw, N - 1);
- for (int j = k + 1; j <= jM; ++j)
- {
- cVl = operator()(k, j);
- if (cVl != 0.0)
- {
- for (int i = k + 1; i <= iM; ++i)
- {
- if (operator()(i, k) != 0.0)
- {
- operator()(i, j) -= operator()(i, k) * cVl;
- }
- }
- }
- }
- }
- return;
- }
- // This function solves Ax=b, then stores x in b
- // The input b is required to be N*m, i.e.,
- // m vectors to be solved.
- template <typename EIGENMAT>
- inline void solve(EIGENMAT &b) const
- {
- int iM;
- for (int j = 0; j <= N - 1; ++j)
- {
- iM = std::min(j + lowerBw, N - 1);
- for (int i = j + 1; i <= iM; ++i)
- {
- if (operator()(i, j) != 0.0)
- {
- b.row(i) -= operator()(i, j) * b.row(j);
- }
- }
- }
- for (int j = N - 1; j >= 0; --j)
- {
- b.row(j) /= operator()(j, j);
- iM = std::max(0, j - upperBw);
- for (int i = iM; i <= j - 1; ++i)
- {
- if (operator()(i, j) != 0.0)
- {
- b.row(i) -= operator()(i, j) * b.row(j);
- }
- }
- }
- return;
- }
- // This function solves ATx=b, then stores x in b
- // The input b is required to be N*m, i.e.,
- // m vectors to be solved.
- template <typename EIGENMAT>
- inline void solveAdj(EIGENMAT &b) const
- {
- int iM;
- for (int j = 0; j <= N - 1; ++j)
- {
- b.row(j) /= operator()(j, j);
- iM = std::min(j + upperBw, N - 1);
- for (int i = j + 1; i <= iM; ++i)
- {
- if (operator()(j, i) != 0.0)
- {
- b.row(i) -= operator()(j, i) * b.row(j);
- }
- }
- }
- for (int j = N - 1; j >= 0; --j)
- {
- iM = std::max(0, j - lowerBw);
- for (int i = iM; i <= j - 1; ++i)
- {
- if (operator()(j, i) != 0.0)
- {
- b.row(i) -= operator()(j, i) * b.row(j);
- }
- }
- }
- }
- };
- #endif // BANDEDSYSTEM_H
|