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

#include <cmath>
#include <cassert>

/***************************************************
 * Constructor
 *
 * Run time:   O( M )
 * Memory:     O( M + cap )
 ***************************************************/

template <int M, int N>
Matrix<M, N>::Matrix( double x, int cap ):
capacity( cap ),
column_index( new int[capacity] ),
off_diagonal( new double[capacity] ) {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] = x;
	}

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

/***************************************************
 * Copy Constructor
 * Make a copy of the matrix A
 * with the same capacity.
 *
 * Run time:   O( M + na )
 * Memory:     O( M + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>::Matrix( Matrix<M, N> const &A ):
capacity( A.capacity ),
column_index( new int[capacity] ),
off_diagonal( new double[capacity] ) {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] = A.diagonal[i];
	}

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

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

/***************************************************
 * Destructor
 *  - Delete those arrays dynamically created
 *    with new[]
 *
 * Run time:   O( 1 )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
Matrix<M, N>::~Matrix() {
	delete [] column_index;
	delete [] off_diagonal;
}

/***************************************************
 * Assignment Operator
 *  - Assign the matrix A to this by making a
 *    deep copy.
 *
 * Run time:   O( M + na )
 * Memory:     O( M + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N> &Matrix<M, N>::operator= ( Matrix<M, N> const &A ) {
	// Check for an assignment of the form A = A

	if ( &A == this ) {
		return *this;
	}

	if ( capacity != A.capacity ) {
		// Delete the current matrix structure.
		delete [] column_index;
		delete [] off_diagonal;

		capacity = A.capacity;
		column_index = new int[capacity];
		off_diagonal = new double[capacity];
	}

	// Copy over the entries of the four arrays

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

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

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

	return *this;
}

/***************************************************
 * Resize the array for storing
 * the off-diagonal entries.
 *
 * Run time:   O( cap )
 * Memory:     O( cap )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::resize( int cap ) {
	if ( capacity != std::max( cap, row_index[M] ) ) {
		increase_capacity( std::max( cap, row_index[M] ) );
	}
}

/***************************************************
 * Assign A(i,j) = value
 *
 * Run time:   Diagonal:      O( 1 )
 *             Off-diagonal:  O( n/M )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::set( int i, int j, double value ) {
	assert( i >= -M && j >= -N && i < M && j < N );

	i += ( i < 0 ) ? M : 0;
	j += ( j < 0 ) ? N : 0;

	// If the entry being accessed is on the diagonal, i.e., i == j,
	// return the corresponding diagonal entry

	if ( i == j ) {
		diagonal[i] = value;
		return;
	}

	if ( value == 0 ) {
		clear( i, j );
		return;
	}

	for ( int k = row_index[i]; k < row_index[i + 1]; ++k ) {
		if ( column_index[k] == j ) {
			off_diagonal[k] = value;
			return;
		}

		if ( column_index[k] > j ) {
			// Allocate new memory

			if ( row_index[M] == capacity ) {
				increase_capacity( capacity + 10 );
			}

			for ( int ell = capacity - 1; ell > k; --ell ) {
				column_index[ell] = column_index[ell - 1];
				off_diagonal[ell] = off_diagonal[ell - 1];
			}

			column_index[k] = j;
			off_diagonal[k] = value;

			for ( int ell = i + 1; ell <= M; ++ell ) {
				++row_index[ell];
			}

			return;
		}
	}

	// At this point, we have passed all entries in row i of the matrix A

	if ( row_index[M] == capacity ) {
		increase_capacity( capacity + 10 );
	}

	for ( int ell = capacity - 1; ell > row_index[i + 1]; --ell ) {
		column_index[ell] = column_index[ell - 1];
		off_diagonal[ell] = off_diagonal[ell - 1];
	}

	column_index[row_index[i + 1]] = j;
	off_diagonal[row_index[i + 1]] = value;

	for ( int ell = i + 1; ell <= M; ++ell ) {
		++row_index[ell];
	}
}

/***************************************************
 * Assign A(i,j) = value
 *        A(j,i) = value
 *
 * Run time:   Diagonal:      O( 1 )
 *             Off-diagonal:  O( n/M )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::set_symmetric( int i, int j, double value ) {
	assert( i >= -M && j >= -N && i < M && j < N );

	i += ( i < 0 ) ? M : 0;
	j += ( j < 0 ) ? N : 0;

	if ( i == j ) {
		set( i , j, value );
	} else if ( i < j ) {
		set( i, j, value );
		set( j, i, value );
	} else {
		set( j, i, value );
		set( i, j, value );
	}
}

/***************************************************
 * Set the entry A(i, j) = 0
 *  - On the diagonal, just set it, and
 *  - Off the diagaonl, remove the appropriate
 *    entry in 'column_index' and 'off_diagonal'
 *    and update 'row_index', if necessary.
 *
 * Run time:   O( M + n )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::clear( int i, int j ) {
	assert( i >= -M && j >= -N && i < M && j < N );

	i += ( i < 0 ) ? M : 0;
	j += ( j < 0 ) ? N : 0;

	if ( i == j ) {
		diagonal[i] = 0.0;
	} else {
		for ( int k = row_index[i]; k < row_index[i + 1]; ++k ) {
			if ( column_index[k] == j ) {
				for ( int ell = i + 1; ell <= M; ++ ell ) {
					--(row_index[ell]);
				}

				for ( int ell = k; ell <= row_index[M]; ++ell ) {
					column_index[ell] = column_index[ell + 1];
					off_diagonal[ell] = off_diagonal[ell + 1];
				}

				return;
			}

			if ( column_index[k] > j ) {
				return;
			}
		}
	}
}

/***************************************************
 * Set all entries in the matrix to 0.
 *
 * Run time:   O( M )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::clear() {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] = 0;
	}

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

/***************************************************
 * Return the value A(i, j)
 *
 * Run time:   Diagonal:      O( 1 )
 *             Off-diagonal:  O( n/M )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
double Matrix<M, N>::operator() ( int i, int j ) const {
	assert( i >= -M && j >= -N && i < M && j < N );

	i += ( i < 0 ) ? M : 0;
	j += ( j < 0 ) ? N : 0;

	if ( i == j ) {
		return diagonal[i];
	}

	// For off-diagonal entries,
	// the entries of row i are located in the
	//        row_index[i], ..., row_index[i + 1] - 1
	// The entries in this block are stored in column order

	for ( int k = row_index[i]; k < row_index[i + 1]; ++k ) {
		// If the column index matches, return the corresponding
		// entry in the array 'off_diagonal'.

		if ( column_index[k] == j ) {
			return off_diagonal[k];
		}

		// If we have passed column j, return 0.

		if ( column_index[k] > j ) {
			return 0;
		}
	}

	// The column entry is greater than the
	// last non-zero entry of row i.

	return 0;
}