add 添加矩阵LU分解

simplify-life · Sep 18, 2023 · 2b610f1 · 2b610f1
1 parent 23fa59a
commit 2b610f1
Show file tree

Hide file tree

Showing 4 changed files with 70 additions and 2 deletions.
diff --git a/include/math/matrix.h b/include/math/matrix.h
@@ -87,7 +87,7 @@ namespace xlib {
             // 计算矩阵的迹
             double trace() const;
 
-
+            bool operator==(const Matrix& mat) const;
     };
 
     std::vector<double> normalGaussianElimination(std::vector<std::vector<double>>& A, std::vector<double>& b);
@@ -99,6 +99,9 @@ namespace xlib {
     std::ostream& operator<<(std::ostream& os, const Matrix& matrix);
 
     std::string toString(const std::vector<std::vector<double>>& matrix);
+
+    //矩阵的LU 分解
+    void luDecomposition(const std::vector<std::vector<double>>& A, std::vector<std::vector<double>>& L, std::vector<std::vector<double>>& U);
 }
 
 #endif /* matrix_h */
diff --git a/main.cpp b/main.cpp
@@ -487,6 +487,23 @@ void testMatrix()
         std::cout << eigenvalue << " ";
     }
     std::cout << std::endl;
+
+    std::vector<std::vector<double>> LU_A = {{2, -1, 3},
+                                          {4, 2, -2},
+                                          {-2, 3, 1}};
+    std::vector<std::vector<double>> LU_L, LU_U;
+
+    luDecomposition(LU_A, LU_L, LU_U);
+
+    Matrix mL = LU_L;
+    Matrix mU = LU_U;
+
+    Matrix mA = LU_A;
+    LOG_I("LU_A:\n%s",mA.toString().c_str());
+    LOG_I("L:\n%s,\nU:\n%s",mL.toString().c_str(),mU.toString().c_str());
+    auto LU_A_NEW = mL*mU;
+    LOG_I("LU_A_NEW:\n%s",LU_A_NEW.toString().c_str());
+    LOG_I("LU_A_NEW %s LU_A",mA==LU_A_NEW?"equals":"not equals");
 }
 
 void testMath()

diff --git a/src/base/XUtf8.cpp b/src/base/XUtf8.cpp
@@ -91,7 +91,7 @@ string XUtf8::utf8ToUnicode(const string& src)
 
 string singleUnicodeToUtf8(const std::string& src){
     std::string str = "0x";
-    str = str.append(src);
+    str.append(src);
     char *tmp;
 //    int place = 0;
     auto utf8code = std::strtol(str.c_str(), &tmp , 16);

diff --git a/src/math/matrix.cpp b/src/math/matrix.cpp
@@ -695,6 +695,20 @@ namespace xlib {
         return result;
     }
 
+    bool Matrix::operator==(const Matrix& mat) const{
+        if (m != mat.m || n != mat.n) {
+            return false;
+        }
+        for (int i = 0; i < m; i++) {
+            for (int j = 0; j < n; j++) {
+                if (a[i][j] != mat.a[i][j]) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
     double Matrix::trace() const {
         double result = 0;
         for (int i = 0; i < n; i++) {
@@ -714,4 +728,38 @@ namespace xlib {
         return os;
     }
 
+    // LU分解函数
+    void luDecomposition(const std::vector<std::vector<double>>& A, std::vector<std::vector<double>>& L, std::vector<std::vector<double>>& U) {
+        int n = A.size();
+        L.resize(n, std::vector<double>(n, 0.0));
+        U.resize(n, std::vector<double>(n, 0.0));
+        for (int i = 0; i < n; i++) {
+            U[0][i] = A[0][i];
+            L[i][0] = A[i][0] / U[0][0];
+        }
+       //计算U的第r行，L的第r列元素
+        for (int r = 1; r < n; r++) {
+            for (int i = r; i < n; i++) {
+                double sLU = 0.0;
+                for (int k = 0; k < r; k++) {
+                    sLU += L[r][k] * U[k][i];
+                }
+                U[r][i] = A[r][i] - sLU;
+                if(i==r){
+                    L[r][r]=1;
+                }
+                else if(r==n){
+                    L[n][n]=1;  
+                } 
+                else{
+                    double sLU = 0.0;
+                    for (int k = 0; k < r; k++) {
+                        sLU += L[i][k] * U[k][r];
+                    }
+                    L[i][r] = (A[i][r] - sLU) / U[r][r];
+                } 
+            }
+        }
+    }
+
 }