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

/***************************************************
 * Add the argument 'x' to all the entries
 * on the diagonal.
 *
 * Run time:   O( min( M, N ) )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
Matrix<M, N> &Matrix<M, N>::operator+= ( double x ) {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] += x;
	}

	return *this;
}

/***************************************************
 * Subtract the argument 'x' from all the entries
 * on the diagonal.
 *
 * Run time:   O( min( M, N ) )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
Matrix<M, N> &Matrix<M, N>::operator-= ( double x ) {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] -= x;
	}

	return *this;
}

/***************************************************
 * Multiply all the entries in the matrix by 'x'.
 * Note that for an non-zero off-diagonal entry
 * A(i,j), * if x*A(i,j) == 0, it must be removed.
 *
 * Run time:   O( M + n )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
Matrix<M, N> &Matrix<M, N>::operator*= ( double x ) {
	if ( x == 0.0 ) {
		clear();

		return *this;
	}

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

	int j = 0;
	int k = 0;

	for ( int i = 1; i <= M; ++i ) {
		for ( ; j < row_index[i]; ++j ) {
			double product = off_diagonal[j]*x;

			if ( product != 0.0 ) {
				column_index[k] = column_index[j];
				off_diagonal[k] = product;

				++k;
			}
		}

		row_index[i] = k;
	}

	return *this;
}

/***************************************************
 * Divide all the entries in the matrix by 'x'.
 * Note that for an non-zero off-diagonal entry
 * A(i,j), * if x*A(i,j) == 0, it must be removed.
 *
 * Run time:   O( M + n )
 * Memory:     O( 1 )
 ***************************************************/

template <int M, int N>
Matrix<M, N> &Matrix<M, N>::operator/= ( double x ) {
	if ( x == 0.0 ) {
		throw division_by_zero();
	}

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

	int j = 0;
	int k = 0;

	for ( int i = 1; i <= M; ++i ) {
		for ( ; j < row_index[i]; ++j ) {
			double ratio = off_diagonal[j]/x;

			if ( ratio != 0.0 ) {
				column_index[k] = column_index[j];
				off_diagonal[k] = ratio;

				++k;
			}
		}

		row_index[i] = k;
	}

	return *this;
}

/***************************************************
 * *********************************************** *
 * *                                             * *
 * *   FRIENDS                                   * *
 * *                                             * *
 * *********************************************** *
 ***************************************************/

/***************************************************
 * Return the matrix A with 'x' added to all the
 * entries on the diagonal.
 *  - A + x*eye(M, N)
 *  - See operator+= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator+ ( Matrix<M, N> const &A, double x ) {
	Matrix<M, N> B( A );

	B += x;

	return B;
}

/***************************************************
 * Return the matrix A with 'x' subtracted from
 * all the entries on the diagonal.
 *  - A - x*eye(M, N)
 *  - See operator-= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator- ( Matrix<M, N> const &A, double x ) {
	Matrix<M, N> B( A );

	B -= x;

	return B;
}

/***************************************************
 * Return the matrix A with all entries multiplied
 * by 'x'.
 *  - x*A
 *  - See operator*= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator* ( Matrix<M, N> const &A, double x ) {
	Matrix<M, N> B( A );

	B *= x;

	return B;
}

/***************************************************
 * Return the matrix A with all entries divided
 * by 'x'.
 *  - A/x
 *  - See operator/= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator/ ( Matrix<M, N> const &u, double x ) {
	Matrix<M, N> w( u );

	w /= x;

	return w;
}

/***************************************************
 * Return the matrix A with 'x' added to all the
 * entries on the diagonal.
 *  - A + x*eye(M, N)
 *  - See operator+= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator+ ( double x, Matrix<M, N> const &A ) {
	Matrix<M, N> B( A );

	B += x;

	return B;
}

/***************************************************
 * Return the matrix of x - A(i,i) on the diagonal
 * and -A(i,j) on the off-diagonal.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator- ( double x, Matrix<M, N> const &A ) {
	Matrix<M, N> B( A.row_index[M] );

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

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

	return B;
}

/***************************************************
 * Return the matrix A with all entries multiplied
 * by 'x'.
 *  - x*A
 *  - See operator*= ( double ) above.
 *
 * Run time:   O( min( M, N ) + na )
 * Memory:     O( min( M, N ) + ca )
 ***************************************************/

template <int M, int N>
Matrix<M, N>  operator* ( double x, Matrix<M, N> const &A ) {
	Matrix<M, N> B( A );

	B *= x;

	return B;
}