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

Don't hold the sparse matrices in a pointer at all.

pull/1170/head
Ryan Pavlik 2021-12-23 11:37:33 -06:00 committed by phkahler
parent 42f5d3ab0d
commit 6edeb66e3d
2 changed files with 22 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/solvespace.h
View File

@ -246,8 +246,8 @@ public:
int m, n; int m, n;
struct { struct {
// This only observes the Expr - does not own them! // This only observes the Expr - does not own them!
std::unique_ptr<Eigen::SparseMatrix<Expr *>> sym; Eigen::SparseMatrix<Expr *> sym;
std::unique_ptr<Eigen::SparseMatrix<double>> num; Eigen::SparseMatrix<double> num;
} A; } A;
Eigen::VectorXd scale; Eigen::VectorXd scale;

40
src/system.cpp
View File

@ -8,7 +8,6 @@
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
#include "solvespace.h" #include "solvespace.h"
#include <memory>
#include <Eigen/Core> #include <Eigen/Core>
#include <Eigen/SparseQR> #include <Eigen/SparseQR>
@ -20,6 +19,7 @@ bool System::WriteJacobian(int tag) {
// Clear all // Clear all
mat.param.clear(); mat.param.clear();
mat.eq.clear(); mat.eq.clear();
mat.A.sym.setZero();
mat.B.sym.clear(); mat.B.sym.clear();
for(Param &p : param) { for(Param &p : param) {
@ -33,8 +33,8 @@ bool System::WriteJacobian(int tag) {
mat.eq.push_back(&e); mat.eq.push_back(&e);
} }
mat.m = mat.eq.size(); mat.m = mat.eq.size();
mat.A.sym.reset(new Eigen::SparseMatrix<Expr *>(mat.m, mat.n)); mat.A.sym.resize(mat.m, mat.n);
mat.A.sym->reserve(Eigen::VectorXi::Constant(mat.n, 10)); mat.A.sym.reserve(Eigen::VectorXi::Constant(mat.n, 10));
// Fill the param id to index map // Fill the param id to index map
std::map<uint32_t, int> paramToIndex; std::map<uint32_t, int> paramToIndex;
@ -64,7 +64,7 @@ bool System::WriteJacobian(int tag) {
pd = pd->FoldConstants(); pd = pd->FoldConstants();
if(pd->IsZeroConst()) if(pd->IsZeroConst())
continue; continue;
mat.A.sym->insert(i, j->second) = pd; mat.A.sym.insert(i, j->second) = pd;
} }
paramsUsed.clear(); paramsUsed.clear();
mat.B.sym.push_back(f); mat.B.sym.push_back(f);
@ -74,25 +74,25 @@ bool System::WriteJacobian(int tag) {
void System::EvalJacobian() { void System::EvalJacobian() {
using namespace Eigen; using namespace Eigen;
mat.A.num.reset(new Eigen::SparseMatrix<double>(mat.m, mat.n)); mat.A.num.setZero();
int size = mat.A.sym->outerSize(); mat.A.num.resize(mat.m, mat.n);
const int size = mat.A.sym.outerSize();
for(int k = 0; k < size; k++) { for(int k = 0; k < size; k++) {
for(SparseMatrix <Expr *>::InnerIterator it(*mat.A.sym, k); it; ++it) { for(SparseMatrix <Expr *>::InnerIterator it(mat.A.sym, k); it; ++it) {
double value = it.value()->Eval(); double value = it.value()->Eval();
if(EXACT(value == 0.0)) continue; if(EXACT(value == 0.0)) continue;
mat.A.num->insert(it.row(), it.col()) = value; mat.A.num.insert(it.row(), it.col()) = value;
} }
} }
mat.A.num->makeCompressed(); mat.A.num.makeCompressed();
} }
bool System::IsDragged(hParam p) { bool System::IsDragged(hParam p) {
hParam *pp; hParam *pp;
for(pp = dragged.First(); pp; pp = dragged.NextAfter(pp)) { const auto b = dragged.begin();
if(p == *pp) return true; const auto e = dragged.end();
} return e != std::find(b, e, p);
return false;
} }
Param *System::GetLastParamSubstitution(Param *p) { Param *System::GetLastParamSubstitution(Param *p) {
@ -217,7 +217,7 @@ int System::CalculateRank() {
using namespace Eigen; using namespace Eigen;
if(mat.n == 0 || mat.m == 0) return 0; if(mat.n == 0 || mat.m == 0) return 0;
SparseQR <SparseMatrix<double>, COLAMDOrdering<int>> solver; SparseQR <SparseMatrix<double>, COLAMDOrdering<int>> solver;
solver.compute(*mat.A.num); solver.compute(mat.A.num);
int result = solver.rank(); int result = solver.rank();
return result; return result;
} }
@ -259,20 +259,20 @@ bool System::SolveLeastSquares() {
} }
} }
int size = mat.A.sym->outerSize(); int size = mat.A.sym.outerSize();
for(int k = 0; k < size; k++) { for(int k = 0; k < size; k++) {
for(SparseMatrix<double>::InnerIterator it(*mat.A.num, k); it; ++it) { for(SparseMatrix<double>::InnerIterator it(mat.A.num, k); it; ++it) {
it.valueRef() *= mat.scale[it.col()]; it.valueRef() *= mat.scale[it.col()];
} }
} }
SparseMatrix <double> AAt = *mat.A.num * mat.A.num->transpose(); SparseMatrix<double> AAt = mat.A.num * mat.A.num.transpose();
AAt.makeCompressed(); AAt.makeCompressed();
VectorXd z(mat.n); VectorXd z(mat.n);
if(!SolveLinearSystem(AAt, mat.B.num, &z)) return false; 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++) { for(int c = 0; c < mat.n; c++) {
mat.X[c] *= mat.scale[c]; mat.X[c] *= mat.scale[c];
@ -560,8 +560,8 @@ void System::Clear() {
param.Clear(); param.Clear();
eq.Clear(); eq.Clear();
dragged.Clear(); dragged.Clear();
mat.A.num.reset(); mat.A.num.setZero();
mat.A.sym.reset(); mat.A.sym.setZero();
} }
void System::MarkParamsFree(bool find) { void System::MarkParamsFree(bool find) {