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

#include <cmath>
#include <cassert>
#include "Heap.h"

/***************************************************
 * The infinity norm of the matrix:
 *  - The maximum absolute row sum of the matrix
 *
 * Run time:   O(M + n)
 * Memory:     O(1)
 ***************************************************/

template <int M, int N>
double Matrix<M, N>::norm() const {
	double value = 0;

	for ( int i = 0; i < M; ++i ) {
		double row_sum;

		if ( i < N ) {
			row_sum = std::fabs( diagonal[i] );
		} else {
			row_sum = 0.0;
		}

		for ( int j = row_index[i]; j < row_index[i + 1]; ++j ) {
			row_sum += std::fabs( off_diagonal[j] );
		}

		value = std::max( value, row_sum );
	}

	return value;
}

template <int M, int N>
double norm( Matrix<M, N> const &A ) {
	return A.norm();
}

/***************************************************
 * The 1- and 2-norms of the matrix, that is
 *  - The maximum absolute column sum of the matrix
 *  - The largest singular value of the matrix
 * respectively.
 *
 * 1-Norm:
 * Run time:   O(M + N + n)
 * Memory:     O(N)
 *
 * 2-Norm:
 * Run time:   O(K(M + N + n))
 * Memory:     O(N)
 ***************************************************/

template <int M, int N>
double Matrix<M, N>::norm( int n ) const {
	if ( n == 1 ) {
		double column_sum[N];

		for ( int i = 0; i < N; ++i ) {
			if ( i < M ) {
				column_sum[i] = std::fabs( diagonal[i] );
			} else {
				column_sum[i] = 0.0;
			}
		}

		for ( int i = 0; i < M; ++i ) {
			for ( int j = row_index[i]; j < row_index[i + 1]; ++j ) {
				column_sum[column_index[j]] += std::fabs( off_diagonal[j] );
			}
		}

		double value = 0;

		for ( int i = 0; i < N; ++i ) {
			value = std::max( value, column_sum[i] );
		}

		return value;
	} else if ( n == 2 ) {
		// This can be more efficient...
		Vector<N> v0;
		double d0;

		Vector<N> v = Vector<N>::random();
		double d = 1e100;
		int i;

		do {
			++i;
			v /= v.norm(2);
			v0 = v;
			d0 = d;
			v = transpose(transpose(*this * v) * *this);
			d = (transpose(v)*v) / (transpose(v0)*v0);
		} while ( std::fabs( d0 - d ) > 1e-20 );

		return std::sqrt( d );
	}

	throw invalid_norm( n );
}

template <int M, int N>
double norm( Matrix<M, N> const &A, int n ) {
	return A.norm( n );
}

/***************************************************
 * The trace of a matrix.
 *  - The sum of the diagonal entries of the
 *    matrix.
 *
 * The sum of the diagonal entries of the matrix.
 * Run time:   O(min(M, N))
 * Memory:     O(1)
 ***************************************************/

template <int M, int N>
double Matrix<M, N>::trace() const {
	double sum = 0;

	for ( int i = 0; i < minMN; ++i ) {
		sum += diagonal[i];
	}

	return sum;
}

template <int M, int N>
double trace( Matrix<M, N> const &A ) {
	return A.trace();
}

/***************************************************
 * The transpose of a matrix.
 *
 * Run time:   O(M + N + na ln(na))
 * Memory:     O(N + na)
 ***************************************************/

template <int M, int N>
Matrix<N, M> transpose( Matrix<M, N> const &A ) {
	Matrix<N, M> B( 0.0, A.capacity );

	for ( int i = 0; i < Matrix<M, N>::minMN; ++i ) {
		B.diagonal[i] = A.diagonal[i];
	}

	Heap hp( A.row_index[M] );

	for ( int i = 0; i < M; ++i ) {
		for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
			hp.push( A.column_index[j], i, A.off_diagonal[j] ); 
		}
	}

	B.row_index[0] = 0;
	int i = 0;

	for ( int j = 0; j < A.row_index[M]; ++j ) {
		int row =      hp.top_row();
		int column =   hp.top_column();
		double value = hp.top_value();

		hp.pop();

		while ( i < row ) {
			++i;
			B.row_index[i + 1] = B.row_index[i];
		}

		++(B.row_index[i + 1]);
		B.column_index[j] = column;
		B.off_diagonal[j] = value;
	}

	while ( i < N - 1 ) {
		++i;
		B.row_index[i + 1] = B.row_index[i];
	}


	return B;
}

