main_2.cpp

#include <iostream>

#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
#include <Eigen/SparseLU>

#include <amgcl/backend/eigen.hpp>
#include <amgcl/coarsening/ruge_stuben.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/damped_jacobi.hpp>
#include <amgcl/profiler.hpp>
#include <amgcl/value_type/static_matrix.hpp>
#include <amgcl/adapter/block_matrix.hpp>

#include <vector>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <ctime>
#include "omp.h"

using namespace Eigen;
Eigen::SparseMatrix<double, RowMajor> createSparseMatrixFromCSR(int col_size, const std::vector<int> &row_offset, const std::vector<int> &col_index, const std::vector<double> &value);
Eigen::SparseMatrix<double, RowMajor> readSparseMatrixFromFile(const std::string &filename);
Eigen::VectorXd readVectorFromFile(const std::string &filename);
void usingAMGCL();
void usingEigenCG();
int main(int argc, char *argv[])
{
    Eigen::initParallel();
    usingAMGCL();
    usingEigenCG();
    return 0;
}
void usingAMGCL()
{
    amgcl::profiler<> prof;

    // Read sparse matrix from MatrixMarket format.
    // In general this should come pre-assembled.
    const std::string A_file = "matrix.txt";
    SparseMatrix<double, RowMajor> A = readSparseMatrixFromFile(A_file);

    const std::string f_file = "prd_15.txt";
    Eigen::VectorXd f = readVectorFromFile(f_file);

    // Zero initial approximation:
    Eigen::VectorXd x = Eigen::VectorXd::Constant(A.rows(), 0.0);

    // Setup the solver:
    typedef amgcl::static_matrix<double, 3, 3> mat_type;
    typedef amgcl::static_matrix<double, 3, 1> rhs_type;
    typedef amgcl::backend::builtin<mat_type> Backend;

    typedef amgcl::make_solver<
        amgcl::amg<
            Backend,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::damped_jacobi>,
        amgcl::solver::bicgstab<Backend>>
        Solver;

    Solver::params prm;
    prm.solver.maxiter = 500;
    prm.solver.tol = 10e-9;
    prm.precond.coarsening.estimate_spectral_radius=true;
    prm.precond.coarsening.power_iters=5;

    prof.tic("setup");
    Solver solve(amgcl::adapter::block_matrix<mat_type>(A), prm);
    prof.toc("setup");
    std::cout << solve << std::endl;

    // Solve the system for the given RHS:
    int iters;
    double error;
    prof.tic("solve");
    auto f_ptr = reinterpret_cast<rhs_type*>(f.data());
    auto x_ptr = reinterpret_cast<rhs_type*>(x.data());
    std::tie(iters, error) = solve(
            amgcl::make_iterator_range(f_ptr, f_ptr + A.rows() / 3),
            amgcl::make_iterator_range(x_ptr, x_ptr + A.rows() / 3));
    prof.toc("solve");

    std::cout << iters << " " << error << std::endl
              << prof << std::endl;
}
void usingEigenCG()
{
    double start, end;

    const std::string A_file = "matrix.txt";
    SparseMatrix<double, RowMajor> A = readSparseMatrixFromFile(A_file);

    const std::string f_file = "prd_15.txt";
    Eigen::VectorXd f = readVectorFromFile(f_file);

    Eigen::ConjugateGradient<Eigen::SparseMatrix<double, RowMajor>, Eigen::Lower | Eigen::Upper> cg_solver;
    cg_solver.setTolerance(1E-15);
    cg_solver.compute(A);

    start = omp_get_wtime();
    Eigen::VectorXd x = cg_solver.solve(f);
    end = omp_get_wtime();

    std::cout << "iteration is " << cg_solver.iterations() << std::endl;
    std::cout << "error is " << cg_solver.error() << std::endl;
    std::cout << "time in solve is          " << (double)(end - start) << " ms.\n";
}

Eigen::SparseMatrix<double, RowMajor> readSparseMatrixFromFile(const std::string &filename)
{
    std::ifstream fin;
    fin.open(filename.c_str());
    if ((bool)fin == false)
    {
        std::cout << "can not open file" << filename << std::endl;
    }
    int row_size;
    int col_size;
    int nonzeros;
    fin >> row_size;
    fin >> col_size;
    fin >> nonzeros;
    std::vector<int> row_offset_array(row_size + 1);
    for (int i = 0; i < row_offset_array.size(); i++)
    {
        fin >> row_offset_array[i];
    }
    std::vector<int> col_index_array(nonzeros);
    for (int i = 0; i < col_index_array.size(); i++)
    {
        fin >> col_index_array[i];
    }
    std::vector<double> value_array(nonzeros);
    for (int i = 0; i < value_array.size(); i++)
    {
        fin >> value_array[i];
    }
    fin.close();
    return createSparseMatrixFromCSR(col_size, row_offset_array, col_index_array, value_array);
}

Eigen::SparseMatrix<double, RowMajor> createSparseMatrixFromCSR(int col_size, const std::vector<int> &row_offset, const std::vector<int> &col_index, const std::vector<double> &value_array)
{
    const int row_size = row_offset.size() - 1;
    SparseMatrix<double, RowMajor> sparse_matrix(row_size, col_size);
    const int nonzeros = col_index.size();

    std::vector<Triplet<double>> triplet(nonzeros);
    for (int row = 0; row < row_size; row++)
    {
        const int start = row_offset[row];
        const int end = row_offset[row + 1];
        for (int offset_index = start; offset_index < end; offset_index++)
        {
            const int col = col_index[offset_index];
            const double value = value_array[offset_index];
            triplet[offset_index] = Triplet<double>(row, col, value);
        }
    }
    sparse_matrix.setFromTriplets(triplet.begin(), triplet.end());
    sparse_matrix.makeCompressed();
    return sparse_matrix;
}
Eigen::VectorXd readVectorFromFile(const std::string &filename)
{
    std::ifstream fin;
    fin.open(filename.c_str());
    if ((bool)fin == false)
    {
        std::cout << "can not open file" << filename << std::endl;
    }
    int size;
    fin >> size;
    VectorXd vector(size);
    for (int i = 0; i < size; i++)
    {
        fin >> vector(i);
    }
    fin.close();
    return vector;
}