int cg_solve(const gVar& A, gVar& x, const gVar& b, int iterationCount)
{
	double rho, newrho, alpha, beta;
	gVar p, q, r;

	r = b - A * x;

	p = r;
	rho = double(dot(r,r));

	double eps = 1e-12 * x.rowCount();

	int i;
	for (i = 0; i < iterationCount && rho > eps; i++)
	{
		q = A * p;
		alpha = rho / double(dot(p,q));

		x += alpha * p;
		r -= alpha * q;

		newrho = double(dot(r,r));
		beta = newrho / rho;

		p = r + beta * p;

		rho = newrho;
	}

	return i;
}

int cg_solve_precond_diag(const gVar& A, const gVar& M, gVar& x, const gVar& b, int iterationCount)
{
	const int correctFreq = 50;

	double rho, newrho, alpha, beta;
	gVar p, q, r, s;

	r = b - A * x;
	p = div(r, M); 

	rho = float(dot(r,p));

	double eps = 1e-12 * x.rowCount();

	int i;
	for (i = 0; i < iterationCount && rho > eps; i++)
	{
		q = A * p;
		alpha = rho / double(dot(p,q));

		x += alpha * p;
		r -= alpha * q;

		s = div(r,M);

		newrho = double(dot(r,s));
		beta = newrho / rho;

		p = s + beta * p;

		rho = newrho;
	}

	return i;
}