diff --git a/src/lib.cpp b/src/lib.cpp index 0c27324..0a772e3 100644 --- a/src/lib.cpp +++ b/src/lib.cpp @@ -212,7 +212,8 @@ default: dbp("bad constraint type %d", sc->type); return; ssys->result = SLVS_RESULT_DIDNT_CONVERGE; break; - case System::REDUNDANT: + case System::REDUNDANT_DIDNT_CONVERGE: + case System::REDUNDANT_OKAY: ssys->result = SLVS_RESULT_INCONSISTENT; break; diff --git a/src/solvespace.h b/src/solvespace.h index 63d6df0..9e83109 100644 --- a/src/solvespace.h +++ b/src/solvespace.h @@ -307,6 +307,7 @@ public: static const double RANK_MAG_TOLERANCE, CONVERGE_TOLERANCE; int CalculateRank(void); + bool TestRank(void); static bool SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS], double B[], int N); bool SolveLeastSquares(void); @@ -323,10 +324,11 @@ public: bool NewtonSolve(int tag); enum { - SOLVED_OKAY = 0, - DIDNT_CONVERGE = 10, - REDUNDANT = 11, - TOO_MANY_UNKNOWNS = 20 + SOLVED_OKAY = 0, + DIDNT_CONVERGE = 10, + REDUNDANT_OKAY = 11, + REDUNDANT_DIDNT_CONVERGE = 12, + TOO_MANY_UNKNOWNS = 20 }; int Solve(Group *g, int *dof, List *bad, bool andFindBad, bool andFindFree); diff --git a/src/system.cpp b/src/system.cpp index 5a83249..f3d200c 100644 --- a/src/system.cpp +++ b/src/system.cpp @@ -173,6 +173,11 @@ int System::CalculateRank(void) { return rank; } +bool System::TestRank(void) { + EvalJacobian(); + return CalculateRank() == mat.m; +} + bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS], double B[], int n) { @@ -195,7 +200,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS], // Don't give up on a singular matrix unless it's really bad; the // assumption code is responsible for identifying that condition, // so we're not responsible for reporting that error. - if(ffabs(max) < 1e-20) return false; + if(ffabs(max) < 1e-20) continue; // Swap row imax with row i for(jp = 0; jp < n; jp++) { @@ -217,7 +222,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS], // We've put the matrix in upper triangular form, so at this point we // can solve by back-substitution. for(i = n - 1; i >= 0; i--) { - if(ffabs(A[i][i]) < 1e-20) return false; + if(ffabs(A[i][i]) < 1e-20) continue; temp = B[i]; for(j = n - 1; j > i; j--) { @@ -399,6 +404,8 @@ int System::Solve(Group *g, int *dof, List *bad, WriteEquationsExceptFor(Constraint::NO_CONSTRAINT, g); int i, j = 0; + bool rankOk; + /* dbp("%d equations", eq.n); for(i = 0; i < eq.n; i++) { @@ -447,20 +454,18 @@ int System::Solve(Group *g, int *dof, List *bad, return System::TOO_MANY_UNKNOWNS; } + rankOk = TestRank(); + // And do the leftovers as one big system if(!NewtonSolve(0)) { goto didnt_converge; } - EvalJacobian(); - - int rank; rank = CalculateRank(); - if(rank != mat.m) { - if(andFindBad) { - FindWhichToRemoveToFixJacobian(g, bad); - } - return System::REDUNDANT; + if(!TestRank()) { + if(andFindBad) FindWhichToRemoveToFixJacobian(g, bad); + return System::REDUNDANT_OKAY; } + // This is not the full Jacobian, but any substitutions or single-eq // solves removed one equation and one unknown, therefore no effect // on the number of DOF. @@ -478,7 +483,7 @@ int System::Solve(Group *g, int *dof, List *bad, p->tag = VAR_DOF_TEST; WriteJacobian(0); EvalJacobian(); - rank = CalculateRank(); + int rank = CalculateRank(); if(rank == mat.m) { p->free = true; } @@ -523,6 +528,6 @@ didnt_converge: } } - return System::DIDNT_CONVERGE; + return rankOk ? System::DIDNT_CONVERGE : System::REDUNDANT_DIDNT_CONVERGE; }