/*
	David Baraff
	School of Computer Science
	Carnegie Mellon University

	Simple test --

		Make m1 be a random 10 x 10 matrix,
		and b a random 10-vector.

		Let m2 = transpose(m2)

		Let m3 = m1 * m2, which is symmetric positive definite
			(probably)

		Save m4 = m3 = m1 * m2

		Solve
			m3 * x = b	(m3 gets nuked by this)

		Compute
			prod = m4 * x

		Compare prod, and b -- should be the same...
*/


#include	"matrix.h"

int	RANDIRANGE(int a, int b)
					{ return (random()%(b-a+1))+a; }
double	RANDRANGE(double a, double b)
{
	return ((RANDIRANGE(0,10000)/10000.) * (b - a)) + a;

}

main()
{
	Matrix	*m1 = MatrixAlloc(10,10),
		*m2 = MatrixAlloc(10,10),
		*m3 = MatrixAlloc(10,10),
		*m4 = MatrixAlloc(10,10);
	Vector	*b = VectorAlloc(10),
		*x = VectorAlloc(10),
		*prod = VectorAlloc(10);

	int	i, j;

	for(i = 0; i < 10; i++)
	{
		V_ELEMENT(b,i) = RANDRANGE(-1.,1.);
		for(j = 0; j < 10; j++)
		{
			M_ELEMENT(m1,i,j) = RANDRANGE(-1.,1.);
		}
	}

	MatrixTranspose(m1,m2);
	MatrixTimesMatrix(m1, m2, m3);
	MatrixTimesMatrix(m1, m2, m4);
	ChSolve(m3, b, x);
	MatrixTimesVector(m4, x, prod);

	for(i = 0; i < 10; i++)
		printf("B = %g,  MX = %g\n", V_ELEMENT(b,i),
		       V_ELEMENT(prod,i));
}