// -*- c++ -*-
//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile$  $Revision$  $State$ 
//     $Author$  $Date$ 
//     $Locker$ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//   This is the file contain the implementations of basic operations.
//   It provides uBLAS implementations of those operations. Therefore,
//   It requires uBLAS package. The tested BOOST release is 1.29.0
//                                                                           
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
//=======================================================================
// Copyright (C) 1997-2003
// Authors: Andrew Lumsdaine <lums@osl.iu.edu> 
//          Lie-Quan Lee     <llee@osl.iu.edu>
//          Gunter Winkler   <guwi17@gmx.de>
//
// This file is part of the Iterative Template Library
//
// You should have received a copy of the License Agreement for the
// Iterative Template Library along with the software;  see the
// file LICENSE.  
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//=======================================================================
//---------------------------------------------------------------------------
//                                                                           
// REVISION HISTORY:                                                         
//                                                                           
// $Log$
//                                                                           
//===========================================================================

#ifndef ITL_UBLAS_INTERFACE_H
#define ITL_UBLAS_INTERFACE_H

/*
  This is to provide an interface for ITL to use boost-ublas. 
  
  Matrix: ublas matrices
  Vector: ublas vector
  Preconditioners: all preconditioners in itl/preconditioners
 */

#include "itl/itl_config.h"
#include "itl/itl_tags.h"
#include "itl/number_traits.h"

#include "boost/numeric/ublas/matrix.hpp"
#include "boost/numeric/ublas/matrix_sparse.hpp"
#include "boost/numeric/ublas/vector.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/io.hpp"
#include "boost/numeric/ublas/operation.hpp"

namespace itl {
  
  //: The vector type used inside of the ITL routines for work space
  template <class Vec>
  struct itl_traits {
    typedef referencing_object_tag   vector_category;
    typedef typename Vec::value_type     value_type;
    typedef typename Vec::size_type      size_type;
  };

  template <class Matrix>
  struct itl_matrix_traits
  {
	typedef typename Matrix::orientation_category orientation;
    typedef typename Matrix::value_type           value_type;
    typedef typename Matrix::size_type            size_type;

	typedef typename boost::numeric::ublas::row_major_tag    row_major;
	typedef typename boost::numeric::ublas::column_major_tag column_major;

