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

template<int M, int N>
template<Orientation D>
Matrix<M, N>  &Matrix<M, N>::operator+= ( Vector<(M < N) ? M : N, D> const &u ) {
	for ( int i = 0; i < minMN; ++i ) {
		diagonal[i] += u.array[i];
	}

	return *this;
}

template<int M, int N>
template<Orientation D>
Matrix<M, N>  &Matrix<M, N>::operator-= ( Vector<(M < N) ? M : N, D> const &u ) {
	for ( int i = 0; i < N; ++i ) {
		diagonal[i] -= u.array[i];
	}

	return *this;
}

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

template<int M, int N, Orientation D>
Matrix<M, N> operator+ ( Matrix<M, N> const &A, Vector<(M < N) ? M : N, D> const &u ) {
	Matrix<M, N> B = A.copy();

	B += u;

	return B;
}

template<int M, int N, Orientation D>
Matrix<M, N> operator+ ( Vector<(M < N) ? M : N, D> const &u, Matrix<M, N> const &A ) {
	Matrix<M, N> B = A.copy();

	B += u;

	return B;
}

template<int M, int N, Orientation D>
Matrix<M, N> operator- ( Matrix<M, N> const &A, Vector<(M < N) ? M : N, D> const &u ) {
	Matrix<M, N> B = A.copy();

	B -= u;

	return B;
}

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

	B += u;

	return B;
}

/***************************************************
 * Calculate the matrix-vector product v = A*u
 *
 * Run time: O(n + N)
 ***************************************************/

template<int N>
Vector<N, COLUMN> &operator*= ( Matrix<N, N> const &A, Vector<N, COLUMN> &u ) {
	double array[N];

	// First calculate A(i,i)*u(i)

	for ( int i = 0; i < N; ++i ) {
		array[i] = A.diagonal[i] * u.array[i];
	}

	// Next, for each off-diagonal entry v(i) += A(i,j)*u(j)

	for ( int i = 0; i < N; ++i ) {
		for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
			array[i] += A.off_diagonal[j] * u.array[A.column_index[j]];
		}
	}

	for ( int i = 0; i < N; ++i ) {
		u.array[i] = array[i];
	}

	return u;
}

template<int M, int N>
Vector<M, COLUMN> operator* ( Matrix<M, N> const &A, Vector<N, COLUMN> const &u ) {
	Vector<M, COLUMN> v;

	// First calculate v(i) = A(i,i)*u(i)

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

	for ( int i = Matrix<M,N>::minMN; i < M; ++i ) {
		v.array[i] = 0.0;
	}

	// Next, for each off-diagonal entry v(i) += A(i,j)*u(j)

	for ( int i = 0; i < M; ++i ) {
		for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
			v.array[i] += A.off_diagonal[j] * u.array[A.column_index[j]];
		}
	}

	return v;
}

template<int N>
Vector<N, ROW> operator*= ( Vector<N, ROW> &u, Matrix<N, N> const &A ) {
	double array[N];

	// First calculate v(i) = A(i,i)*u(i)

	for ( int i = 0; i < N; ++i ) {
		array[i] = A.diagonal[i] * u.array[i];
	}

	// Next, for each off-diagonal entry v(i) += A(i,j)*u(j)

	for ( int i = 0; i < N; ++i ) {
		for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
			array[A.column_index[j]] += u.array[i] * A.off_diagonal[j];
		}
	}

	for ( int i = 0; i < N; ++i ) {
		u.array[i] = array[i];
	}

	return u;
}

template<int M, int N>
Vector<N, ROW> operator* ( Vector<M, ROW> const &u, Matrix<M, N> const &A ) {
	Vector<M, ROW> v;

	// First calculate v(i) = A(i,i)*u(i)

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

	for ( int i = Matrix<M,N>::minMN; i < N; ++i ) {
		v.array[i] = 0.0;
	}

	// Next, for each off-diagonal entry v(i) += A(i,j)*u(j)

	for ( int i = 0; i < M; ++i ) {
		for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
			v.array[A.column_index[j]] += u.array[i] * A.off_diagonal[j];
		}
	}

	return v;
}

/***************************************************
 * Calculate A*v by ignoring all off-diagonal
 * entries of A.
 *
 * Run time:   O( M )             DIAGONAL
 *             O( M + n )         All others
 * Memory:     O( M )
 ***************************************************/

template<int M, int N>
Vector<M, COLUMN> matrix_vector_product( Matrix<M, N> const &A, Vector<N, COLUMN> const &u, Product prod = FULL ) {
	Vector<M, COLUMN> v;

	if ( prod == OFF_DIAGONAL || prod == STRICT_UPPER || prod == STRICT_LOWER ) {
		for ( int i = 0; i < Matrix<M,N>::minMN; ++i ) {
			v.array[i] = 0.0;
		}
	} else {
		for ( int i = 0; i < Matrix<M,N>::minMN; ++i ) {
			v.array[i] = A.diagonal[i] * u.array[i];
		}
	}

	for ( int i = Matrix<M,N>::minMN; i < M; ++i ) {
		v.array[i] = 0.0;
	}

	if ( prod == FULL || prod == OFF_DIAGONAL ) {
		// For each off-diagonal entry v(i) += A(i,j)*u(j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
				v.array[i] += A.off_diagonal[j] * u.array[A.column_index[j]];
			}
		}
	} else if ( prod == LOWER || prod == STRICT_LOWER ) {
		// For each off-diagonal entry v(i) += A(i,j)*u(j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i]; j < A.row_index[i + 1] && A.column_index[j] < i; ++j ) {
				v.array[i] += A.off_diagonal[j] * u.array[A.column_index[j]];
			}
		}
	} else if ( prod == UPPER || prod == STRICT_UPPER ) {
		// For each off-diagonal entry v(i) += A(i,j)*u(j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i + 1] - 1; j >= A.row_index[i] && A.column_index[j] > i; --j ) {
				v.array[i] += A.off_diagonal[j] * u.array[A.column_index[j]];
			}
		}
	}


	return v;
}

/***************************************************
 * Calculate v*A by ignoring all off-diagonal
 * entries of A.
 *
 * Run time:   O( N )             DIAGONAL
 *             O( N + n )         All others
 * Memory:     O( N )
 ***************************************************/

template<int M, int N>
Vector<N, ROW> vector_matrix_product( Vector<M, ROW> const &u, Matrix<M, N> const &A, Product prod = FULL ) {
	Vector<N, ROW> v;

	if ( prod == OFF_DIAGONAL || prod == STRICT_UPPER || prod == STRICT_LOWER ) {
		// Set v(i) = 0

		for ( int i = 0; i < Matrix<M,N>::minMN; ++i ) {
			v.array[i] = 0.0;
		}
	} else {
		for ( int i = 0; i < Matrix<M,N>::minMN; ++i ) {
			v.array[i] = A.diagonal[i] * u.array[i];
		}
	}

	for ( int i = Matrix<M,N>::minMN; i < N; ++i ) {
		v.array[i] = 0.0;
	}

	if ( prod == FULL || prod == OFF_DIAGONAL ) {
		// For each off-diagonal entry v(i) += u(j)*A(i,j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i]; j < A.row_index[i + 1]; ++j ) {
				v.array[A.column_index[j]] += u.array[i] * A.off_diagonal[j];
			}
		}
	} else if ( prod == LOWER || prod == STRICT_LOWER ) {
		// For each off-diagonal entry v(i) += u(j)*A(i,j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i]; j < A.row_index[i + 1] && A.column_index[j] < i; ++j ) {
				v.array[A.column_index[j]] += u.array[i] * A.off_diagonal[j];
			}
		}
	} else if ( prod == UPPER || prod == STRICT_UPPER ) {
		// For each off-diagonal entry v(i) += u(j)*A(i,j)

		for ( int i = 0; i < M; ++i ) {
			for ( int j = A.row_index[i + 1] - 1; j >= A.row_index[i] && A.column_index[j] > i; --j ) {
				v.array[A.column_index[j]] += u.array[i] * A.off_diagonal[j];
			}
		}
	}

	return v;
}