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

/***************************************************
 * Increases the arrays for storing the off-
 * diagonal entries.
 *
 * Run time:   O( n )
 * Memory:     O( cap )
 ***************************************************/

template <int M, int N>
void Matrix<M, N>::increase_capacity( int cap ) {
	assert( cap > row_index[M] );

	// Store the current values in temporary variables

	int    *tmp_column_index = column_index;
	double *tmp_off_diagonal = off_diagonal;

	// Allocate new memory

	capacity = cap;
	column_index = new int[capacity];
	off_diagonal = new double[capacity];

	// Copy all off-diagonal entries

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

	// Delete the old arrays

	delete [] tmp_column_index;
	delete [] tmp_off_diagonal;
}

/***************************************************
 * Return the number of off-diagonal entries which
 * will result in non-zero values if the matrices
 * are either added (true) or subtracted (false).
 *
 * Run time:   O( na + nb )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
int Matrix<M, N>::count_nonzero_entries( Matrix const &A, Matrix const &B, bool sum ) {
	int count = 0;

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

		while ( aj < A.row_index[i + 1] && bj < B.row_index[i + 1] ) {
			if ( A.column_index[aj] == B.column_index[bj] ) {
				if (
					( sum && (A.off_diagonal[aj] != -B.off_diagonal[bj] )) ||
					(!sum && (A.off_diagonal[aj] !=  B.off_diagonal[bj] ))
				) {
					++count;
				}

				++aj;
				++bj;
			} else if ( A.column_index[aj] < B.column_index[bj] ) {
				++count;
				++aj;
			} else {
				++count;
				++bj;
			}
		}

		count += (A.row_index[i + 1] - aj + B.row_index[i + 1] - bj);
	}

	return count;
}