	typedef typename Matrix::iterator1            row_iterator;
	typedef typename Matrix::iterator2            column_iterator;
	typedef typename Matrix::const_iterator1      const_row_iterator;
	typedef typename Matrix::const_iterator2      const_column_iterator;
  };
  /*
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::const_column_iterator 
  column_begin(const Matrix& A)
  { return A.begin2(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::const_row_iterator 
  row_begin(const Matrix& A)
  { return A.begin1(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::const_column_iterator 
  column_end(const Matrix& A)
  { return A.end2(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::const_row_iterator 
  row_end(const Matrix& A)
  { return A.end1(); }
  */
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::column_iterator 
  column_begin(Matrix& A)
  { return A.begin2(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::row_iterator 
  row_begin(Matrix& A)
  { return A.begin1(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::column_iterator 
  column_end(Matrix& A)
  { return A.end2(); }
  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::row_iterator 
  row_end(Matrix& A)
  { return A.end1(); }

  template <class Iterator>
  inline typename Iterator::size_type row_index(const Iterator& it)
  { return it.index1(); }
  template <class Iterator>
  inline typename Iterator::size_type column_index(const Iterator& it)
  { return it.index2(); }

  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::size_type nrows(const Matrix& A)
  { return A.size1(); }

  template <class Matrix>
  inline typename itl_matrix_traits<Matrix>::size_type ncols(const Matrix& A)
  { return A.size2(); }

  template <class Matrix>
  inline typename Matrix::size_type nnz(const Matrix& A)
  { return A.non_zeros(); }

  template <class Vec>
  inline 
  typename itl::number_traits< typename Vec::value_type >::magnitude_type
  two_norm(const Vec& v) {
    return boost::numeric::ublas::norm_2(v);
  }

  //deal with the case when b is a handle, 
  template <class VecA, class VecB>
  inline void copy(const boost::numeric::ublas::vector_expression<VecA>& a, VecB& b) {
	// assume: a is no alias of b
	b.assign( a );
  }
  
  template <class Matrix, class VecX, class VecY>
  inline void trans_mult(const boost::numeric::ublas::matrix_expression<Matrix>& A, 
						 const boost::numeric::ublas::vector_expression<VecX>& x, 
						 VecY& y) {
	y = boost::numeric::ublas::prod( boost::numeric::ublas::trans(A), x);
  }

  template <class Matrix, class VecX, class VecY>
  inline void mult(const boost::numeric::ublas::matrix_expression<Matrix>& A, 
				   const boost::numeric::ublas::vector_expression<VecX>& x, 
				   VecY& y) {
	//	y = boost::numeric::ublas::prod (A, x);
	boost::numeric::ublas::axpy_prod (A, x, y, true);
  }

  template <class Matrix, class VecX, class VecY, class VecZ>
  inline void mult(const boost::numeric::ublas::matrix_expression<Matrix>& A, 
				   const boost::numeric::ublas::vector_expression<VecX>& x, 
				   const boost::numeric::ublas::vector_expression<VecY>& y, 
				   VecZ& z) 
  {
	copy(y, z);
	boost::numeric::ublas::axpy_prod(A, x, z, false);
	//z = y + boost::numeric::ublas::prod(A, x);
  }

  template <class VecA, class VecB>
  inline typename VecA::value_type dot(const boost::numeric::ublas::vector_expression<VecA>& a, 
									   const boost::numeric::ublas::vector_expression<VecB>& b) 
  {
    return boost::numeric::ublas::inner_prod(a, b);
  }
  
  template <class VecA, class VecB>
  inline typename VecA::value_type  dot_conj(const boost::numeric::ublas::vector_expression<VecA>& a, 
											 const boost::numeric::ublas::vector_expression<VecB>& b) 
  {
    return boost::numeric::ublas::inner_prod(a, boost::numeric::ublas::conj(b));
  }
  

  template <class VecX, class VecY>
  inline void add(const boost::numeric::ublas::vector_expression<VecX>& x, 
				  VecY& y) {
	// I we knew that x is no alias of y, then y.plus_assign( x ) would be faster
	y += x;
  }

  template <class VecX, class VecY, class VecZ>
  inline void add(const boost::numeric::ublas::vector_expression<VecX>& x, 
				  const boost::numeric::ublas::vector_expression<VecY>& y, 
				  VecZ& z) {
	z = x + y;
  }

  template <class VecX, class VecY, class VecZ, class VecR>
  inline void add(const boost::numeric::ublas::vector_expression<VecX>& x, 
				  const boost::numeric::ublas::vector_expression<VecY>& y, 
				  const boost::numeric::ublas::vector_expression<VecZ>& z, 
				  VecR& r) {
	r = x+y+z;
  }

  // (v * t) [i] = v [i] * t
  template<class E1, class T2>
  typename boost::numeric::ublas::vector_binary_scalar2_traits<
	E1, 
	const T2, 
	boost::numeric::ublas::scalar_multiplies<
	  typename E1::value_type, 
	  T2
	  >,
	boost::numeric::ublas::scalar_reference<const T2>
	>::result_type
  scaled(const boost::numeric::ublas::vector_expression<E1> &v, const T2 &t)
  {
	return v*t;
  }

  template <class Vec, class T>
  inline void scale(const Vec& v, T t) {
	v *= t;
  }

  template <class VecA>
  inline void ele_mult(const boost::numeric::ublas::vector_expression<VecA>& x, 
					   const boost::numeric::ublas::vector_expression<VecA>& y, 
					   VecA& z) {
	// assume: size(x) == size(y) == size(z)
	typename VecA::size_type i = 0;
	typename VecA::size_type n = size(x);
	for (i=0; i<n; ++i) z[i] = x()(i) * y()(i);
  }

  template <class Vector>
  inline typename Vector::size_type size(const boost::numeric::ublas::vector_expression<Vector>& x) 
  {
    return x ().size();
  }


  template <class Vector, class Size>
  inline void resize(Vector& x, const Size& sz) 
  {
    x.resize(sz);
  }


  //used inside of GCR and GMRES algorithm

  template <class T>
  class internal_matrix_traits {
  public:
    typedef typename boost::numeric::ublas::matrix<T>  Matrix;
  };
  template <class T>
  class internal_vector_traits {
  public:
    typedef typename boost::numeric::ublas::vector<T>  Vector;
  };

  
  template <class Hessenberg, class Vec>
  inline void upper_tri_solve(const Hessenberg& hh, Vec& rs, int i) 
  {
	using boost::numeric::ublas::matrix;
	using boost::numeric::ublas::matrix_range;
	using boost::numeric::ublas::range;
	using boost::numeric::ublas::solve;

	inplace_solve( matrix_range<Hessenberg>(A,range(0,i),range(0,i)), 
				   rs, 
				   boost::numeric::ublas::upper_tag() );
  }

  // FIXME: adapt tri-solves to current syntax

  template <class Matrix, class Vec>
  inline void tri_solve(const Matrix& M, Vec& rhs, upper_tag) 
  {
	inplace_solve(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::upper>(M), rhs, boost::numeric::ublas::upper_tag());
  }

  template <class Matrix, class Vec>
  inline void tri_solve(const Matrix& M, Vec& rhs, unit_upper_tag) 
  {
	inplace_solve(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::unit_upper>(M), rhs, boost::numeric::ublas::unit_upper_tag());
  }

  template <class Matrix, class Vec>
  inline void tri_solve(const Matrix& M, Vec& rhs, lower_tag) 
  {
	// inplace_solve(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::lower>(M), rhs, boost::numeric::ublas::lower_tag());
	inplace_solve(M, rhs, boost::numeric::ublas::lower_tag());
  }

  template <class Matrix, class Vec>
  inline void tri_solve(const Matrix& M, Vec& rhs, unit_lower_tag) 
  {
	inplace_solve(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::unit_lower>(M), rhs, boost::numeric::ublas::unit_lower_tag());
  }

  template <class Matrix, class Vec>
  inline void tri_trans_solve(const Matrix& M, Vec& rhs, upper_tag) 
  {
	inplace_solve(trans(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::upper>(M)), rhs, boost::numeric::ublas::lower_tag());
	// inplace_solve(trans(M), rhs, boost::numeric::ublas::lower_tag() );
  }

  template <class Matrix, class Vec>
  inline void tri_trans_solve(const Matrix& M, Vec& rhs, unit_upper_tag) 
  {
	inplace_solve(trans(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::unit_upper>(M)), rhs, boost::numeric::ublas::unit_lower_tag());
	// inplace_solve(trans(M), rhs, boost::numeric::ublas::unit_lower_tag() );
  }

  template <class Matrix, class Vec>
  inline void tri_trans_solve(const Matrix& M, Vec& rhs, lower_tag) 
  {
	inplace_solve(trans(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::lower>(M)), rhs, boost::numeric::ublas::upper_tag());
	// inplace_solve(trans(M), rhs, boost::numeric::ublas::upper_tag() );
  }

  template <class Matrix, class Vec>
  inline void tri_trans_solve(const Matrix& M, Vec& rhs, unit_lower_tag) 
  {
	inplace_solve(trans(typename boost::numeric::ublas::triangular_adaptor<const Matrix, boost::numeric::ublas::unit_lower>(M)), rhs, boost::numeric::ublas::unit_upper_tag());
	// inplace_solve(trans(M), rhs, boost::numeric::ublas::unit_upper_tag() );
  }

  template <class DATA>
  void debug_print(const char* text, const DATA& data)
  {
#ifndef NDEBUG
	std::cout << text << data << endl;
#endif
	return;
  }

}

//for classical gram schmidt 
//#include "itl/interface/detail/mtl_classical_gram_schmidt.h"

#endif /*ITL_MTL_INTERFACE_H*/