// -*- c++ -*- //
/** \file main.cpp  \brief main routines */
/***************************************************************************
 -   begin                : 2003-09-05
 -   copyright            : (C) 2003 by Gunter Winkler
 -   email                : guwi17@gmx.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <iostream>

#include <boost/timer.hpp>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>


#ifndef NDEBUG
const size_t default_size = 100;
#else
const size_t default_size = 1000;
#endif

const size_t subsize = 3;

using namespace boost::numeric::ublas ;

using std::cout;
using std::endl;

namespace boost { namespace numeric { namespace ublas {

    /** \brief compute residual r = b - Ax
     *
     */
    template<class V, class T1, class IA1, class TA1, class E2>
    void 
    residual ( const compressed_matrix<T1, row_major, 0, IA1, TA1> & A,
               const vector_expression<E2> &x,
               const vector_expression<E2> &b,
               V &r, row_major_tag) {

        BOOST_UBLAS_CHECK( x().size() >= A.size2(), external_logic() );
        BOOST_UBLAS_CHECK( b().size() >= A.size1(), external_logic() );
        BOOST_UBLAS_CHECK( r.size() >= A.size1(), external_logic() );

        typedef typename V::size_type size_type;
        typedef typename V::value_type value_type;
        
        for (size_type i = 0; i < A.size1 (); ++ i) {
            size_type begin = A.index1_data () [i];
            size_type end = A.index1_data () [i + 1];
            value_type t (b() (i));
            for (size_type j = begin; j < end; ++ j)
                t -= A.value_data () [j] * x () (A.index2_data () [j]);
            r (i) = t;
        }
        return;
    }

}}}

template <class MAT>
void matrix_fill(MAT& A)
{
  typedef typename MAT::size_type  size_type;

  for (size_type i=0; i<A.size1(); ++i)	{
	if (i>=1) { A.push_back(i,i-1,-1.0); }
	A.push_back(i,i,4);
	if (i+1<A.size2()) { A.push_back(i,i+1,-1.0); }
  }
}

int main(int argc, char *argv[])
{

  typedef size_t  size_type;

  size_type size = default_size;

  if (argc > 1)
	size = ::atoi (argv [1]);

  cout << "size " << size << endl;

  compressed_matrix<double, row_major> A(size,size);
  matrix_fill(A);

  vector<double>                       b( scalar_vector<double>(size, 2.0) );
  b(0)=3; 
  b(size-1)=3;

  vector<double>                       x( scalar_vector<double>(size, 0.99) );

  {
	boost::timer t;
	t.restart();
	vector<double>                       r(size);
	r = prod(A, x) - b;
	cout << "r=Ax-b : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
  }

  {
	boost::timer t;
	t.restart();
	vector<double>                       r(size);
	noalias(r) = b - prod(A, x);
	cout << "noalias(r)=b-Ax : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
  }

  {
	boost::timer t;
	t.restart();
	vector<double>                       r(b - prod(A, x));
	cout << "vector<double> r(b - prod(A, x)) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
  }

  {
	boost::timer t;
	t.restart();
	vector<double>                       r(b);
	axpy_prod(A,-x,r,false);
	cout << "axpy_prod(A,-x,r,false) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
  }

  {
	boost::timer t;
	t.restart();
	vector<double>                       r(size);
	residual(A,x,b,r,row_major_tag());
	cout << "residual(A,x,b,r,row_major_tag()) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
  }

  return EXIT_SUCCESS;

}