diff --git a/CustomMath_lib/InfinityNorm/InfinityNorm.cpp b/CustomMath_lib/InfinityNorm/InfinityNorm.cpp index dd5c52d..819ac6f 100644 --- a/CustomMath_lib/InfinityNorm/InfinityNorm.cpp +++ b/CustomMath_lib/InfinityNorm/InfinityNorm.cpp @@ -18,6 +18,18 @@ lld MatrixNorm_Infinity(const Matrix &A) { return max; } +lld MatrixNorm_E(const Matrix &A) { + lld sum = 0; + for (int i = 0; i < A.rows; i++) { + for (int j = 0; j < A.cols; j++) { + if (i != j) { + sum += std::abs(A.matrix[i][j]); + } + } + } + return sum; +} + lld VectorNorm_Infinity(const Vector &x) { lld max = 0; for (int i = 0; i < x.size; i++) { diff --git a/CustomMath_lib/InfinityNorm/InfinityNorm.h b/CustomMath_lib/InfinityNorm/InfinityNorm.h index fae4e2b..7ecb93f 100644 --- a/CustomMath_lib/InfinityNorm/InfinityNorm.h +++ b/CustomMath_lib/InfinityNorm/InfinityNorm.h @@ -10,6 +10,9 @@ /// Calculate infinity norm of a matrix lld MatrixNorm_Infinity(const Matrix &A); +/// Calculate E(A) which is the non-diagonal norm of A +lld MatrixNorm_E(const Matrix &A); + /// Calculate infinity norm of a vector lld VectorNorm_Infinity(const Vector &x); diff --git a/CustomMath_lib/JacobiMethod/JacobiMethod.cpp b/CustomMath_lib/JacobiMethod/JacobiMethod.cpp index 24c1c08..fe748b8 100644 --- a/CustomMath_lib/JacobiMethod/JacobiMethod.cpp +++ b/CustomMath_lib/JacobiMethod/JacobiMethod.cpp @@ -3,3 +3,151 @@ // #include "JacobiMethod.h" + +Matrix JacobiRotationMatrix(ull n, ull p, ull q, lld c, lld s) { + auto P = Matrix::identity(n); + P.matrix[p][p] = c; + P.matrix[q][q] = c; + P.matrix[p][q] = -s; + P.matrix[q][p] = s; + + return P; +} + +typedef struct { + ull p; + ull q; +} JacobiIterationPivot; + +JacobiIterationPivot JacobiIterationPivotFinder(const Matrix &A) { + auto n = A.rows; + ull p = 0, q = 1; // here i < j is always assumed + lld pivot = std::abs(A.matrix[0][1]); + for (ull i = 0; i < n; i++) { + for (ull j = i + 1; j < n; j++) { + if (std::abs(A.matrix[i][j]) > std::abs(pivot)) { + pivot = A.matrix[i][j]; + p = i; + q = j; + } + } + } + return {p, q}; +} + +std::optional JacobiIterationPivotFinder(const Matrix &A, lld threshold) { + auto n = A.rows; + for (ull i = 0; i < n; i++) { + for (ull j = i + 1; j < n; j++) { + if (std::abs(A.matrix[i][j]) > threshold) { + return {JacobiIterationPivot{i, j}}; + } + } + } + return std::nullopt; +} + +void JacobiIteration(ull n, Matrix &A, Matrix &P, JacobiIterationPivot pivot) { + auto p = pivot.p; + auto q = pivot.q; + + auto a_pp = A.matrix[p][p]; + auto a_qq = A.matrix[q][q]; + auto a_pq = A.matrix[p][q]; + + auto rho = (a_pp - a_qq) / (2 * a_pq); + + auto t = SIGN(rho) / (std::abs(rho) + std::sqrt(1 + rho * rho)); // tan(theta) + auto c = 1 / std::sqrt(1 + t * t); // cos(theta) + auto s = c * t; // sin(theta) + + + auto J = JacobiRotationMatrix(n, p, q, c, s); + P = P * J; + + // manually calculate A = P^T A P to speed up + A.matrix[p][p] = c * c * a_pp + s * s * a_qq + 2 * c * s * a_pq; // beta_pp + A.matrix[q][q] = s * s * a_pp + c * c * a_qq - 2 * c * s * a_pq; // beta_qq + A.matrix[p][q] = 0; // beta_pq + A.matrix[q][p] = 0; // beta_qp + + + // update beta_ip and beta_iq + for (ull i = 0; i < n; i++) { + if (i != p && i != q) { + auto a_ip = A.matrix[i][p]; + auto a_iq = A.matrix[i][q]; + A.matrix[i][p] = c * a_ip + s * a_iq; + A.matrix[p][i] = A.matrix[i][p]; + A.matrix[i][q] = c * a_iq - s * a_ip; + A.matrix[q][i] = A.matrix[i][q]; + } + } +} + +JIterationMethodOutput JacobiMethod(const Matrix &matrix, JacobiMethodType type) { + auto n = matrix.rows; + auto A = matrix; + auto P = Matrix::identity(n); + + if (type == JACOBI_METHOD_TYPE_CLASSIC) { + auto iteration_count = 0; + ITERATION_METHOD_TIMING_START + + while (true) { + auto pivot = JacobiIterationPivotFinder(A); + JacobiIteration(n, A, P, pivot); + iteration_count++; + + if (MatrixNorm_E(A) < 1e-10) { + + ITERATION_METHOD_TIMING_END + + return {A, P, iteration_count, ITERATION_METHOD_RETURN_DURATION}; + } + } + } else if (type == JACOBI_METHOD_TYPE_LOOP) { + auto iteration_count = 0; + ITERATION_METHOD_TIMING_START + + while (true) { + for (ull i = 0; i < n - 1; i++) { + for (ull j = i + 1; j < n; j++) { + auto pivot = JacobiIterationPivot{i, j}; + JacobiIteration(n, A, P, pivot); + iteration_count++; + + if (MatrixNorm_E(A) < 1e-10) { + + ITERATION_METHOD_TIMING_END + + return {A, P, iteration_count, ITERATION_METHOD_RETURN_DURATION}; + } + } + } + } + } else if (type == JACOBI_METHOD_TYPE_THRESHOLD) { + auto delta = MatrixNorm_E(A); + auto iteration_count = 0; + + ITERATION_METHOD_TIMING_START + + while (delta > 1e-10) { + while (true) { + auto pivot = JacobiIterationPivotFinder(A, delta); + if (pivot.has_value()) { + iteration_count++; + JacobiIteration(n, A, P, pivot.value()); + } else { + break; + } + } + + delta /= 10; + } + + ITERATION_METHOD_TIMING_END + + return {A, P, iteration_count, ITERATION_METHOD_RETURN_DURATION}; + } +} \ No newline at end of file diff --git a/CustomMath_lib/JacobiMethod/JacobiMethod.h b/CustomMath_lib/JacobiMethod/JacobiMethod.h index 6dee044..177331c 100644 --- a/CustomMath_lib/JacobiMethod/JacobiMethod.h +++ b/CustomMath_lib/JacobiMethod/JacobiMethod.h @@ -8,11 +8,18 @@ #include "CustomMath_lib.h" typedef struct { - Matrix H; + Matrix A; Matrix P; } JacobiMethodOutput; +using JIterationMethodOutput = IterationMethodOutput; +enum JacobiMethodType { + JACOBI_METHOD_TYPE_CLASSIC, JACOBI_METHOD_TYPE_LOOP, JACOBI_METHOD_TYPE_THRESHOLD +}; +/// P^T A P = H, returns H, P +/// A must be symmetric +JIterationMethodOutput JacobiMethod(const Matrix &matrix, JacobiMethodType type = JACOBI_METHOD_TYPE_CLASSIC); #endif //NUMERICAL_ALGEBRA_JACOBIMETHOD_H diff --git a/CustomMath_lib/Matrix/MatrixProperties.cpp b/CustomMath_lib/Matrix/MatrixProperties.cpp index 11bc2da..bc9007e 100644 --- a/CustomMath_lib/Matrix/MatrixProperties.cpp +++ b/CustomMath_lib/Matrix/MatrixProperties.cpp @@ -5,14 +5,13 @@ #include "Matrix.h" void Matrix::print() { + std::cout << std::fixed << std::setprecision(6); + std::cout << "["; for (ull i = 0; i < this->rows; i++) { std::cout << "["; for (ull j = 0; j < this->cols; j++) { - std::cout << this->matrix[i][j]; - if (j != this->cols - 1) { - std::cout << ", "; - } + std::cout << "\t" << this->matrix[i][j] << "\t"; } std::cout << "]"; if (i != this->rows - 1) { diff --git a/CustomMath_lib/base.h b/CustomMath_lib/base.h index 66d06fe..dc89f63 100644 --- a/CustomMath_lib/base.h +++ b/CustomMath_lib/base.h @@ -15,12 +15,14 @@ #include "string" #include "random" #include "chrono" +#include "optional" #define ull unsigned long long #define lld long double #define MAX(a, b) ((a) > (b) ? (a) : (b)) #define MIN(a, b) ((a) < (b) ? (a) : (b)) +#define SIGN(a) ((a) > 0 ? 1 : ((a) < 0 ? -1 : 0)) bool cmp(lld a, lld b); diff --git a/main.cpp b/main.cpp index 4b542ca..98cc8f3 100644 --- a/main.cpp +++ b/main.cpp @@ -1,105 +1,18 @@ #include "CustomMath_lib.h" -void par_1_each(Vector &a) { - auto max_root = MaxRootForPolynomial(a); - - std::cout << "Polynomial is "; - a.print(); - std::cout << "Max root for polynomial is " << std::setprecision(10) << max_root.result << std::endl; - std::cout << "Iteration times is " << max_root.iteration_count << std::endl; - std::cout << "Time cost is " << max_root.time_cost.count() << " microseconds" << std::endl; -} - - -void par_1() { - std::cout << "------ Q 6.1 ------" << std::endl; - - auto a_list = std::vector{ // - Vector(std::vector{1, -5, 3}), // - Vector(std::vector{0, -3, -1}), // - Vector(std::vector{101, 208.01, 10891.01, 9802.08, 79108.9, -99902, 790, -1000}) // - }; - - for (auto &a: a_list) { - par_1_each(a); - } -} - -void hessenberg_test() { - auto m = Matrix(std::vector>{ // - {1.0, 2.9, 3.8, 4.7}, // - {5.6, 6.5, 7.4, 8.3}, // - {9.2, 10.1, 11.0, 12.9}, // - {13.8, 14.7, 15.6, 16.5} // - }); - m.print(); - Matrix u; - HessenbergMethod_Inplace(m, u); - m.print(); - u.print(); - auto test = u.transpose() * m * u; - test.print(); -} - -void par_2_2() { - std::cout << std::endl << "------ Q 6.2(2) ------" << std::endl; - -// auto coefficients = Vector(std::vector{1, -5, 3}); - auto coefficients = Vector(41); - coefficients.array[41 - 3 + 1] = 1; - coefficients.array[41 - 1] = 1; - - auto r = AllRootsForPolynomial(coefficients); - std::cout << "Roots for polynomial are:" << std::endl; - r.result.print(); - std::cout << "Iteration times is " << r.iteration_count << std::endl; - std::cout << "Time cost is " << r.time_cost.count() << " microseconds" << std::endl; -} - -void par_2_3() { - std::cout << std::endl << "------ Q 6.2(3) ------" << std::endl; +int main() { + for (auto type: {JACOBI_METHOD_TYPE_CLASSIC, JACOBI_METHOD_TYPE_LOOP, JACOBI_METHOD_TYPE_THRESHOLD}) { + auto m = Matrix("[5 1 -2; 1 2 0; -2 0 -10]"); + auto r = JacobiMethod(m, type); -// auto m = Matrix(std::vector>{ -// {1.0, 2.9, 3.8, 4.7}, -// {5.6, 6.5, 7.4, 8.3}, -// {9.2, 10.1, 11.0, 12.9}, -// {13.8, 14.7, 15.6, 16.5} -// }); -// auto m = Matrix(std::vector>{ -// {2,3,4,5,6}, -// {4,4,5,6,7}, -// {0,3,6,7,8}, -// {0,0,2,8,9}, -// {0,0,0,1,10} -// }); - auto x_list = std::vector{0.9, 1.0, 1.1}; - for (auto x: x_list) { - std::cout << "x = " << x << std::endl; - auto m = Matrix(std::vector>{{9.1, 3.0, 2.6, 4.0}, - {4.2, 5.3, 4.7, 1.6}, - {3.2, 1.7, 9.4, x}, - {6.1, 4.9, 3.5, 6.2}}); + r.result.A.print(); + r.result.P.print(); - auto h = QRMethod(m); - auto k = AllEigenValues(h.result); - std::cout << "Eigen values are:" << std::endl; - print_llc(k); - std::cout << "Iteration times is " << h.iteration_count << std::endl; - std::cout << "Time cost is " << h.time_cost.count() << " microseconds" << std::endl; + (r.result.P.transpose() * r.result.A * r.result.P).print(); - std::cout << std::endl; + std::cout << r.iteration_count << std::endl; + std::cout << r.time_cost.count() << std::endl; } -} - -void par_2() { - par_2_2(); - par_2_3(); -} - -int main() { -// hessenberg_test(); - par_1(); - par_2(); return 0; } \ No newline at end of file