Skip to content

Commit

Permalink
moving L and U to state
Browse files Browse the repository at this point in the history
  • Loading branch information
@K20shores
K20shores committed Oct 20, 2023
1 parent ddeaaa3 commit f81d44d
Show file tree
Hide file tree
Showing 10 changed files with 725 additions and 670 deletions.
10 changes: 5 additions & 5 deletions include/micm/solver/linear_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,6 @@ namespace micm
std::vector<std::pair<std::size_t, std::size_t>> Uij_xj_;

LuDecompositionPolicy lu_decomp_;
SparseMatrixPolicy<T> lower_matrix_;
SparseMatrixPolicy<T> upper_matrix_;

public:
/// @brief default constructor
Expand All @@ -71,15 +69,17 @@ namespace micm
const std::function<LuDecompositionPolicy(const SparseMatrixPolicy<T>&)> create_lu_decomp);

/// @brief Decompose the matrix into upper and lower triangular matrices
void Factor(const SparseMatrixPolicy<T>& matrix);
void Factor(const SparseMatrixPolicy<T>& matrix, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix);

/// @brief Solve for x in Ax = b
template<template<class> class MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<T>> || !VectorizableSparse<SparseMatrixPolicy<T>>)
void Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x);
void Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix);
template<template<class> class MatrixPolicy>
requires(VectorizableDense<MatrixPolicy<T>> && VectorizableSparse<SparseMatrixPolicy<T>>)
void Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x);
void Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix);

std::pair<SparseMatrixPolicy<T>, SparseMatrixPolicy<T>> GetLUMatrices(const SparseMatrixPolicy<T>& matrix, T initial_value) const;
};

} // namespace micm
Expand Down
54 changes: 30 additions & 24 deletions include/micm/solver/linear_solver.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,55 +72,56 @@ namespace micm
lu_decomp_(create_lu_decomp(matrix))
{
auto lu = lu_decomp_.GetLUMatrices(matrix, initial_value);
lower_matrix_ = std::move(lu.first);
upper_matrix_ = std::move(lu.second);
for (std::size_t i = 0; i < lower_matrix_[0].size(); ++i)
auto lower_matrix = std::move(lu.first);
auto upper_matrix = std::move(lu.second);
for (std::size_t i = 0; i < lower_matrix[0].size(); ++i)
{
std::size_t nLij = 0;
for (std::size_t j_id = lower_matrix_.RowStartVector()[i]; j_id < lower_matrix_.RowStartVector()[i + 1]; ++j_id)
for (std::size_t j_id = lower_matrix.RowStartVector()[i]; j_id < lower_matrix.RowStartVector()[i + 1]; ++j_id)
{
std::size_t j = lower_matrix_.RowIdsVector()[j_id];
std::size_t j = lower_matrix.RowIdsVector()[j_id];
if (j >= i)
break;
Lij_yj_.push_back(std::make_pair(lower_matrix_.VectorIndex(0, i, j), j));
Lij_yj_.push_back(std::make_pair(lower_matrix.VectorIndex(0, i, j), j));
++nLij;
}
// There must always be a non-zero element on the diagonal
nLij_Lii_.push_back(std::make_pair(nLij, lower_matrix_.VectorIndex(0, i, i)));
nLij_Lii_.push_back(std::make_pair(nLij, lower_matrix.VectorIndex(0, i, i)));
}
for (std::size_t i = upper_matrix_[0].size() - 1; i != static_cast<std::size_t>(-1); --i)
for (std::size_t i = upper_matrix[0].size() - 1; i != static_cast<std::size_t>(-1); --i)
{
std::size_t nUij = 0;
for (std::size_t j_id = upper_matrix_.RowStartVector()[i]; j_id < upper_matrix_.RowStartVector()[i + 1]; ++j_id)
for (std::size_t j_id = upper_matrix.RowStartVector()[i]; j_id < upper_matrix.RowStartVector()[i + 1]; ++j_id)
{
std::size_t j = upper_matrix_.RowIdsVector()[j_id];
std::size_t j = upper_matrix.RowIdsVector()[j_id];
if (j <= i)
continue;
Uij_xj_.push_back(std::make_pair(upper_matrix_.VectorIndex(0, i, j), j));
Uij_xj_.push_back(std::make_pair(upper_matrix.VectorIndex(0, i, j), j));
++nUij;
}
// There must always be a non-zero element on the diagonal
nUij_Uii_.push_back(std::make_pair(nUij, upper_matrix_.VectorIndex(0, i, i)));
nUij_Uii_.push_back(std::make_pair(nUij, upper_matrix.VectorIndex(0, i, i)));
}
};


