// Author:  Douglas Wilhelm Harder
// Copyright (c) 2009 by Douglas Wilhelm Harder.  All rights reserved.

#include "Sparse.h"
#include <iostream>
using namespace std;

int main() {
	int const M = 3;
	int const N = 3;
	
	Matrix<M, N> A(3);
	Matrix<M, N> B(4);
	Matrix<M, N> C(5);
	Matrix<M, N> D;
	
	A.set(0,2, 0.5);
	A.set(2,0,-0.7);
	
	B.set(0,2,0.3);
	B.set(2,0,0.2);
	B.set(2,1,0.4);
	B.set(1,2,0.1);
	B.set(0,1,0.5);
	B.set(1,0,0.6);
	
	cout << "A = " << A << endl;
	A.matlab( "A" );
	cout << "B = " << B << endl;
	B.matlab( "B" );
	cout << "C = " << C << endl;
	C.matlab( "C" );

	Vector<N> v;

	v(0) =  0.3;
	v(1) = -0.5;
	v(2) =  0.8;

	D += v;
	cout << "D = " << D << endl;
	D.matlab( "D" );

	A.set(0,2, 1.5);
	A.set(2,0,-1.5);

	cout << "A = " << A << endl;
	A.matlab( "A" );

	cout.precision( 10 );

	cout << "norm( A, 1 ) = " << norm( A, 1 ) << endl;
	cout << "norm( A, 2 ) = " << norm( A, 2 ) << endl;
	cout << "norm( A, inf ) = " << norm( A ) << endl;

	cout << "norm( B, 1 ) = " << norm( B, 1 ) << endl;
	cout << "norm( B, 2 ) = " << norm( B, 2 ) << endl;
	cout << "norm( B, inf ) = " << norm( B ) << endl;

	cout << "norm( C, 1 ) = " << norm( C, 1 ) << endl;
	cout << "norm( C, 2 ) = " << norm( C, 2 ) << endl;
	cout << "norm( C, inf ) = " << norm( C ) << endl;

	cout << "norm( D, 1 ) = " << norm( D, 1 ) << endl;
	cout << "norm( D, 2 ) = " << norm( D, 2 ) << endl;
	cout << "norm( D, inf ) = " << norm( D ) << endl;

	Matrix<M, N> E( A );

	cout << E << endl;

	Matrix<M, N> F = -A;

	cout << F << endl;

	cout << "A = " << A << endl;
	cout << "2*A = " << 2*A << endl;
	cout << "B = " << B << endl;
	cout << "2*B = " << 2*B << endl;
	cout << "C = " << C << endl;
	cout << "2*C = " << 2*C << endl;
	cout << "D = " << D << endl;
	cout << "2*D = " << 2*D << endl;

	cout << "trace( A ) = " << trace( A ) << endl;
	cout << "trace( B ) = " << trace( B ) << endl;
	cout << "trace( C ) = " << trace( C ) << endl;
	cout << "trace( D ) = " << trace( D ) << endl;

	cout << ">> A = " << A << endl;
	cout << ">> F = " << F << endl;

	cout << "A + A = " << endl << A + A << endl;
	cout << "A + B = " << endl << A + B << endl;
	cout << "A + C = " << endl << A + C << endl;
	cout << "A + D = " << endl << A + D << endl;
	cout << "A + E = " << endl << A + E << endl;
	cout << "A + F = " << endl << A + F << endl;
	cout << "B + C = " << endl << B + C << endl;
	cout << "B + D = " << endl << B + D << endl;
	cout << "B + E = " << endl << B + E << endl;
	cout << "B + F = " << endl << B + F << endl;
	cout << "C + D = " << endl << C + D << endl;
	cout << "C + E = " << endl << C + E << endl;
	cout << "C + F = " << endl << C + F << endl;
	cout << "D + E = " << endl << D + E << endl;
	cout << "D + F = " << endl << D + F << endl;
	cout << "E + F = " << endl << E + F << endl;

	Matrix<M, N> G = A + F;
	G.matlab( "G" );

	A.details();
	E.details();

	cout << ">> A = " << A << endl;
	cout << ">> B = " << B << endl;
	cout << ">> C = " << C << endl;
	cout << ">> D = " << D << endl;
	cout << ">> E = " << E << endl;
	cout << ">> F = " << F << endl;

	cout << "A - A = " << endl << A - A << endl;
	cout << "A - B = " << endl << A - B << endl;
	cout << "A - C = " << endl << A - C << endl;
	cout << "A - D = " << endl << A - D << endl;
	cout << "A - E = " << endl << A - E << endl;
	cout << "A - F = " << endl << A - F << endl;
	cout << "B - C = " << endl << B - C << endl;
	cout << "B - D = " << endl << B - D << endl;
	cout << "B - E = " << endl << B - E << endl;
	cout << "B - F = " << endl << B - F << endl;
	cout << "C - D = " << endl << C - D << endl;
	cout << "C - E = " << endl << C - E << endl;
	cout << "C - F = " << endl << C - F << endl;
	cout << "D - E = " << endl << D - E << endl;
	cout << "D - F = " << endl << D - F << endl;
	cout << "E - F = " << endl << E - F << endl;
	cout << "F - F = " << endl << F - F << endl;

	cout << ">> F = " << F << endl;

	Vector<3> u;

	u = jacobi( A, v, 1000, 1e-40 );
	cout << "norm( Au - v ) = " << norm( A*u - v ) << endl;
	u = jacobi( B, v, 1000, 1e-40 );
	cout << "norm( Bu - v ) = " << norm( B*u - v ) << endl;
	u = jacobi( C, v, 1000, 1e-40 );
	cout << "norm( Cu - v ) = " << norm( C*u - v ) << endl;
	u = jacobi( D, v, 1000, 1e-40 );
	cout << "norm( Du - v ) = " << norm( D*u - v ) << endl;
	u = jacobi( E, v, 1000, 1e-40 );
	cout << "norm( Eu - v ) = " << norm( E*u - v ) << endl;
	u = jacobi( F, v, 1000, 1e-40 );
	cout << "norm( Fu - v ) = " << norm( F*u - v ) << endl;

	u = gauss_seidel( A, v, 1000, 1e-40 );
	cout << "norm( Au - v ) = " << norm( A*u - v ) << endl;
	u = gauss_seidel( B, v, 1000, 1e-40 );
	cout << "norm( Bu - v ) = " << norm( B*u - v ) << endl;
	u = gauss_seidel( C, v, 1000, 1e-40 );
	cout << "norm( Cu - v ) = " << norm( C*u - v ) << endl;
	u = gauss_seidel( D, v, 1000, 1e-40 );
	cout << "norm( Du - v ) = " << norm( D*u - v ) << endl;
	u = gauss_seidel( E, v, 1000, 1e-40 );
	cout << "norm( Eu - v ) = " << norm( E*u - v ) << endl;
	u = gauss_seidel( F, v, 1000, 1e-40 );
	cout << "norm( Fu - v ) = " << norm( F*u - v ) << endl;

	return 0;
}