From a31782e1ea5e18bf81e3dc135452d0bfe25c35e6 Mon Sep 17 00:00:00 2001 From: Jonathan Westhues Date: Thu, 26 Jun 2008 01:34:26 -0800 Subject: [PATCH] If any equations are soluble in a single variable, then handle them separately, even before doing the rank test. In some cases, that's a huge speedup. [git-p4: depot-paths = "//depot/solvespace/": change = 1811] --- expr.cpp | 49 +++++++++++++++++++++++ expr.h | 3 ++ solvespace.h | 4 +- system.cpp | 108 ++++++++++++++++++++++++++++----------------------- wishlist.txt | 3 +- 5 files changed, 113 insertions(+), 54 deletions(-) diff --git a/expr.cpp b/expr.cpp index 86a4cd9..817aac3 100644 --- a/expr.cpp +++ b/expr.cpp @@ -444,7 +444,48 @@ void Expr::Substitute(hParam oldh, hParam newh) { if(c >= 2) b->Substitute(oldh, newh); } +//----------------------------------------------------------------------------- +// If the expression references only one parameter that appears in pl, then +// return that parameter. If no param is referenced, then return NO_PARAMS. +// If multiple params are referenced, then return MULTIPLE_PARAMS. +//----------------------------------------------------------------------------- +const hParam Expr::NO_PARAMS = { 0 }; +const hParam Expr::MULTIPLE_PARAMS = { 1 }; +hParam Expr::ReferencedParams(ParamList *pl) { + if(op == PARAM) { + if(pl->FindByIdNoOops(x.parh)) { + return x.parh; + } else { + return NO_PARAMS; + } + } + if(op == PARAM_PTR) oops(); + int c = Children(); + if(c == 0) { + return NO_PARAMS; + } else if(c == 1) { + return a->ReferencedParams(pl); + } else if(c == 2) { + hParam pa, pb; + pa = a->ReferencedParams(pl); + pb = b->ReferencedParams(pl); + if(pa.v == NO_PARAMS.v) { + return pb; + } else if(pb.v == NO_PARAMS.v) { + return pa; + } else if(pa.v == pb.v) { + return pa; // either, doesn't matter + } else { + return MULTIPLE_PARAMS; + } + } else oops(); +} + + +//----------------------------------------------------------------------------- +// Routines to pretty-print an expression. Mostly for debugging. +//----------------------------------------------------------------------------- static char StringBuffer[4096]; void Expr::App(char *s, ...) { @@ -490,6 +531,14 @@ p: } } + +//----------------------------------------------------------------------------- +// A parser; convert a string to an expression. Infix notation, with the +// usual shift/reduce approach. I had great hopes for user-entered eq +// constraints, but those don't seem very useful, so right now this is just +// to provide calculator type functionality wherever numbers are entered. +//----------------------------------------------------------------------------- + #define MAX_UNPARSED 1024 static Expr *Unparsed[MAX_UNPARSED]; static int UnparsedCnt, UnparsedP; diff --git a/expr.h b/expr.h index d555d57..b4775b2 100644 --- a/expr.h +++ b/expr.h @@ -78,6 +78,9 @@ public: Expr *FoldConstants(void); void Substitute(hParam oldh, hParam newh); + static const hParam NO_PARAMS, MULTIPLE_PARAMS; + hParam ReferencedParams(ParamList *pl); + void ParamsToPointers(void); void App(char *str, ...); diff --git a/solvespace.h b/solvespace.h index 5221661..ef2dc36 100644 --- a/solvespace.h +++ b/solvespace.h @@ -168,8 +168,6 @@ public: // The corresponding parameter for each column hParam param[MAX_UNKNOWNS]; - bool bound[MAX_UNKNOWNS]; - // We're solving AX = B int m, n; struct { @@ -197,7 +195,7 @@ public: double B[], int N); bool SolveLeastSquares(void); - void WriteJacobian(int eqTag, int paramTag); + void WriteJacobian(int tag); void EvalJacobian(void); void WriteEquationsExceptFor(hConstraint hc, hGroup hg); diff --git a/system.cpp b/system.cpp index ca1e6bd..19c552f 100644 --- a/system.cpp +++ b/system.cpp @@ -1,12 +1,12 @@ #include "solvespace.h" -void System::WriteJacobian(int eqTag, int paramTag) { +void System::WriteJacobian(int tag) { int a, i, j; j = 0; for(a = 0; a < param.n; a++) { Param *p = &(param.elem[a]); - if(p->tag != paramTag) continue; + if(p->tag != tag) continue; mat.param[j] = p->h; j++; } @@ -15,7 +15,7 @@ void System::WriteJacobian(int eqTag, int paramTag) { i = 0; for(a = 0; a < eq.n; a++) { Equation *e = &(eq.elem[a]); - if(e->tag != eqTag) continue; + if(e->tag != tag) continue; mat.eq[i] = e->h; Expr *f = e->e->DeepCopyWithParamsAsPointers(¶m, &(SS.param)); @@ -152,10 +152,6 @@ bool System::Tol(double v) { int System::GaussJordan(void) { int i, j; - for(j = 0; j < mat.n; j++) { - mat.bound[j] = false; - } - // Now eliminate. i = 0; int rank = 0; @@ -201,8 +197,7 @@ int System::GaussJordan(void) { mat.A.num[is][j] = 0; } - // And mark this as a bound variable. - mat.bound[j] = true; + // And this is a bound variable, so rank goes up by one. rank++; // Move on to the next row, since we just used this one to @@ -314,7 +309,7 @@ bool System::SolveLeastSquares(void) { } bool System::NewtonSolve(int tag) { - WriteJacobian(tag, tag); + WriteJacobian(tag); if(mat.m > mat.n) return false; int iter = 0; @@ -396,7 +391,7 @@ void System::FindWhichToRemoveToFixJacobian(Group *g) { WriteEquationsExceptFor(c->h, g->h); eq.ClearTags(); - WriteJacobian(0, 0); + WriteJacobian(0); EvalJacobian(); int rank = GaussJordan(); @@ -423,26 +418,43 @@ void System::Solve(Group *g) { dbp(" param %08x at %.3f", param.elem[i].h.v, param.elem[i].val); } */ + // All params and equations are assigned to group zero. param.ClearTags(); eq.ClearTags(); SolveBySubstitution(); - WriteJacobian(0, 0); + // Before solving the big system, see if we can find any equations that + // are soluble alone. This can be a huge speedup. We don't know whether + // the system is consistent yet, but if it isn't then we'll catch that + // later. + int alone = 1; + for(i = 0; i < eq.n; i++) { + Equation *e = &(eq.elem[i]); + if(e->tag != 0) continue; + + hParam hp = e->e->ReferencedParams(¶m); + if(hp.v == Expr::NO_PARAMS.v) continue; + if(hp.v == Expr::MULTIPLE_PARAMS.v) continue; + + Param *p = param.FindById(hp); + if(p->tag != 0) continue; // let rank test catch inconsistency + + e->tag = alone; + p->tag = alone; + if(!NewtonSolve(alone)) { + // Failed to converge, bail out early + goto didnt_converge; + } + alone++; + } + + // Now write the Jacobian for what's left, and do a rank test; that + // tells us if the system is inconsistently constrained. + WriteJacobian(0); EvalJacobian(); -/* - for(i = 0; i < mat.m; i++) { - dbp("function %d: %s", i, mat.B.sym[i]->Print()); - } - dbp("m=%d", mat.m); - for(i = 0; i < mat.m; i++) { - for(j = 0; j < mat.n; j++) { - dbp("A(%d,%d) = %.10f;", i+1, j+1, mat.A.num[i][j]); - } - } */ - int rank = GaussJordan(); if(rank != mat.m) { FindWhichToRemoveToFixJacobian(g); @@ -451,35 +463,33 @@ void System::Solve(Group *g) { return; } -/* dbp("bound states:"); - for(j = 0; j < mat.n; j++) { - dbp(" param %08x: %d", mat.param[j], mat.bound[j]); - } */ + // And do the leftovers as one big system + if(!NewtonSolve(0)) { + goto didnt_converge; + } - bool ok = NewtonSolve(0); - - if(ok) { - // System solved correctly, so write the new values back in to the - // main parameter table. - for(i = 0; i < param.n; i++) { - Param *p = &(param.elem[i]); - double val; - if(p->tag == VAR_SUBSTITUTED) { - val = param.FindById(p->substd)->val; - } else { - val = p->val; - } - Param *pp = SS.GetParam(p->h); - pp->val = val; - pp->known = true; + // System solved correctly, so write the new values back in to the + // main parameter table. + for(i = 0; i < param.n; i++) { + Param *p = &(param.elem[i]); + double val; + if(p->tag == VAR_SUBSTITUTED) { + val = param.FindById(p->substd)->val; + } else { + val = p->val; } - if(g->solved.how != Group::SOLVED_OKAY) { - g->solved.how = Group::SOLVED_OKAY; - TextWindow::ReportHowGroupSolved(g->h); - } - } else { - g->solved.how = Group::DIDNT_CONVERGE; + Param *pp = SS.GetParam(p->h); + pp->val = val; + pp->known = true; + } + if(g->solved.how != Group::SOLVED_OKAY) { + g->solved.how = Group::SOLVED_OKAY; TextWindow::ReportHowGroupSolved(g->h); } + return; + +didnt_converge: + g->solved.how = Group::DIDNT_CONVERGE; + TextWindow::ReportHowGroupSolved(g->h); } diff --git a/wishlist.txt b/wishlist.txt index 5618416..9ae75d4 100644 --- a/wishlist.txt +++ b/wishlist.txt @@ -9,10 +9,9 @@ auto-generate circles and faces when lathing copy the section geometry to other end when sweeping cylindrical faces -partitioned subsystems in the solver - long term ----- incremental regen of entities? +my own plane poly triangulation code