template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Factor(const SparseMatrixPolicy<T>& matrix)
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Factor(const SparseMatrixPolicy<T>& matrix, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix)
{
lu_decomp_.template Decompose<T, SparseMatrixPolicy>(matrix, lower_matrix_, upper_matrix_);
lu_decomp_.template Decompose<T, SparseMatrixPolicy>(matrix, lower_matrix, upper_matrix);
}

template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
template<template<class> class MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<T>> || !VectorizableSparse<SparseMatrixPolicy<T>>)
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x)
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix)
{
for (std::size_t i_cell = 0; i_cell < b.size(); ++i_cell)
{
auto b_cell = b[i_cell];
auto x_cell = x[i_cell];
const std::size_t lower_grid_offset = i_cell * lower_matrix_.FlatBlockSize();
const std::size_t upper_grid_offset = i_cell * upper_matrix_.FlatBlockSize();
const std::size_t lower_grid_offset = i_cell * lower_matrix.FlatBlockSize();
const std::size_t upper_grid_offset = i_cell * upper_matrix.FlatBlockSize();
auto& y_cell = x_cell; // Alias x for consistency with equations, but to reuse memory
{
auto b_elem = b_cell.begin();
Expand All @@ -131,10 +132,10 @@ namespace micm
*y_elem = *(b_elem++);
for (std::size_t i = 0; i < nLij_Lii.first; ++i)
{
*y_elem -= lower_matrix_.AsVector()[lower_grid_offset + (*Lij_yj).first] * y_cell[(*Lij_yj).second];
*y_elem -= lower_matrix.AsVector()[lower_grid_offset + (*Lij_yj).first] * y_cell[(*Lij_yj).second];
++Lij_yj;
}
*(y_elem++) /= lower_matrix_.AsVector()[lower_grid_offset + nLij_Lii.second];
*(y_elem++) /= lower_matrix.AsVector()[lower_grid_offset + nLij_Lii.second];
}
}
{
Expand All @@ -151,12 +152,12 @@ namespace micm

for (std::size_t i = 0; i < nUij_Uii.first; ++i)
{
*x_elem -= upper_matrix_.AsVector()[upper_grid_offset + (*Uij_xj).first] * x_cell[(*Uij_xj).second];
*x_elem -= upper_matrix.AsVector()[upper_grid_offset + (*Uij_xj).first] * x_cell[(*Uij_xj).second];
++Uij_xj;
}

// don't iterate before the beginning of the vector
*(x_elem) /= upper_matrix_.AsVector()[upper_grid_offset + nUij_Uii.second];
*(x_elem) /= upper_matrix.AsVector()[upper_grid_offset + nUij_Uii.second];
if (x_elem != x_cell.begin())
{
--x_elem;
Expand All @@ -169,7 +170,7 @@ namespace micm
template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
template<template<class> class MatrixPolicy>
requires(VectorizableDense<MatrixPolicy<T>> && VectorizableSparse<SparseMatrixPolicy<T>>)
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x)
inline void LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::Solve(const MatrixPolicy<T>& b, MatrixPolicy<T>& x, SparseMatrixPolicy<T>& lower_matrix, SparseMatrixPolicy<T>& upper_matrix)
{
const std::size_t n_cells = b.GroupVectorSize();
// Loop over groups of blocks
Expand All @@ -178,9 +179,9 @@ namespace micm
auto b_group = std::next(b.AsVector().begin(), i_group * b.GroupSize());
auto x_group = std::next(x.AsVector().begin(), i_group * x.GroupSize());
auto L_group =
std::next(lower_matrix_.AsVector().begin(), i_group * lower_matrix_.GroupSize(lower_matrix_.FlatBlockSize()));
std::next(lower_matrix.AsVector().begin(), i_group * lower_matrix.GroupSize(lower_matrix.FlatBlockSize()));
auto U_group =
std::next(upper_matrix_.AsVector().begin(), i_group * upper_matrix_.GroupSize(upper_matrix_.FlatBlockSize()));
std::next(upper_matrix.AsVector().begin(), i_group * upper_matrix.GroupSize(upper_matrix.FlatBlockSize()));
auto y_group = x_group; // Alias x for consistency with equations, but to reuse memory
{
auto b_elem = b_group;
Expand Down Expand Up @@ -229,4 +230,9 @@ namespace micm
}
}

template<typename T, template<class> class SparseMatrixPolicy, class LuDecompositionPolicy>
std::pair<SparseMatrixPolicy<T>, SparseMatrixPolicy<T>> LinearSolver<T, SparseMatrixPolicy, LuDecompositionPolicy>::GetLUMatrices(const SparseMatrixPolicy<T>& matrix, T initial_value) const {
return lu_decomp_.GetLUMatrices(matrix, initial_value);
}

} // namespace micm
2 changes: 1 addition & 1 deletion include/micm/solver/rosenbrock.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,7 @@ namespace micm
/// @param singular indicates if the matrix is singular
/// @param number_densities constituent concentration (molec/cm^3)
/// @param rate_constants Rate constants for each process (molecule/cm3)^(n-1) s-1
void LinearFactor(double& H, const double gamma, bool& singular, const MatrixPolicy<double>& number_densities, SolverStats& stats, SparseMatrixPolicy<double> jacobian);
void LinearFactor(double& H, const double gamma, bool& singular, const MatrixPolicy<double>& number_densities, SolverStats& stats, State<MatrixPolicy, SparseMatrixPolicy>& state);

protected:
/// @brief Computes the scaled norm of the vector errors
Expand Down
16 changes: 12 additions & 4 deletions include/micm/solver/rosenbrock.inl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -536,12 +536,19 @@ namespace micm
inline State<MatrixPolicy, SparseMatrixPolicy> RosenbrockSolver<MatrixPolicy, SparseMatrixPolicy, LinearSolverPolicy>::GetState() const
{
auto state = State<MatrixPolicy, SparseMatrixPolicy>{ state_parameters_ };

state.jacobian_ = build_jacobian<SparseMatrixPolicy>(
state_parameters_.nonzero_jacobian_elements_,
state_parameters_.number_of_grid_cells_,
system_.StateSize()
);

auto lu = linear_solver_.GetLUMatrices(state.jacobian_, 1.0e-30);
auto lower_matrix = std::move(lu.first);
auto upper_matrix = std::move(lu.second);
state.lower_matrix_ = lower_matrix;
state.upper_matrix_ = upper_matrix;

return state;
}

Expand Down Expand Up @@ -611,7 +618,7 @@ namespace micm
{
bool is_singular{ false };
// Form and factor the rosenbrock ode jacobian
TIMED_METHOD(stats.total_linear_factor_time, time_it, LinearFactor, H, parameters_.gamma_[0], is_singular, Y, stats, state.jacobian_);
TIMED_METHOD(stats.total_linear_factor_time, time_it, LinearFactor, H, parameters_.gamma_[0], is_singular, Y, stats, state);
if (is_singular)
{
result.state_ = SolverState::RepeatedlySingularMatrix;
Expand Down Expand Up @@ -647,7 +654,7 @@ namespace micm
K[stage].ForEach([&](double& iKstage, double& iKj) { iKstage += HC * iKj; }, K[j]);
}
temp.AsVector().assign(K[stage].AsVector().begin(), K[stage].AsVector().end());
TIMED_METHOD(stats.total_linear_solve_time, time_it, linear_solver_.template Solve<MatrixPolicy>, temp, K[stage]);
TIMED_METHOD(stats.total_linear_solve_time, time_it, linear_solver_.template Solve<MatrixPolicy>, temp, K[stage], state.lower_matrix_, state.upper_matrix_);
stats.solves += 1;
}

Expand Down Expand Up @@ -792,17 +799,18 @@ namespace micm
bool& singular,
const MatrixPolicy<double>& number_densities,
SolverStats& stats,
SparseMatrixPolicy<double> jacobian)
State<MatrixPolicy, SparseMatrixPolicy>& state)
{
// TODO: invesitage this function. The fortran equivalent appears to have a bug.
// From my understanding the fortran do loop would only ever do one iteration and is equivalent to what's below
auto jacobian = state.jacobian_;
uint64_t n_consecutive = 0;
singular = true;
while (true)
{
double alpha = 1 / (H * gamma);
AlphaMinusJacobian(jacobian, alpha);
linear_solver_.Factor(jacobian);
linear_solver_.Factor(jacobian, state.lower_matrix_, state.upper_matrix_);
singular = false; // TODO This should be evaluated in some way
stats.decompositions += 1;
if (!singular)
Expand Down
2 changes: 2 additions & 0 deletions include/micm/solver/state.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,8 @@ namespace micm
std::map<std::string, std::size_t> variable_map_;
std::map<std::string, std::size_t> custom_rate_parameter_map_;
std::vector<std::string> variable_names_{};
SparseMatrixPolicy<double> lower_matrix_;
SparseMatrixPolicy<double> upper_matrix_;

/// @brief
State();
Expand Down
Loading

0 comments on commit f81d44d

Please sign in to comment.