12 #include <Eigen/Sparse>
57 LinearSolver(
int num_rows_,
int num_variables_,
int num_rhs_,
bool lsq_)
59 assert(num_variables_ > 0);
60 assert(num_rhs_ <= 4);
87 std::vector<EigenVectorX>
b;
88 std::vector<EigenVectorX>
x;
103 return new LinearSolver(num_rows, num_columns, num_rhs,
false);
108 return new LinearSolver(num_rows, num_columns, num_rhs,
true);
120 solver->
variable[index].value[rhs] = value;
125 return solver->
variable[index].value[rhs];
130 if (!solver->
variable[index].locked) {
132 solver->
variable[index].locked =
true;
138 if (solver->
variable[index].locked) {
140 solver->
variable[index].locked =
false;
151 for (
int j = 0; j < num_rhs; j++)
152 solver->
x[j][
v->index] =
v->value[j];
164 for (
int j = 0; j < num_rhs; j++)
165 v->value[j] = solver->
x[j][
v->index];
199 for (
int i = 0; i < solver->
num_rhs; i++) {
200 solver->
b[i].setZero(m);
201 solver->
x[i].setZero(n);
246 solver->
b[rhs][index] += value;
248 else if (!solver->
variable[index].locked) {
249 index = solver->
variable[index].index;
250 solver->
b[rhs][index] += value;
259 if (solver->
m == 0 || solver->
n == 0)
268 solver->
M.resize(solver->
m, solver->
n);
274 solver->
MtM = solver->
M.transpose() * solver->
M;
284 sparseLU->compute(
M);
285 result = (sparseLU->info() == Eigen::Success);
292 for (
int rhs = 0; rhs < solver->
num_rhs; rhs++) {
300 std::vector<LinearSolver::Coeff> &
a = variable->
a;
302 for (
int j = 0; j <
a.size(); j++)
303 b[
a[j].index] -=
a[j].value * variable->
value[rhs];
310 solver->
x[rhs] = solver->
sparseLU->solve(Mtb);
317 if (solver->
sparseLU->info() != Eigen::Success)
326 for (
int rhs = 0; rhs < solver->
num_rhs; rhs++)
327 solver->
b[rhs].setZero(solver->
m);
336 std::cout <<
"A:" << solver->
M << std::endl;
338 for (
int rhs = 0; rhs < solver->
num_rhs; rhs++)
339 std::cout <<
"b " << rhs <<
":" << solver->
b[rhs] << std::endl;
341 if (solver->
MtM.rows() && solver->
MtM.cols())
342 std::cout <<
"AtA:" << solver->
MtM << std::endl;
ATTR_WARN_UNUSED_RESULT const BMVert * v
void EIG_linear_solver_print_matrix(LinearSolver *solver)
Eigen::VectorXd EigenVectorX
void EIG_linear_solver_variable_set(LinearSolver *solver, int rhs, int index, double value)
LinearSolver * EIG_linear_solver_new(int num_rows, int num_columns, int num_rhs)
void EIG_linear_solver_right_hand_side_add(LinearSolver *solver, int rhs, int index, double value)
void EIG_linear_solver_variable_unlock(LinearSolver *solver, int index)
Eigen::SparseMatrix< double, Eigen::ColMajor > EigenSparseMatrix
Eigen::SparseLU< EigenSparseMatrix > EigenSparseLU
LinearSolver * EIG_linear_least_squares_solver_new(int num_rows, int num_columns, int num_rhs)
void EIG_linear_solver_delete(LinearSolver *solver)
Eigen::Triplet< double > EigenTriplet
static void linear_solver_variables_to_vector(LinearSolver *solver)
double EIG_linear_solver_variable_get(LinearSolver *solver, int rhs, int index)
static void linear_solver_vector_to_variables(LinearSolver *solver)
void EIG_linear_solver_matrix_add(LinearSolver *solver, int row, int col, double value)
bool EIG_linear_solver_solve(LinearSolver *solver)
static void linear_solver_ensure_matrix_construct(LinearSolver *solver)
void EIG_linear_solver_variable_lock(LinearSolver *solver, int index)
struct LinearSolver LinearSolver
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
@ STATE_VARIABLES_CONSTRUCT
std::vector< Variable > variable
LinearSolver(int num_rows_, int num_variables_, int num_rhs_, bool lsq_)
std::vector< EigenVectorX > b
std::vector< EigenVectorX > x
std::vector< EigenTriplet > Mtriplets