//
// Big Vector and Sparse Matrix Classes
//

#include <float.h>

#include "vec3n.h"

float conjgrad_lasterror;
float conjgrad_epsilon = 0.1f;
int conjgrad_loopcount;
int conjgrad_looplimit = 100;

/*EXPORTVAR(conjgrad_lasterror);
EXPORTVAR(conjgrad_epsilon  );
EXPORTVAR(conjgrad_loopcount);
EXPORTVAR(conjgrad_looplimit);
*/

int ConjGradient(float3N &X, float3Nx3N &A, float3N &B)
{
	// Solves for unknown X in equation AX=B
	conjgrad_loopcount = 0;
	int n = B.count;
	float3N q(n), d(n), tmp(n), r(n);
	r = B - Mul(tmp, A, X);  // just use B if X known to be zero
	d = r;
	float s = dot(r, r);
	float starget = s * squared(conjgrad_epsilon);
	while (s > starget && conjgrad_loopcount++ < conjgrad_looplimit)
	{
		Mul(q, A, d);  // q = A*d;
		float a = s / dot(d, q);
		X = X + d * a;
		if (conjgrad_loopcount % 50 == 0)
		{
			r = B - Mul(tmp, A, X);
		}
		else
		{
			r = r - q * a;
		}
		float s_prev = s;
		s = dot(r, r);
		d = r + d * (s / s_prev);
	}
	conjgrad_lasterror = s;
	return conjgrad_loopcount < conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
}

int ConjGradientMod(float3N &X, float3Nx3N &A, float3N &B, int3 hack)
{
	// obsolete!!!
	// Solves for unknown X in equation AX=B
	conjgrad_loopcount = 0;
	int n = B.count;
	float3N q(n), d(n), tmp(n), r(n);
	r = B - Mul(tmp, A, X);  // just use B if X known to be zero
	r[hack[0]] = r[hack[1]] = r[hack[2]] = float3(0, 0, 0);
	d = r;
	float s = dot(r, r);
	float starget = s * squared(conjgrad_epsilon);
	while (s > starget && conjgrad_loopcount++ < conjgrad_looplimit)
	{
		Mul(q, A, d);  // q = A*d;
		q[hack[0]] = q[hack[1]] = q[hack[2]] = float3(0, 0, 0);
		float a = s / dot(d, q);
		X = X + d * a;
		if (conjgrad_loopcount % 50 == 0)
		{
			r = B - Mul(tmp, A, X);
			r[hack[0]] = r[hack[1]] = r[hack[2]] = float3(0, 0, 0);
		}
		else
		{
			r = r - q * a;
		}
		float s_prev = s;
		s = dot(r, r);
		d = r + d * (s / s_prev);
		d[hack[0]] = d[hack[1]] = d[hack[2]] = float3(0, 0, 0);
	}
	conjgrad_lasterror = s;
	return conjgrad_loopcount < conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
}

static inline void filter(float3N &V, const float3Nx3N &S)
{
	for (int i = 0; i < S.blocks.count; i++)
	{
		V[S.blocks[i].r] = V[S.blocks[i].r] * S.blocks[i].m;
	}
}

int ConjGradientFiltered(float3N &X, const float3Nx3N &A, const float3N &B, const float3Nx3N &S)
{
	// Solves for unknown X in equation AX=B
	conjgrad_loopcount = 0;
	int n = B.count;
	float3N q(n), d(n), tmp(n), r(n);
	r = B - Mul(tmp, A, X);  // just use B if X known to be zero
	filter(r, S);
	d = r;
	float s = dot(r, r);
	float starget = s * squared(conjgrad_epsilon);
	while (s > starget && conjgrad_loopcount++ < conjgrad_looplimit)
	{
		Mul(q, A, d);  // q = A*d;
		filter(q, S);
		float a = s / dot(d, q);
		X = X + d * a;
		if (conjgrad_loopcount % 50 == 0)
		{
			r = B - Mul(tmp, A, X);
			filter(r, S);
		}
		else
		{
			r = r - q * a;
		}
		float s_prev = s;
		s = dot(r, r);
		d = r + d * (s / s_prev);
		filter(d, S);
	}
	conjgrad_lasterror = s;
	return conjgrad_loopcount < conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable
}

// test big vector math library:
static void testfloat3N()
{
	float3N a(2), b(2), c(2);
	a[0] = float3(1, 2, 3);
	a[1] = float3(4, 5, 6);
	b[0] = float3(10, 20, 30);
	b[1] = float3(40, 50, 60);
	//      c =  a+b+b * 10.0f;
	//      float d = dot(a+b,-b);
	int k;
	k = 0;
}
class dotest
{
public:
	dotest() { testfloat3N(); }
} do_test_at_program_startup;