amgcl_w_my_class.cpp

#include "my_class.cpp"



int main()
{
 
	//std::cout << iters << " " << error << std::endl;

	my_solver<double> *system = new my_solver<double>(50, 7);
	std::complex<double> *x = new std::complex<double>(50);
	std::complex<double> *rhs = new std::complex<double>(50);
	for (int p = 0; p < 50; p++)
	{
		x[p] = p;
		rhs[p] = 0.1 * p + std::pow(p, 2.0);
	}


	system->gmres_amgcl(x, rhs);
	return 0;

}

my_class.cpp

#include "my_class.h"
namespace amgcl {
	namespace backend {

		template <class Alpha, class T, class Vector1, class Beta, class Vector2>
		struct spmv_impl<Alpha, my_solver<T>, Vector1, Beta, Vector2> {
			static void apply(Alpha alpha, const my_solver<T> &A, const Vector1 &x, Beta beta, Vector2 &y)
			{
				double omega = A.omega;
				std::vector<std::complex<T>> t(y.size());
				// This should implement operation y = beta * y + alpha * A * x
				// you can probably adjust implementation of your `lhs_gmres_solver_complex()` function,
				// but right now you would need a temporary vector to do this:
				//std::vector<std::complex<T>> t(y.size());
				A.lhs_gmres_complex(&t[0], &x[0], omega);
				for (int i = 0; i < y.size(); ++i) {
					y[i] = beta * y[i] + alpha * t[i];
					//y[i] = 0 * y[i] + 1. * lhs[i];
				}
			}
		};

		template <class T, class Vector1, class Vector2, class Vector3>
		struct residual_impl<
			my_solver<T>, Vector1, Vector2, Vector3
		>
		{
			static void apply(
				const Vector1 &rhs,
				const my_solver<T> &A,
				const Vector2 &x,
				Vector3       &res
			)
			{
				double omega = A.omega;
				A.lhs_gmres_complex(&res[0], &x[0], omega);
				for (int p = 0; p < x.size(); p++)
				{
					res[p] = rhs[p] - res[p];
				}
			}
		};
	} // namespace backend
} //namespace amgcl
template <typename T>
double my_solver<T>::get_omega() { return omega; }

template <typename T>
int my_solver<T>::lhs_gmres_complex(std::complex<T>* w, const std::complex<T>* v, double omega) const
{
	for (int i = 0; i < size_in; i++)
	{
		w[i] = 2. * v[i] - omega * i / 2.;
	}

	return 0;
}


template <typename T>
int my_solver<T>::gmres_amgcl(std::complex<double> *x, std::complex<double> *rhs)
{
	std::cout << "Hello World!\n";

	std::vector<std::complex<double>> b_in(rhs, rhs + this->size_in);
	std::vector<std::complex<double>> x_in(this->size_in, std::complex<double>(0, 0));
	std::vector<int> ja(this->size_in);
	std::vector<int> ia(this->size_in + 1);
	std::vector<std::complex<double>> a(this->size_in, std::complex<double>(1, 0));

	ia[0] = 0;
	for (int i = 0; i < this->size_in; ++i)
	{
		ja[i] = i;
		ia[i + 1] = i + 1;
	}

	typedef amgcl::backend::builtin< std::complex<double> > Backend;
	typedef amgcl::make_solver<
		// Use AMG as preconditioner:
		amgcl::amg<
		Backend,
		amgcl::coarsening::smoothed_aggregation,
		amgcl::relaxation::spai0
		>,
		// And BiCGStab as iterative solver:
		amgcl::solver::fgmres<Backend>
	> Solver;


	Solver solve(std::tie(this->size_in, ia, ja, a));

	int    iters;
	double error;
	std::tie(iters, error) = solve(this, b_in, x_in);
	x = &x_in[0];
  return 0;
}

my_class.h

#pragma once

#include <iostream>
#include <complex>
#include <cstdlib>
#include <vector>

#include <amgcl/util.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/solver/fgmres.hpp>
#include <amgcl/value_type/complex.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/adapter/reorder.hpp>


template <typename T>
class my_solver
{

public:
	typedef std::complex<T> value_type;
	double omega;
	int size_in;

	my_solver(int size_in, double omega) : size_in(size_in), omega(omega) {}
	
	double get_omega();
	int lhs_gmres_complex(std::complex<T>* w, const std::complex<T>* v, double omega) const;


	int gmres_amgcl(std::complex<double> *x, std::complex<double> *rhs);
};