zcy
/
solvespace
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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]
solver
Jonathan Westhues 2008-06-26 01:34:26 -08:00
parent a47a77c37c
commit a31782e1ea
5 changed files with 113 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
expr.cpp
View File

@ -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;

3
expr.h
View File

@ -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, ...);

4
solvespace.h
View File

@ -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);

108
system.cpp
View File

@ -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(&param, &(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(&param);
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);
}

3
wishlist.txt
View File

@ -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