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

Little optimizations in the solver; don't write the Jacobian 

redundantly, more zero partial avoidance, slight speedup to linear
system solve.

[git-p4: depot-paths = "//depot/solvespace/": change = 1940]
solver
Jonathan Westhues 2009-04-18 19:55:46 -08:00
parent 16de9a485a
commit bab13b821f
4 changed files with 25 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
expr.cpp
View File

@ -364,10 +364,10 @@ Expr *Expr::PartialWrt(hParam p) {
} }
} }
DWORD Expr::ParamsUsed(void) { QWORD Expr::ParamsUsed(void) {
DWORD r = 0; QWORD r = 0;
if(op == PARAM) r |= (1 << (x.parh.v % 31)); if(op == PARAM) r |= ((QWORD)1 << (x.parh.v % 61));
if(op == PARAM_PTR) r |= (1 << (x.parp->h.v % 31)); if(op == PARAM_PTR) r |= ((QWORD)1 << (x.parp->h.v % 61));
int c = Children(); int c = Children();
if(c >= 1) r |= a->ParamsUsed(); if(c >= 1) r |= a->ParamsUsed();
@ -375,6 +375,16 @@ DWORD Expr::ParamsUsed(void) {
return r; return r;
} }
bool Expr::DependsOn(hParam p) {
if(op == PARAM) return (x.parh.v == p.v);
if(op == PARAM_PTR) return (x.parp->h.v == p.v);
int c = Children();
if(c == 1) return a->DependsOn(p);
if(c == 2) return a->DependsOn(p) || b->DependsOn(p);
return false;
}
bool Expr::Tol(double a, double b) { bool Expr::Tol(double a, double b) {
return fabs(a - b) < 0.001; return fabs(a - b) < 0.001;
} }

3
expr.h
View File

@ -77,7 +77,8 @@ public:
Expr *PartialWrt(hParam p); Expr *PartialWrt(hParam p);
double Eval(void); double Eval(void);
DWORD ParamsUsed(void); QWORD ParamsUsed(void);
bool DependsOn(hParam p);
static bool Tol(double a, double b); static bool Tol(double a, double b);
Expr *FoldConstants(void); Expr *FoldConstants(void);
void Substitute(hParam oldh, hParam newh); void Substitute(hParam oldh, hParam newh);

3
solvespace.h
View File

@ -43,6 +43,7 @@ inline double WRAP_SYMMETRIC(double v, double n) {
#define isforname(c) (isalnum(c) || (c) == '_' || (c) == '-' || (c) == '#') #define isforname(c) (isalnum(c) || (c) == '_' || (c) == '-' || (c) == '#')
typedef unsigned __int64 QWORD;
typedef signed long SDWORD; typedef signed long SDWORD;
typedef signed short SWORD; typedef signed short SWORD;
@ -201,7 +202,7 @@ void MakePathAbsolute(char *base, char *path);
class System { class System {
public: public:
#define MAX_UNKNOWNS 200 #define MAX_UNKNOWNS 1000
EntityList entity; EntityList entity;
ParamList param; ParamList param;

12
system.cpp
View File

@ -24,11 +24,13 @@ void System::WriteJacobian(int tag) {
Expr *f = e->e->DeepCopyWithParamsAsPointers(&param, &(SS.param)); Expr *f = e->e->DeepCopyWithParamsAsPointers(&param, &(SS.param));
f = f->FoldConstants(); f = f->FoldConstants();
// Hash table (31 bits) to accelerate generation of zero partials. // Hash table (61 bits) to accelerate generation of zero partials.
DWORD scoreboard = f->ParamsUsed(); QWORD scoreboard = f->ParamsUsed();
for(j = 0; j < mat.n; j++) { for(j = 0; j < mat.n; j++) {
Expr *pd; Expr *pd;
if(scoreboard & (1 << (mat.param[j].v % 31))) { if(scoreboard & ((QWORD)1 << (mat.param[j].v % 61)) &&
f->DependsOn(mat.param[j]))
{
pd = f->PartialWrt(mat.param[j]); pd = f->PartialWrt(mat.param[j]);
pd = pd->FoldConstants(); pd = pd->FoldConstants();
pd = pd->DeepCopyWithParamsAsPointers(&param, &(SS.param)); pd = pd->DeepCopyWithParamsAsPointers(&param, &(SS.param));
@ -225,7 +227,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
for(ip = i+1; ip < n; ip++) { for(ip = i+1; ip < n; ip++) {
temp = A[ip][i]/A[i][i]; temp = A[ip][i]/A[i][i];
for(jp = 0; jp < n; jp++) { for(jp = i; jp < n; jp++) {
A[ip][jp] -= temp*(A[i][jp]); A[ip][jp] -= temp*(A[i][jp]);
} }
B[ip] -= temp*B[i]; B[ip] -= temp*B[i];
@ -291,7 +293,6 @@ bool System::SolveLeastSquares(void) {
} }
bool System::NewtonSolve(int tag) { bool System::NewtonSolve(int tag) {
WriteJacobian(tag);
if(mat.m > mat.n) return false; if(mat.m > mat.n) return false;
int iter = 0; int iter = 0;
@ -438,6 +439,7 @@ void System::Solve(Group *g, bool andFindFree) {
e->tag = alone; e->tag = alone;
p->tag = alone; p->tag = alone;
WriteJacobian(alone);
if(!NewtonSolve(alone)) { if(!NewtonSolve(alone)) {
// Failed to converge, bail out early // Failed to converge, bail out early
goto didnt_converge; goto didnt_converge;