amgcl with sparse matrix conditioner and self-defined matrix-vector product

amgcl_w_my_class.cpp

// amgcl_w_my_class.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include "my_solver.cpp"


int main()
{
    std::cout << "Hello World!\n";
	int N = 100;
	my_solver* test = new my_solver(N, 20.);
	std::complex<double> * b_in = (std::complex<double>*) malloc(N * sizeof(std::complex<double>));
	std::complex<double> * x_in = (std::complex<double>*) malloc(N * sizeof(std::complex<double>));

	test->gmres_amgcl_complex<double>(b_in, x_in);

	return 0;

}

my_solver.cpp

#include "my_solver.h"

namespace amgcl {
	namespace backend {

		template <class Alpha, class Vector1, class Beta, class Vector2>
		struct spmv_impl<Alpha, /* this is my struct: */my_solver, Vector1, Beta, Vector2> {
			static void apply(Alpha alpha, const my_solver &A, const Vector1 &x, Beta beta, Vector2 &y)
			{
				double omega = A.get_omega();
				typedef double T;
				std::complex<T>* v = &x[0];
				std::complex<T>* lhs = (std::complex<T>*) malloc(y.size() * sizeof(std::complex<T>));
				// 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>(lhs, v, omega);
				for (int i = 0; i < y.size(); ++i) {
					//y[i] = beta * y[i] + alpha * lhs[i];
					y[i] = 0 * y[i] + 1. * lhs[i];
				}
				delete[]lhs;
			}
		};
		template <class my_solver, class Vector1, class Vector2, class Vector3>
		struct residual_impl<
			my_solver, Vector1, Vector2, Vector3
		>
		{
			static void apply(
				Vector1 const &rhs,
				my_solver  const &A,
				Vector2 const &x,
				Vector3       &res
			)
			{
				double omega = A.get_oemga();
				typedef double T;
				std::complex<T> *x_in = &x[0];
				std::complex<T>* lhs_complex = (std::complex<T>*)malloc(rhs.size() * sizeof(std::complex<T>));
				A.lhs_gmres_complex<T>(lhs_complex, x, omega);
				for (int p = 0; p < x.size(); p++)
				{
					res[p] = rhs[p] - lhs_complex[p];
				}
				delete[]lhs_complex;
			}
		};
	} // namespace backend
} //namespace amgcl


my_solver::my_solver(int size_in, double omega) : size_in(size_in), omega(omega) {}


double my_solver::get_omega() { return this->omega;}

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

	return 0;
}

template<typename T>
int my_solver::gmres_amgcl_complex(
	std::complex<T> * b_in,
	std::complex<T> *x_in
	)
{
	typedef amgcl::backend::builtin<std::complex<T>> 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;

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

	ia[0] = 0;
	for (int i = 0; i < this->size_in; i++)
	{
		ja[i] = i;
		ia[i + 1] = i;
	}
	// constructor with jac
	Solver solve(std::tie(size_in, ia, ja, a));


	int    iters;
	double error;
	std::tie(iters, error) = solve(*this, b_in, x_in);
	return 0;
}

my_solver.h

#pragma once

#include <complex>
#include <cstdlib>

#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>
class my_solver
{
	typedef std::complex<double> value_type;
	private:
		double omega;
		int size_in;
    public:
		double get_omega();


		my_solver(int size_in, double omega);
		template <typename T>
		int lhs_gmres_complex(std::complex<T> *w, std::complex<T> *v, double omega);
		
		template<typename T>
		int gmres_amgcl_complex(
			std::complex<T> * b_in,
			std::complex<T> *x_in
			);

};