Re: ITL-Boost interface

// -*- 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 <[email protected]> 
//          Lie-Quan Lee     <[email protected]>
//          Gunter Winkler   <[email protected]>
//
// 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/vector.hpp"
#include "boost/numeric/ublas/triangular.hpp"
#include "boost/numeric/ublas/io.hpp"

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

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

  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 VecA& a, VecB& b) {
	b = a;
  }
  
  template <class Matrix, class VecX, class VecY>
  inline void mult(const Matrix& A, const VecX& x, VecY& y) {
	y = boost::numeric::ublas::prod (A, x);
  }

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

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

  template <class VecX, class VecY>
  inline void add(const VecX& x, VecY& y) {
	y += x;
  }

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

  template <class VecX, class VecY, class VecZ, class VecR>
  inline void add(const VecX& x, const VecY& y, const 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, T2, boost::numeric::ublas::scalar_multiplies<typename E1::value_type, 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 VecA& x, const 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 Vector& x) 
  {
    return x.size();
  }


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


  template <class Matrix, class VecX, class VecY>
  inline void trans_mult(const Matrix& A, const VecX& x, const VecY& y) {
	y = boost::numeric::ublas::prod( boost::numeric::ublas::trans(A), x);
  }

  //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;

	rs = solve( matrix_range<Hessenberg>(A,range(0,i),range(0,i)), 
				rs, 
				boost::numeric::ublas::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*/

#ifndef ITL_DIAGONAL_PRECOND_H
#define ITL_DIAGONAL_PRECOND_H

/*
 * changes for uBlas:
 *   new interface function ele_mult(x,y,z):  z = x .* y
 *   new interface type internal_vector_traits::Vector
 */

namespace itl {

  template <class V>
  struct D_precond {
    D_precond(const V& v) : diag(v) {}
    V diag;
  };


 template <class V, class VectorX, class VectorY>
  inline void solve(const D_precond<V>& M, const VectorX& x, 
					VectorY& y) {
    ele_mult(M.diag, x, y);
  }

  template <class V, class VectorX, class VectorY>
  inline void trans_solve(const D_precond<V>& M, const VectorX& x, 
						  VectorY& y) {
    ele_mult(M.diag, x, y);
  }

  template <class Matrix>
  struct diagonal_precond {

    typedef typename itl_traits<Matrix>::value_type     T;
	typedef typename internal_vector_traits<T>::Vector  Vector;

    diagonal_precond(const Matrix& A): diag( nrows(A) ) 
	{
      typedef typename itl_traits<Matrix>::size_type size_type;
	  size_type n = nrows(A);
      for (size_type i=0; i < n; ++i) diag[i] = T(1) / A(i, i);
    }

    typedef D_precond< Vector > Precond;
    inline Precond operator()() { return Precond(diag); }

    Vector   diag;
  };


}


#endif //ITL_DIAGONAL_PRECOND_H