From b056b772aa3c81bb2c98177cba02fd36d829a457 Mon Sep 17 00:00:00 2001 From: BoykoNeov Date: Mon, 27 Jul 2026 22:23:26 +0300 Subject: [PATCH 1/2] Solve the Newton step without forming the normal equations. SolveLeastSquares() computed the minimum norm Newton step as x = A'*(A*A')^-1*B, forming A*A' explicitly and factoring it with a rank-revealing sparse QR. That squares the condition number of the Jacobian. Our equations mix dimensionless quantities with lengths and with areas, so the spread of magnitudes in A already grows with the physical size of the sketch; squaring it pushes the smallest pivot of A*A' below the threshold Eigen uses to call a column linearly dependent, which is proportional to the largest column norm. Eigen's rank-truncated solve then silently zeroes that component of the step, so the residual of one equation can never be driven to zero, Newton's method stalls, and a perfectly solvable sketch is reported as having incompatible constraints. On the file from the bug report, two 4 m lines constrained perpendicular with a point on line, the transition is exact: with the second line pinned at 4 m and the first at 2880 mm, the fifth pivot of A*A' is 3.657e-7 against a threshold of 3.6557e-7 and the sketch solves; at 2881 mm the pivot falls just below the threshold, the rank drops from 5 to 4, and the solve fails, with the step norm collapsing to 1e-13 while one residual stays pinned at 0.019. Factor A' = Q*R directly instead. Then A*A' = P*R'*R*P', so the same z = (A*A')^-1*B comes from two triangular solves against R, and the rank decision is taken on pivots that scale like A rather than like A*A'. In that same sweep, the longest first line that solves goes from 2.88 m to 76.9 km. Rank determination for redundant constraints is untouched; it runs in TestRank(), against A itself. Over-constrained sketches are in fact reported better than before: duplicating a constraint in that file used to be reported as redundant below about 3 m but as unsolvable above it, and is now reported as redundant at every size, out to 100 m. The wide case, more equations than unknowns, keeps using the normal equations, since Eigen's sparse QR wants a matrix that is at least as tall as it is wide. Such a system is redundant anyway, and it gives the same result as before. A system with no equations or no unknowns now returns a zero step, rather than multiplying an uninitialized vector by a matrix of a mismatched size. Also teach the debug tool to load a file and report how each group solved, and to load a file and save it back, which is how a linked part gets re-solved without a GUI. Fixes the original report and mesr's assembly from #1354. Co-Authored-By: Claude Opus 5 (1M context) --- src/system.cpp | 56 +++++++++++++++++++++++++++++++---- test/debugtool.cpp | 73 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 123 insertions(+), 6 deletions(-) diff --git a/src/system.cpp b/src/system.cpp index 89d315938..ded9e9eab 100644 --- a/src/system.cpp +++ b/src/system.cpp @@ -291,13 +291,59 @@ bool System::SolveLeastSquares() { } } - SparseMatrix 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 At = mat.A.num.transpose(); + At.makeCompressed(); + + SparseQR, COLAMDOrdering> 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 R = solver.matrixR().topLeftCorner(rank, rank); + VectorXd w = R.triangularView().transpose() + .solve(rhs.topRows(rank)); + v.topRows(rank) = R.triangularView().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 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]; diff --git a/test/debugtool.cpp b/test/debugtool.cpp index aea6b8475..41d973a88 100644 --- a/test/debugtool.cpp +++ b/test/debugtool.cpp @@ -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 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. )"); } From b473a13fef3d4e9d4152ead951a81c2d76bed4f7 Mon Sep 17 00:00:00 2001 From: BoykoNeov Date: Mon, 27 Jul 2026 22:23:39 +0300 Subject: [PATCH 2/2] Add a regression test for a sketch with large dimensions. The fixture is the file from the bug report, unmodified: two 4 m lines constrained perpendicular, with the start of the second one on the first. It has to be the file as saved, not a canonical re-save, because the failure depends on the stored parameter values being the ones from before the perpendicular constraint was added; re-saving it stores the solution, and then there is no Newton step left to take. The test checks the geometry as well as the solve result, since the constraints could be satisfied by more than one configuration. Co-Authored-By: Claude Opus 5 (1M context) --- test/CMakeLists.txt | 1 + .../large_dimensions/perpendicular_4m.slvs | 381 ++++++++++++++++++ test/constraint/large_dimensions/test.cpp | 42 ++ 3 files changed, 424 insertions(+) create mode 100644 test/constraint/large_dimensions/perpendicular_4m.slvs create mode 100644 test/constraint/large_dimensions/test.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0dafd42db..51d1cbd80 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,6 +54,7 @@ set(testsuite_SOURCES constraint/equal_radius/test.cpp constraint/where_dragged/test.cpp constraint/comment/test.cpp + constraint/large_dimensions/test.cpp request/arc_of_circle/test.cpp request/circle/test.cpp request/cubic/test.cpp diff --git a/test/constraint/large_dimensions/perpendicular_4m.slvs b/test/constraint/large_dimensions/perpendicular_4m.slvs new file mode 100644 index 000000000..bde75e04f --- /dev/null +++ b/test/constraint/large_dimensions/perpendicular_4m.slvs @@ -0,0 +1,381 @@ +±²³SolveSpaceREVa + + +Group.h.v=00000001 +Group.type=5000 +Group.name=#references +Group.color=ff000000 +Group.skipFirst=0 +Group.predef.swapUV=0 +Group.predef.negateU=0 +Group.predef.negateV=0 +Group.visible=1 +Group.suppress=0 +Group.relaxConstraints=0 +Group.allowRedundant=0 +Group.allDimsReference=0 +Group.remap={ +} +AddGroup + +Group.h.v=00000002 +Group.type=5001 +Group.order=1 +Group.name=sketch-in-plane +Group.activeWorkplane.v=80020000 +Group.color=ff000000 +Group.subtype=6000 +Group.skipFirst=0 +Group.predef.q.w=1.00000000000000000000 +Group.predef.origin.v=00010001 +Group.predef.swapUV=0 +Group.predef.negateU=0 +Group.predef.negateV=0 +Group.visible=1 +Group.suppress=0 +Group.relaxConstraints=0 +Group.allowRedundant=0 +Group.allDimsReference=0 +Group.remap={ +} +AddGroup + +Param.h.v.=00010010 +AddParam + +Param.h.v.=00010011 +AddParam + +Param.h.v.=00010012 +AddParam + +Param.h.v.=00010020 +Param.val=1.00000000000000000000 +AddParam + +Param.h.v.=00010021 +AddParam + +Param.h.v.=00010022 +AddParam + +Param.h.v.=00010023 +AddParam + +Param.h.v.=00020010 +AddParam + +Param.h.v.=00020011 +AddParam + +Param.h.v.=00020012 +AddParam + +Param.h.v.=00020020 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020021 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020022 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00020023 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00030010 +AddParam + +Param.h.v.=00030011 +AddParam + +Param.h.v.=00030012 +AddParam + +Param.h.v.=00030020 +Param.val=0.50000000000000000000 +AddParam + +Param.h.v.=00030021 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00030022 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00030023 +Param.val=-0.50000000000000000000 +AddParam + +Param.h.v.=00040010 +Param.val=-873.94336616523173688620 +AddParam + +Param.h.v.=00040011 +Param.val=888.35750558391430331540 +AddParam + +Param.h.v.=00040013 +Param.val=-310.42323046797423558019 +AddParam + +Param.h.v.=00040014 +Param.val=-3071.74919205236892594257 +AddParam + +Param.h.v.=00050010 +Param.val=-943.11295456966990968795 +AddParam + +Param.h.v.=00050011 +Param.val=-1154.79416602316041462473 +AddParam + +Param.h.v.=00050013 +Param.val=-4920.05203201180302130524 +AddParam + +Param.h.v.=00050014 +Param.val=-1583.69459073004077254154 +AddParam + +Param.h.v.=40000002 +Param.val=0.56702743975097524842 +AddParam + +Request.h.v=00000001 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000002 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000003 +Request.type=100 +Request.group.v=00000001 +Request.construction=0 +AddRequest + +Request.h.v=00000004 +Request.type=200 +Request.workplane.v=80020000 +Request.group.v=00000002 +Request.construction=0 +AddRequest + +Request.h.v=00000005 +Request.type=200 +Request.workplane.v=80020000 +Request.group.v=00000002 +Request.construction=0 +AddRequest + +Entity.h.v=00010000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00010001 +Entity.normal.v=00010020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00010001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00010020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00010001 +Entity.actNormal.w=1.00000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00020001 +Entity.normal.v=00020020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00020020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00020001 +Entity.actNormal.w=0.50000000000000000000 +Entity.actNormal.vx=0.50000000000000000000 +Entity.actNormal.vy=0.50000000000000000000 +Entity.actNormal.vz=0.50000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=00030001 +Entity.normal.v=00030020 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030001 +Entity.type=2000 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00030020 +Entity.type=3000 +Entity.construction=0 +Entity.point[0].v=00030001 +Entity.actNormal.w=0.50000000000000000000 +Entity.actNormal.vx=-0.50000000000000000000 +Entity.actNormal.vy=-0.50000000000000000000 +Entity.actNormal.vz=-0.50000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040000 +Entity.type=11000 +Entity.construction=0 +Entity.point[0].v=00040001 +Entity.point[1].v=00040002 +Entity.workplane.v=80020000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040001 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-873.94336616523173688620 +Entity.actPoint.y=888.35750558391430331540 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00040002 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-310.42323046797423558019 +Entity.actPoint.y=-3071.74919205236892594257 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050000 +Entity.type=11000 +Entity.construction=0 +Entity.point[0].v=00050001 +Entity.point[1].v=00050002 +Entity.workplane.v=80020000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050001 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-943.11295456966990968795 +Entity.actPoint.y=-1154.79416602316041462473 +Entity.actVisible=1 +AddEntity + +Entity.h.v=00050002 +Entity.type=2001 +Entity.construction=0 +Entity.workplane.v=80020000 +Entity.actPoint.x=-4920.05203201180302130524 +Entity.actPoint.y=-1583.69459073004077254154 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020000 +Entity.type=10000 +Entity.construction=0 +Entity.point[0].v=80020002 +Entity.normal.v=80020001 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020001 +Entity.type=3010 +Entity.construction=0 +Entity.point[0].v=80020002 +Entity.actNormal.w=1.00000000000000000000 +Entity.actVisible=1 +AddEntity + +Entity.h.v=80020002 +Entity.type=2012 +Entity.construction=1 +Entity.actVisible=1 +AddEntity + +Constraint.h.v=00000002 +Constraint.type=42 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valP.v=40000002 +Constraint.ptA.v=00050001 +Constraint.entityA.v=00040000 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +AddConstraint + +Constraint.h.v=00000003 +Constraint.type=30 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valA=4000.00000000000000000000 +Constraint.ptA.v=00040001 +Constraint.ptB.v=00040002 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +Constraint.disp.offset.x=497.85777038281713657852 +Constraint.disp.offset.y=-105.42002325418791031097 +AddConstraint + +Constraint.h.v=00000004 +Constraint.type=30 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.valA=4000.00000000000000000000 +Constraint.ptA.v=00050001 +Constraint.ptB.v=00050002 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +Constraint.disp.offset.x=-27.53621872585551599855 +Constraint.disp.offset.y=264.80569179149119918293 +AddConstraint + +Constraint.h.v=00000005 +Constraint.type=122 +Constraint.group.v=00000002 +Constraint.workplane.v=80020000 +Constraint.entityA.v=00050000 +Constraint.entityB.v=00040000 +Constraint.other=0 +Constraint.other2=0 +Constraint.reference=0 +AddConstraint + diff --git a/test/constraint/large_dimensions/test.cpp b/test/constraint/large_dimensions/test.cpp new file mode 100644 index 000000000..a242d728e --- /dev/null +++ b/test/constraint/large_dimensions/test.cpp @@ -0,0 +1,42 @@ +#include "solvespace.h" + +#include "harness.h" + +// The two line segments of the sketch, each dimensioned to 4 m, constrained +// perpendicular to each other with the start of the second one on the first. +static const hEntity LINE_A_A = { 0x00040001 }, LINE_A_B = { 0x00040002 }; +static const hEntity LINE_B_A = { 0x00050001 }, LINE_B_B = { 0x00050002 }; + +static Vector Dir(hEntity a, hEntity b) { + return SK.GetEntity(b)->PointGetNum().Minus(SK.GetEntity(a)->PointGetNum()); +} + +// A sketch whose dimensions are large enough that the entries of the Jacobian +// span several orders of magnitude. Solving the Newton step through the normal +// equations squared that spread, which pushed a perfectly good pivot under the +// threshold below which a column counts as linearly dependent; the step then +// came back with that component zeroed, and the solver reported these entirely +// compatible constraints as incompatible. See #1354. +TEST_CASE(perpendicular_4m) { + // Unlike most fixtures in this suite, perpendicular_4m.slvs is the + // reporter's file exactly as it was uploaded, and it must stay that way. + // Re-saving it from SolveSpace writes out the *solved* parameter values, + // so the sketch loads already at its solution, the first Newton step is + // the zero step, and the bug this test is here for cannot happen. Do not + // canonicalize this file. + CHECK_LOAD("perpendicular_4m.slvs"); + + for(Group &g : SK.group) { + CHECK_TRUE(g.IsSolvedOkay()); + } + + // And it solved it correctly, not just to something. + Vector da = Dir(LINE_A_A, LINE_A_B), db = Dir(LINE_B_A, LINE_B_B); + CHECK_EQ_EPS(da.Magnitude(), 4000.0); + CHECK_EQ_EPS(db.Magnitude(), 4000.0); + CHECK_EQ_EPS(da.WithMagnitude(1).Dot(db.WithMagnitude(1)), 0.0); + + // The start of the second line lies on the first one. + Vector onLine = Dir(LINE_A_A, LINE_B_A); + CHECK_EQ_EPS(onLine.Cross(da).Magnitude() / da.Magnitude(), 0.0); +}