Skip to content

Instantly share code, notes, and snippets.

@ddemidov
Created August 14, 2022 15:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ddemidov/74655be8a8a57aae861396c5a1a7ba40 to your computer and use it in GitHub Desktop.
Save ddemidov/74655be8a8a57aae861396c5a1a7ba40 to your computer and use it in GitHub Desktop.
amgcl without boost
#include <vector>
#include <tuple>
#define AMGCL_NO_BOOST
#include <amgcl/make_solver.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/solver/cg.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/profiler.hpp>
// Generates matrix for poisson problem in a unit cube.
template <typename PointerType, typename ColumnType, typename ValueType, typename RhsType>
int sample_problem(
ptrdiff_t n,
std::vector<PointerType> &ptr,
std::vector<ColumnType> &col,
std::vector<ValueType> &val,
std::vector<RhsType> &rhs,
double anisotropy = 1.0
)
{
ptrdiff_t n3 = n * n * n;
ptr.clear();
col.clear();
val.clear();
rhs.clear();
ptr.reserve(n3 + 1);
col.reserve(n3 * 7);
val.reserve(n3 * 7);
rhs.reserve(n3);
ValueType one = amgcl::math::identity<ValueType>();
double hx = 1;
double hy = hx * anisotropy;
double hz = hy * anisotropy;
ptr.push_back(0);
for(ptrdiff_t k = 0, idx = 0; k < n; ++k) {
for(ptrdiff_t j = 0; j < n; ++j) {
for (ptrdiff_t i = 0; i < n; ++i, ++idx) {
if (k > 0) {
col.push_back(idx - n * n);
val.push_back(-1.0/(hz * hz) * one);
}
if (j > 0) {
col.push_back(idx - n);
val.push_back(-1.0/(hy * hy) * one);
}
if (i > 0) {
col.push_back(idx - 1);
val.push_back(-1.0/(hx * hx) * one);
}
col.push_back(idx);
val.push_back((2 / (hx * hx) + 2 / (hy * hy) + 2 / (hz * hz)) * one);
if (i + 1 < n) {
col.push_back(idx + 1);
val.push_back(-1.0/(hx * hx) * one);
}
if (j + 1 < n) {
col.push_back(idx + n);
val.push_back(-1.0/(hy * hy) * one);
}
if (k + 1 < n) {
col.push_back(idx + n * n);
val.push_back(-1.0/(hz * hz) * one);
}
rhs.push_back( amgcl::math::constant<RhsType>(1.0) );
ptr.push_back( static_cast<PointerType>(col.size()) );
}
}
}
return n3;
}
//---------------------------------------------------------------------------
template <typename PointerType, typename ColumnType, typename ValueType, typename RhsType>
void solve(
// Matrix in CRS format:
const std::vector<PointerType> &ptr,
const std::vector<ColumnType> &col,
const std::vector<ValueType> &val,
// RHS:
const std::vector<RhsType> &rhs,
// Initial solution on input, solution on output:
std::vector<RhsType> &x
)
{
typedef amgcl::backend::builtin<ValueType, ColumnType, PointerType> Backend;
typedef amgcl::make_solver<
amgcl::amg<
Backend,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::spai0
>,
amgcl::solver::cg<Backend>
> Solver;
PointerType n = ptr.size() - 1;
Solver solve(std::tie(n, ptr, col, val));
std::cout << solve << std::endl;
// solve the problem:
int iters;
double error;
std::tie(iters, error) = solve(rhs, x);
std::cout << "Iterations: " << iters << std::endl;
std::cout << "Error : " << error << std::endl;
}
//---------------------------------------------------------------------------
int main(int argc, char *argv[]) {
int n = argc > 1 ? std::stoi(argv[1]) : 32;
std::vector<ptrdiff_t> ptr;
std::vector<int> col;
std::vector<double> val, rhs;
int n3 = sample_problem(n, ptr, col, val, rhs);
std::vector<double> x(n3, 0.0);
solve(ptr, col, val, rhs, x);
}
cmake_minimum_required(VERSION 3.11)
project(hello)
find_package(amgcl)
add_executable(hello amgcl.cpp)
target_link_libraries(hello amgcl::amgcl)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment