#ifndef _sparse_inner_prod_h_
#define _sparse_inner_prod_h_



#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/vector_of_vector.hpp>

template <class T, std::size_t IB, class IA, class TA>
BOOST_UBLAS_INLINE
T sparse_vector_inner_prod(const ublas::compressed_vector<T, IB, IA, TA> &lhs,
		                     const ublas::compressed_vector<T, IB, IA, TA> &rhs)
{
	typedef T result_type;
	typedef ublas::compressed_vector<T, IB, IA, TA> vector_type;
	typedef IA index_array_type;
	typedef TA value_array_type;

	const index_array_type &index1 = lhs.index_data();
	const value_array_type &value1 = lhs.value_data();
	size_t nnz1(lhs.nnz());

	const index_array_type &index2 = rhs.index_data();
	const value_array_type &value2 = rhs.value_data();
	size_t nnz2(rhs.nnz());

	result_type retval(0);
	if((nnz1 != 0) && (nnz2 != 0))
	{
		//
		// Multiply the existing elements in the first vector by any common elements
		// in the second vector.
		//
		size_t last_pos(0);
		for(size_t i(0); i<nnz1; ++i)
		{
			//
			// We resume our search from the position at which we found our last match.
			//
			for(size_t j(last_pos); j<nnz2; ++j)
			{
				//
				// If they are the same index, perform the element-element multiplication.
				//
				if(index1[i] == index2[j])
				{
					retval += value1[i] * value2[j];
					last_pos = j+1;
					break;
				}
				//
				// If the index in the lhs is less than the current index of the rhs we can
				// stop searching...its not here.
				//
				else if(index1[i] < index2[j])
				{
					last_pos = j;
					break;
				}
			}
		}
	}

	return(retval);
}

template <class T, class L, class V, class M>
BOOST_UBLAS_INLINE
void gvov_ata_prod(const ublas::generalized_vector_of_vector<T, L, V> &A, ublas::matrix_expression<M> &result)
{
	typedef T value_type;

	value_type val;
	for(size_t i(0); i<A.size2(); ++i)
	{
		//
		// Take advantage of the symmetric structure of A'A.
		//
		for(size_t j(i); j<A.size2(); ++j)
		{
			//
			// Perform the inner product of the columns inserting into the result iff the
			// value is not zero.
			//
			val = sparse_vector_inner_prod((A.data())[i], (A.data())[j]);
			if(fabs(val) > 1.0e-07)
			{
				(result())(i, j) = val;
				(result())(j, i) = val;
			}
		}
	}
}

#endif