/***************************************************
 * Calculate -A
 *
 * Run time:   O( M + n )
 * Memory:     O( M + n )
 ***************************************************/

template <int M, int N>
Matrix<M, N> Matrix<M, N>::operator- () const {
	Matrix<M, N> A;

	for ( int i = 0; i < minMN; ++i ) {
		A.diagonal[i] = -diagonal[i];
	}

	for ( int i = 0; i <= M; ++i ) {
		A.row_index[i] = row_index[i];
	}

	for ( int i = 0; i < row_index[M]; ++i ) {
		A.column_index[i] =  column_index[i];
		A.off_diagonal[i] = -off_diagonal[i];
	}

	return A;
}

/***************************************************
 * Validate the structure of the matrix.
 *  - Check that the internal structure is self-
 *    consistent:
 *     - The capacity is large enough
 *     - row_index is monotonically increasing
 *     - column_index is strictly monotonically
 *       increasing on each interval
 *         row_index[i], ..., row_index[i + 1] - 1
 *       for i = 0, ..., M - 1
 *     - diagonal entries are not stored in the
 *       off-diagonal
 *     - zeros are not stored in the off-diagonal
 *
 * Run time:   O( M + n )
 * Memory:     O( M + n )
 ***************************************************/

template <int M, int N>
bool Matrix<M, N>::validate() const {
	int is_valid = true;

	if ( capacity < row_index[M] ) {
		std::cerr << "Capacity " << capacity
			<< " too small for the number of entries:  " << row_index[M] << std::endl;

		is_valid = false;
	}

	if ( row_index[0] != 0 ) {
		std::cerr << "Invalid row_index[0] == " << row_index[0] << std::endl;

		is_valid = false;
	}

	for ( int i = 0; i < M; ++i ) {
		if ( row_index[i] > row_index[i + 1] ) {
			std::cerr << "row_index[" << i << "] == "
				<< row_index[i] << " > row_index[" << (i + 1)
				<< "] == " << row_index[i + 1] << std::endl;

			is_valid = false;
		}

		for ( int j = row_index[i] + 1; j < row_index[i + 1]; ++j ) {
			if ( column_index[j - 1] > column_index[j] ) {
				std::cerr << "Within row_index[" << i << "] == "
					<< row_index[i] << " and  row_index[" << (i + 1)
					<< "] == " << row_index[i + 1] << ":" << std::endl
					<< "column_index[" << (j - 1) << "] ==
					" << column_index[j - 1] << " > column_index[" << j
					<< "] == " << column_index[j] << std::endl;

				is_valid = false;
			}

			if ( column_index[j - 1] == column_index[j] ) {
				std::cerr << "Within row_index[" << i << "] == "
					<< row_index[i] << " and  row_index[" << (i + 1)
					<< "] == " << row_index[i + 1] << ":" << std::endl
					<< "column_index[" << (j - 1) << "] == " << column_index[j - 1]
					<< " equals column_index[" << j << "] == " << column_index[j]
					<< std::endl;

				is_valid = false;
			}

			if ( column_index[j] == i ) {
				std::cerr << "Within row_index[" << i << "] == " << row_index[i]
					<< " and  row_index[" << (i + 1) << "] == " << row_index[i + 1]
					<< ":" << std::endl << "column_index[" << j << "] equals i == "
					<< i << " (which would be on the diagonal)" << std::endl;

				is_valid = false;
			}

			if ( off_diagonal[j] == 0 ) {
				std::cerr << "Within row_index[" << i << "] == " << row_index[i]
					<< " and  row_index[" << (i + 1) << "] == " << row_index[i + 1]
					<< ":" << std::endl << "column_index[" << j << "] == "
					<< column_index[j] << " contains a zero entry at off_diagonal["
					<< j << "]" << std::endl;

				is_valid = false;
			}
		}
	}

	return is_valide;
}