solvespace
diff --git a/‎src/system.cpp‎
Lines changed: 51 additions & 5 deletions b/‎src/system.cpp‎
Lines changed: 51 additions & 5 deletions
diff --git a/‎test/debugtool.cpp‎
Lines changed: 72 additions & 1 deletion b/‎test/debugtool.cpp‎
Lines changed: 72 additions & 1 deletion
@@ -291,13 +291,59 @@ bool System::SolveLeastSquares() {
         }
     }
 
-    SparseMatrix<double> AAt = mat.A.num * mat.A.num.transpose();
-    AAt.makeCompressed();
-    VectorXd z(mat.n);
+    if(mat.m == 0 || mat.n == 0) {
+        mat.X = VectorXd::Zero(mat.n);
+        return true;
+    }
+
+    if(mat.m <= mat.n) {
+        // We want the minimum norm solution of the underdetermined system
+        // A*x = B, which is x = A'*(A*A')^-1*B. Forming the normal equations
+        // A*A' explicitly squares the condition number of A; since our
+        // Jacobian mixes dimensionless quantities with lengths, areas and
+        // volumes, its condition number already grows with the size of the
+        // sketch, and squaring it is what makes large models fail. Eigen's
+        // rank-revealing QR then drops pivots below a threshold proportional
+        // to the largest column norm, and its rank-truncated solve silently
+        // zeroes the corresponding component of the Newton step, so the
+        // residual of one equation can never be driven to zero and we report
+        // unsolvable constraints for a perfectly solvable sketch.
+        //
+        // Factor A' instead, as A'*P = Q*R. Then A*A' = P*R'*R*P', so
+        //     z = (A*A')^-1*B = P*R^-1*R'^-1*P'*B,
+        // which is two triangular solves, and the rank decision is taken on
+        // R (whose pivots scale like A, not like A*A').
+        SparseMatrix<double> At = mat.A.num.transpose();
+        At.makeCompressed();
+
+        SparseQR<SparseMatrix<double>, COLAMDOrdering<int>> solver;
+        solver.compute(At);
+        if(solver.info() != Success) return false;
+
+        const int rank = (int)solver.rank();
+        VectorXd rhs = solver.colsPermutation().transpose() * mat.B.num;
+        VectorXd v   = VectorXd::Zero(mat.m);
+        if(rank > 0) {
+            SparseMatrix<double> R = solver.matrixR().topLeftCorner(rank, rank);
+            VectorXd w = R.triangularView<Upper>().transpose()
+                             .solve(rhs.topRows(rank));
+            v.topRows(rank) = R.triangularView<Upper>().solve(w);
+        }
+        VectorXd z = solver.colsPermutation() * v;
 
-    if(!SolveLinearSystem(AAt, mat.B.num, &z)) return false;
+        mat.X = mat.A.num.transpose() * z;
+    } else {
+        // More equations than unknowns; A' is wide, which Eigen's sparse QR
+        // does not handle, so fall back to the normal equations. A system
+        // like that is redundant anyway, and gets reported as such.
+        SparseMatrix<double> AAt = mat.A.num * mat.A.num.transpose();
+        AAt.makeCompressed();
+        VectorXd z(mat.m);
+
+        if(!SolveLinearSystem(AAt, mat.B.num, &z)) return false;
 
-    mat.X = mat.A.num.transpose() * z;
+        mat.X = mat.A.num.transpose() * z;
+    }
 
     for(int c = 0; c < mat.n; c++) {
         mat.X[c] *= scale[c];
 
@@ -4,15 +4,82 @@
 // Copyright 2017 whitequark
 //-----------------------------------------------------------------------------
 
+#include "solvespace.h"
 #include "expr.h"
 #include "platform/platform.h"
 
 using namespace SolveSpace;
 
+static const char *SolveResultName(SolveResult how) {
+    switch(how) {
+        case SolveResult::OKAY:                     return "OKAY";
+        case SolveResult::DIDNT_CONVERGE:           return "DIDNT_CONVERGE";
+        case SolveResult::REDUNDANT_OKAY:           return "REDUNDANT_OKAY";
+        case SolveResult::REDUNDANT_DIDNT_CONVERGE: return "REDUNDANT_DIDNT_CONVERGE";
+        case SolveResult::TOO_MANY_UNKNOWNS:        return "TOO_MANY_UNKNOWNS";
+    }
+    return "?";
+}
+
+// Load a file, regenerate it, and report how every group solved. Useful to
+// reproduce solver failures without a GUI.
+static int CmdSolve(const std::string &filename) {
+    SS.Init();
+    SS.showToolbar = false;
+    SS.checkClosedContour = false;
+
+    if(!SS.LoadFromFile(Platform::Path::From(filename))) {
+        fprintf(stderr, "cannot load: %s\n", filename.c_str());
+        return 1;
+    }
+    SS.AfterNewFile();
+
+    int failed = 0;
+    for(Group &g : SK.group) {
+        bool ok = g.IsSolvedOkay();
+        if(!ok) failed++;
+        fprintf(stderr, "group %08x %-24s %-24s dof=%d%s\n", g.h.v,
+                g.DescriptionString().c_str(), SolveResultName(g.solved.how),
+                g.solved.dof, ok ? "" : "  FAILED");
+        for(int i = 0; i < g.solved.remove.n; i++) {
+            Constraint *c = SK.constraint.FindByIdNoOops(g.solved.remove[i]);
+            if(!c) continue;
+            fprintf(stderr, "    bad constraint %08x %s\n", c->h.v,
+                    c->DescriptionString().c_str());
+        }
+    }
+    fprintf(stderr, "%s: %d group(s) failed to solve\n", filename.c_str(), failed);
+    return failed == 0 ? 0 : 1;
+}
+
+// Load a file, regenerate it, and write it back out. This is what happens when
+// a part that other files link to is opened, edited and saved again.
+static int CmdResave(const std::string &filename) {
+    SS.Init();
+    SS.showToolbar = false;
+    SS.checkClosedContour = false;
+
+    Platform::Path path = Platform::Path::From(filename);
+    if(!SS.LoadFromFile(path)) {
+        fprintf(stderr, "cannot load: %s\n", filename.c_str());
+        return 1;
+    }
+    SS.AfterNewFile();
+    if(!SS.SaveToFile(path)) {
+        fprintf(stderr, "cannot save: %s\n", filename.c_str());
+        return 1;
+    }
+    return 0;
+}
+
 int main(int argc, char **argv) {
     std::vector<std::string> args = Platform::InitCli(argc, argv);
 
-    if(args.size() == 3 && args[1] == "expr") {
+    if(args.size() == 3 && args[1] == "solve") {
+        return CmdSolve(args[2]);
+    } else if(args.size() == 3 && args[1] == "resave") {
+        return CmdResave(args[2]);
+    } else if(args.size() == 3 && args[1] == "expr") {
         std::string expr = args[2], err;
         Expr *e = Expr::Parse(expr.c_str(), &err);
         if(e == NULL) {
@@ -28,6 +95,10 @@ int main(int argc, char **argv) {
 Commands:
     expr [expr]
         Evaluate an expression.
+    solve [file.slvs]
+        Load a file and report how each group solved.
+    resave [file.slvs]
+        Load a file, regenerate it, and save it back.
 )");
     }