feat: add JacobiMethod impl

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++) {

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);
 

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<JacobiIterationPivot> 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};
+    }
+}

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<JacobiMethodOutput>;
 
+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

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) {

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);
 

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>{ //
-            Vector(std::vector<lld>{1, -5, 3}), //
-            Vector(std::vector<lld>{0, -3, -1}), //
-            Vector(std::vector<lld>{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<std::vector<lld>>{ //
-            {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<lld>{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<std::vector<lld>>{
-//            {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<std::vector<lld>>{
-//            {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<lld>{0.9, 1.0, 1.1};
-    for (auto x: x_list) {
-        std::cout << "x = " << x << std::endl;
-        auto m = Matrix(std::vector<std::vector<lld>>{{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;
 }