Merge pull request #56 from gaganjyot/EigenLinearSolver

Gauss Jordan to Eigen colPivHouseHolderQr for Linear System Solver
LibreCAD · May 23, 2016 · a7ece7a · a7ece7a
2 parents 90ffbdf + 980db53
commit a7ece7a
Show file tree

Hide file tree

Showing 3 changed files with 37 additions and 10 deletions.
diff --git a/lckernel/cad/math/intersectionhandler.cpp b/lckernel/cad/math/intersectionhandler.cpp
@@ -7,16 +7,26 @@ std::vector<geo::Coordinate> Intersection::LineLine(const Equation& l1,
     std::vector<lc::geo::Coordinate> ret;
     const auto &m1 = l1.Matrix();
     const auto &m2 = l2.Matrix();
-    std::vector<std::vector<double>> ce = {
-        // D, E, F
-        {m1(2,0) + m1(0,2), m1(2,1) + m1(1,2), -m1(2,2)},
-        {m2(2,0) + m2(0,2), m2(2,1) + m2(1,2), -m2(2,2)}
-    };
-
-    std::vector<double> sn(2, 0.);
-    if (Math::linearSolver(ce, sn)) {
-        ret.emplace_back(sn[0], sn[1]);
-    }
+//    std::vector<std::vector<double>> ce = {
+//        // D, E, F
+//        {m1(2,0) + m1(0,2), m1(2,1) + m1(1,2), -m1(2,2)},
+//        {m2(2,0) + m2(0,2), m2(2,1) + m2(1,2), -m2(2,2)}
+//    };
+
+//    std::vector<double> sn(2, 0.);
+//    if (Math::linearSolver(ce, sn)) {
+//        ret.emplace_back(sn[0], sn[1]);
+//    }
+//    return ret;
+    Eigen::Matrix2d M;
+    Eigen::Vector2d V;
+
+    M << m1(2,0) + m1(0,2), m1(2,1) + m1(1,2),
+         m2(2,0) + m2(0,2), m2(2,1) + m2(1,2);
+    V << -m1(2,2), -m2(2,2);
+
+    Eigen::Vector2d R = M.colPivHouseholderQr().solve(V);
+    ret.emplace_back(R[0], R[1]);
     return ret;
 }
 

diff --git a/lckernel/cad/math/intersectionhandler.h b/lckernel/cad/math/intersectionhandler.h
@@ -2,6 +2,7 @@
 
 #include "cad/geometry/geocoordinate.h"
 #include "cad/math/equation.h"
+#include "eigen3/Eigen/Dense"
 
 namespace lc {
     namespace maths {

diff --git a/unittest/testmatrices.h b/unittest/testmatrices.h
@@ -2,6 +2,7 @@
 #include "cad/math/quadratic_math.h"
 #include "cad/math/intersectionhandler.h"
 #include "cad/geometry/geovector.h"
+#include "cad/math/lcmath.h"
 #include <iostream>
 
 using namespace std;
@@ -116,3 +117,18 @@ TEST(QM, LineQuad) {
         EXPECT_DOUBLE_EQ(x[i].y(), y[i].y()) << "Y differs at index " << i;
     }
 }
+
+TEST(Maths, LinearSolver) {
+    std::vector<std::vector<double>> a {{1,-1, -1}, {3,1, 9}};
+    std::vector<double> res1;
+    Eigen::Matrix2d M;
+    Eigen::Vector2d V;
+
+    M << 1, -1, 3, 1;
+    V << -1, 9;
+    lc::Math::linearSolver(a, res1);
+    Eigen::Vector2d res2 = M.colPivHouseholderQr().solve(V);
+
+    EXPECT_DOUBLE_EQ(res1[0], res2[0]) << "X differs";
+    EXPECT_DOUBLE_EQ(res1[1], res2[1]) << "Y differs";
+}