Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79595529
sparse_solver_eigen.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Aug 27, 06:53
Size
2 KB
Mime Type
text/x-c
Expires
Thu, Aug 29, 06:53 (2 d)
Engine
blob
Format
Raw Data
Handle
20181951
Attached To
rAKA akantu
sparse_solver_eigen.cc
View Options
#include "sparse_solver_eigen.hh"
#include "dof_manager_default.hh"
#include "solver_vector_default.hh"
#include "sparse_matrix_aij.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
SparseSolverEigen::SparseSolverEigen(DOFManagerDefault & dof_manager,
const ID & matrix_id, const ID & id)
: SparseSolver(dof_manager, matrix_id, id), dof_manager(dof_manager) {
AKANTU_DEBUG_ASSERT(communicator.getNbProc() == 1,
"This solver does not work in parallel");
}
/* -------------------------------------------------------------------------- */
void SparseSolverEigen::setA() {
const auto & Aij =
aka::as_type<SparseMatrixAIJ>(dof_manager.getMatrix(matrix_id));
if (this->last_value_release == Aij.getRelease() and
this->last_value_release != -1) {
return;
}
using Triplet = Eigen::Triplet<Real>;
std::vector<Triplet> triplets;
triplets.reserve(Aij.getNbNonZero());
for (auto && [i, j, a] : zip(Aij.getIRN(), Aij.getJCN(), Aij.getA())) {
triplets.emplace_back(i - 1, j - 1, a);
if (Aij.getMatrixType() == _symmetric and i != j) {
triplets.emplace_back(j - 1, i - 1, a);
}
}
A.resize(Aij.size(), Aij.size());
A.setZero();
A.setFromTriplets(triplets.begin(), triplets.end());
}
/* -------------------------------------------------------------------------- */
void SparseSolverEigen::solve() {
const auto & Aij =
aka::as_type<SparseMatrixAIJ>(dof_manager.getMatrix(matrix_id));
auto & b_ = aka::as_type<SolverVectorDefault>(this->dof_manager.getResidual())
.getGlobalVector();
Array<Real> x_(b_.size() * b_.getNbComponent());
VectorProxy<Real> x(x_.data(), x_.size() * x_.getNbComponent());
VectorProxy<Real> b(b_.data(), b_.size() * b_.getNbComponent());
setA();
if (this->last_profile_release != Aij.getProfileRelease()) {
solver.analyzePattern(A);
this->last_profile_release = Aij.getProfileRelease();
}
if (this->last_value_release != Aij.getValueRelease()) {
solver.factorize(A);
this->last_value_release = Aij.getValueRelease();
}
x = solver.solve(b);
aka::as_type<SolverVectorDefault>(this->dof_manager.getSolution())
.setGlobalVector(x_);
this->dof_manager.splitSolutionPerDOFs();
}
} // namespace akantu
Event Timeline
Log In to Comment