//          Copyright Gunter Winkler 2004 - 2007.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// authors: Gunter Winkler <guwi17 at gmx dot de>
//          Konstantin Kutzkow

#include "boost/numeric/ublas/vector_expression.hpp"

using boost::numeric::ublas::vector_expression;
using boost::numeric::ublas::vector_container;

#include <algorithm>

/** \brief The class representing a sparse vector denoted as v.

 * \param index: a dense vector of size nnz containing the indexes of
 * the original nonzero elements in v.

 * \param value: a dense vector of size nnz containing the real of the
 * original nonzero elements in v.
 */
class working_row 
{
private:
  typedef double      value_type;

  typedef std::vector<unsigned int>  index_array_type;
  typedef std::vector<value_type>    value_array_type;

  index_array_type index;
  value_array_type value;

public:
  working_row() :
    index(), value()
  { }

  working_row(unsigned int size) :
    index(), value(size, 0.0)
  { }

  template < class V >
  working_row &
  assign(vector_expression<V> const & ve) 
  { 
    index.reserve(ve().nnz());
    value.resize(ve().size());
    std::fill( value.begin(), value.end(), 0.0 );
    copy_from_sparse(ve());
    return *this;
  }

  // assign temporary wr. wr must not be used afterwards.
  void assign_temporary(working_row & wr) 
  {
    swap(index, wr.index);
    swap(value, wr.value);
  }

  void resize(unsigned int const new_size) 
  {
    value.resize(new_size);
    clear();
  }

  void clear() 
  {
    index.resize(0);
    std::fill( value.begin(), value.end(), 0.0 );
  }

  template < class V >
  void copy_from_sparse(vector_expression<V> const & ve) 
  {
    // precondition: value array is set to zero.
    index.resize(ve().nnz());
    assert( ve().size() == value.size() );
    typename V::const_iterator v_iter = ve().begin();
    typename V::const_iterator v_end  = ve().end();
    for (unsigned int i=0; v_iter != v_end; ++v_iter, ++i) {
      index[i] = v_iter.index();
      value[ v_iter.index() ] = *v_iter;
    }
    // postcondition: wr == ve
  }

  template < class V >
  void copy_to_sparse(vector_container<V> & ve) const 
  {
    ve().clear();
    typename index_array_type::const_iterator ia_iter = index.begin();
    typename index_array_type::const_iterator ia_end  = index.end();
    for ( ; ia_iter != ia_end; ++ia_iter ) {
      ve().insert_element( *ia_iter, value[ *ia_iter ] );
    }
    // postcondition: ve == wr
  }

  template < class V >
  void zero_and_copy_to_sparse(vector_container<V> & ve) 
  {
    ve().clear();
    typename index_array_type::const_iterator ia_iter = index.begin();
    typename index_array_type::const_iterator ia_end  = index.end();
    for ( ; ia_iter != ia_end; ++ia_iter ) {
      ve().insert_element( *ia_iter, value[ *ia_iter ] );
      value[ *ia_iter ] = 0;
    }
    // postcondition: value[i] == 0.0 for 0<=i<value.size()
  }

  template < class V >
  void add_sparse(typename V::value_type const & alpha, 
                  vector_expression<V> const & ve)
  {
    assert( ve().size() == value.size() );
    typename V::const_iterator v_iter = ve().begin();
    typename V::const_iterator v_end  = ve().end();
    if ( alpha == 0.0 ) return;
    for ( ; v_iter != v_end; ++v_iter) {
      value_type t = value[ v_iter.index() ];
      if ( 0.0 == t ) {
        t = alpha * (*v_iter);
        if (0.0 != t) index.push_back( v_iter.index() );
      } else {
        t += alpha * (*v_iter);
      }
      value[ v_iter.index() ] = t;
    }
  }

};

/** \brief The function for the converting of a sparse vector to
 * element of the class working_row.
 */
template <class Vector>
void scatter(Vector& x, working_row& wr)
{
  wr.assign(x);
}

/** \brief The function for the converting of an element of the class
 * working_row to a sparse vector 
 */
template <class Vector>
void gather(Vector& x, working_row const & wr)
{
  wr.copy_to_sparse(x);
}

/** \brief The function for the adding of a sparse vector to an element of the class working_row
 */
template <class Vector>
void sp_axpy(double const & a, Vector const & y, working_row& wr)
{
  wr.add_sparse(a,y);
}