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

#ifdef CA_UWATERLOO_ALUMNI_DWHARDER_SPARSE_SYSTEMS

template <int M, int N>
class Matrix {
	private:
		int capacity;
		double diagonal[( M < N ) ? M : N];
		int row_index[M + 1];

		int *column_index;
		double *off_diagonal;

		static int const minMN;

		static int count_nonzero_entries( Matrix const &, Matrix const &, bool );
		void increase_capacity( int );

	public:
		Matrix( double = 0.0, int = 0 );
		Matrix( Matrix const & );
		~Matrix();

		Matrix &operator=( Matrix const & );

		double operator()( int, int ) const;

		// Special Constructors

		// ...

		// Matrix Functions
		double norm() const;
		double norm( int ) const;
		double trace() const;
		Matrix operator- () const;
		bool validate() const;

		// Matrix Mechanics

		void set( int, int, double );
		void set_symmetric( int, int, double );

		void clear( int, int );
		void clear();

		void resize( int = 0 );

		// Matrix-Matrix Operations

		Matrix &operator+=( Matrix const & );
		Matrix &operator-=( Matrix const & );

		// Matrix-Vector Operations

		template <Orientation D>
		Matrix &operator+=( Vector<M < N ? M : N, D> const & );

		template <Orientation D>
		Matrix &operator-=( Vector<M < N ? M : N, D> const & );

		// Matrix-Scalar Operations

		Matrix &operator+=( double );
		Matrix &operator-=( double );
		Matrix &operator*=( double );
		Matrix &operator/=( double );

		// Matrix Queries

		bool is_lower_triangular() const;
		bool is_upper_triangular() const;

		bool is_strict_lower_triangular() const;
		bool is_strict_upper_triangular() const;

		bool is_diagonal() const;
		bool is_symmetric() const;
		bool is_antisymmetric() const;

		// Matrix Streaming/Printing
		void details() const;
		void matlab( std::string const & ) const;
		void matlab_dense( std::string const & ) const;

	template <int Mt, int Nt>
	friend Matrix<Mt, Nt> operator+( Matrix<Mt, Nt> const &, Matrix<Mt, Nt> const & );

	// Matrix Functions

	template <int Mt, int Nt>
	friend Matrix<Nt, Mt> transpose( Matrix<Mt, Nt> const & );

	// Matrix-Vector Operators

	template <int Mt, int Nt>
	friend Vector<Nt, ROW> operator*( Vector<Mt, ROW> const &, Matrix<Mt, Nt> const & );

	template <int Mt, int Nt>
	friend Vector<Mt, COLUMN> operator*( Matrix<Mt, Nt> const &, Vector<Nt, COLUMN> const & );

	template <int Nt>
	friend Vector<Nt, COLUMN> &operator*=( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> & );

	template <int Nt>
	friend Vector<Nt, ROW> operator*=( Vector<Nt, ROW> &, Matrix<Nt, Nt> const & );

	// Matrix-Scalar Functions

	template <int Mt, int Nt>
	friend Matrix<Mt, Nt> operator-( double, Matrix<Mt, Nt> const & );

	// Matrix Solvers

	template <int Nt>
	friend Vector<Nt, COLUMN> jacobi( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> const &, int, double );

	template <int Nt>
	friend Vector<Nt, COLUMN> gauss_seidel( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> const &, int, double );

	template <int Nt>
	friend Vector<Nt, COLUMN> diagonal_solve( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> const & );

	template <int Nt>
	friend Vector<Nt, COLUMN> backward_substitution( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> const &, bool );

	template <int Nt>
	friend Vector<Nt, COLUMN> forward_substitution( Matrix<Nt, Nt> const &, Vector<Nt, COLUMN> const &, bool );

	template <int Mt, int Nt>
	friend Vector<Mt, COLUMN> matrix_vector_product( Matrix<Mt, Nt> const &, Vector<Nt, COLUMN> const &, Product );

	template <int Mt, int Nt>
	friend Vector<Nt, ROW> vector_matrix_product( Vector<Mt, ROW> const &, Matrix<Mt, Nt> const &, Product );

	// Matrix Streaming/Printing
	template <int Mt, int Nt>
	friend std::ostream &operator<<( std::ostream &, Matrix<Mt, Nt> const & );
};

template <int M, int N>
int const Matrix<M, N>::minMN = ( M < N ) ? M : N;

#include "matrix_mechanics.cpp"
#include "matrix_functions.cpp"
#include "matrix_private_functions.cpp"
#include "matrix_queries.cpp"
#include "matrix_solvers.cpp"
#include "matrix_scalar_operations.cpp"
#include "matrix_vector_operations.cpp"
#include "matrix_matrix_operations.cpp"
#include "matrix_streaming.cpp"
#include "projection_methods.cpp"

#endif