/*
	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));
}