viennacl/linalg/amg.hpp

Go to the documentation of this file.

#ifndef VIENNACL_LINALG_AMG_HPP_
#define VIENNACL_LINALG_AMG_HPP_

/* =========================================================================
   Copyright (c) 2010-2013, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.
   Portions of this software are copyright by UChicago Argonne, LLC.

                            -----------------
                  ViennaCL - The Vienna Computing Library
                            -----------------

   Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
               
   (A list of authors and contributors can be found in the PDF manual)

   License:         MIT (X11), see file LICENSE in the base directory
============================================================================= */

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/triangular.hpp> 
#include <vector>
#include <cmath>
#include "viennacl/forwards.h"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/linalg/direct_solve.hpp"

#include "viennacl/linalg/detail/amg/amg_base.hpp"
#include "viennacl/linalg/detail/amg/amg_coarse.hpp"
#include "viennacl/linalg/detail/amg/amg_interpol.hpp"

#include <map>

#ifdef VIENNACL_WITH_OPENMP
 #include <omp.h>
#endif

#include "viennacl/linalg/detail/amg/amg_debug.hpp"

#define VIENNACL_AMG_COARSE_LIMIT 50
#define VIENNACL_AMG_MAX_LEVELS 100

namespace viennacl
{
  namespace linalg
  {    
    typedef detail::amg::amg_tag          amg_tag;
    
    
    
    template <typename InternalType1, typename InternalType2>
    void amg_setup(InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
    {
      typedef typename InternalType1::value_type SparseMatrixType;
      typedef typename InternalType2::value_type PointVectorType;     
      typedef typename SparseMatrixType::value_type ScalarType;   
      typedef typename SparseMatrixType::iterator1 InternalRowIterator;
      typedef typename SparseMatrixType::iterator2 InternalColIterator;
      
      unsigned int i, iterations, c_points, f_points;
      detail::amg::amg_slicing<InternalType1,InternalType2> Slicing;
      
      // Set number of iterations. If automatic coarse grid construction is chosen (0), then set a maximum size and stop during the process.
      iterations = tag.get_coarselevels();
      if (iterations == 0)
        iterations = VIENNACL_AMG_MAX_LEVELS;
             
      // For parallel coarsenings build data structures (number of threads set automatically).
      if (tag.get_coarse() == VIENNACL_AMG_COARSE_RS0 || tag.get_coarse() == VIENNACL_AMG_COARSE_RS3)
        Slicing.init(iterations);
      
      for (i=0; i<iterations; ++i)
      {  
        // Initialize Pointvector on level i and construct points.
        Pointvector[i] = PointVectorType(A[i].size1());
        Pointvector[i].init_points();
        
        // Construct C and F points on coarse level (i is fine level, i+1 coarse level).
        detail::amg::amg_coarse (i, A, Pointvector, Slicing, tag);

        // Calculate number of C and F points on level i.
        c_points = Pointvector[i].get_cpoints();
        f_points = Pointvector[i].get_fpoints();

        #if defined (VIENNACL_AMG_DEBUG) //or defined(VIENNACL_AMG_DEBUGBENCH) 
        std::cout << "Level " << i << ": ";
        std::cout << "No of C points = " << c_points << ", ";
        std::cout << "No of F points = " << f_points << std::endl;
        #endif
        
        // Stop routine when the maximal coarse level is found (no C or F point). Coarsest level is level i.
        if (c_points == 0 || f_points == 0)
          break;
          
        // Construct interpolation matrix for level i.
        detail::amg::amg_interpol (i, A, P, Pointvector, tag);
        
        // Compute coarse grid operator (A[i+1] = R * A[i] * P) with R = trans(P).
        detail::amg::amg_galerkin_prod(A[i], P[i], A[i+1]);
        
        // Test triple matrix product. Very slow for large matrix sizes (ublas).
        // test_triplematprod(A[i],P[i],A[i+1]);
        
        Pointvector[i].delete_points();
        
        #ifdef VIENNACL_AMG_DEBUG
        std::cout << "Coarse Grid Operator Matrix:" << std::endl;
        printmatrix (A[i+1]);
        #endif  
        
        // If Limit of coarse points is reached then stop. Coarsest level is level i+1.
        if (tag.get_coarselevels() == 0 && c_points <= VIENNACL_AMG_COARSE_LIMIT)
        {
          tag.set_coarselevels(i+1);
          return;
        }
      }
      tag.set_coarselevels(i);  
    }
    
    template <typename MatrixType, typename InternalType1, typename InternalType2>
    void amg_init(MatrixType const & mat, InternalType1 & A, InternalType1 & P, InternalType2 & Pointvector, amg_tag & tag)
    {
      typedef typename MatrixType::value_type ScalarType;
      typedef typename InternalType1::value_type SparseMatrixType;
      
      if (tag.get_coarselevels() > 0)
      {
        A.resize(tag.get_coarselevels()+1);
        P.resize(tag.get_coarselevels());
        Pointvector.resize(tag.get_coarselevels());
      }
      else
      {
        A.resize(VIENNACL_AMG_MAX_LEVELS+1);
        P.resize(VIENNACL_AMG_MAX_LEVELS);
        Pointvector.resize(VIENNACL_AMG_MAX_LEVELS);
      }
      
      // Insert operator matrix as operator for finest level.
      SparseMatrixType A0 (mat);
      A.insert_element (0, A0);  
    }
    
    template <typename InternalType1, typename InternalType2>
    void amg_transform_cpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup, amg_tag & tag)
    {  
      typedef typename InternalType1::value_type MatrixType;
      
      // Resize internal data structures to actual size.
      A.resize(tag.get_coarselevels()+1);
      P.resize(tag.get_coarselevels());
      R.resize(tag.get_coarselevels());
      
      // Transform into matrix type.
      for (unsigned int i=0; i<tag.get_coarselevels()+1; ++i)
      {
        A[i].resize(A_setup[i].size1(),A_setup[i].size2(),false);
        A[i] = A_setup[i];
      }
      for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
      {
        P[i].resize(P_setup[i].size1(),P_setup[i].size2(),false);
        P[i] = P_setup[i];
      }
      for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
      {
        R[i].resize(P_setup[i].size2(),P_setup[i].size1(),false);
        P_setup[i].set_trans(true);
        R[i] = P_setup[i];
        P_setup[i].set_trans(false);
      }
    }
    
    template <typename InternalType1, typename InternalType2>
    void amg_transform_gpu (InternalType1 & A, InternalType1 & P, InternalType1 & R, InternalType2 & A_setup, InternalType2 & P_setup, amg_tag & tag)
    {
      typedef typename InternalType1::value_type MatrixType;
      typedef typename InternalType2::value_type::value_type ScalarType;
      
      // Resize internal data structures to actual size.
      A.resize(tag.get_coarselevels()+1);
      P.resize(tag.get_coarselevels());
      R.resize(tag.get_coarselevels());
      
      // Copy to GPU using the internal sparse matrix structure: std::vector<std::map>.
      for (unsigned int i=0; i<tag.get_coarselevels()+1; ++i)
      {
        //A[i].resize(A_setup[i].size1(),A_setup[i].size2(),false);
        viennacl::copy(*(A_setup[i].get_internal_pointer()),A[i]);
      }
      for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
      {
        //P[i].resize(P_setup[i].size1(),P_setup[i].size2(),false);
        viennacl::copy(*(P_setup[i].get_internal_pointer()),P[i]);
        //viennacl::copy((boost::numeric::ublas::compressed_matrix<ScalarType>)P_setup[i],P[i]);
      }
      for (unsigned int i=0; i<tag.get_coarselevels(); ++i)
      {
        //R[i].resize(P_setup[i].size2(),P_setup[i].size1(),false);
        P_setup[i].set_trans(true);
        viennacl::copy(*(P_setup[i].get_internal_pointer()),R[i]);
        P_setup[i].set_trans(false);
      }
    }
    
    template <typename InternalVectorType, typename SparseMatrixType>
    void amg_setup_apply (InternalVectorType & result, InternalVectorType & rhs, InternalVectorType & residual, SparseMatrixType const & A, amg_tag const & tag)
    {        
      typedef typename InternalVectorType::value_type VectorType;
      
      result.resize(tag.get_coarselevels()+1);
      rhs.resize(tag.get_coarselevels()+1);
      residual.resize(tag.get_coarselevels());
            
      for (unsigned int level=0; level < tag.get_coarselevels()+1; ++level)
      {
        result[level] = VectorType(A[level].size1());
        result[level].clear();
        rhs[level] = VectorType(A[level].size1());
        rhs[level].clear();
      }
      for (unsigned int level=0; level < tag.get_coarselevels(); ++level)
      {
        residual[level] = VectorType(A[level].size1());
        residual[level].clear();
      }
    }
    template <typename ScalarType, typename SparseMatrixType>
    void amg_lu(boost::numeric::ublas::compressed_matrix<ScalarType> & op, boost::numeric::ublas::permutation_matrix<> & Permutation, SparseMatrixType const & A)
    {
      typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
      typedef typename SparseMatrixType::const_iterator2 ConstColIterator;
    
      // Copy to operator matrix. Needed 
      op.resize(A.size1(),A.size2(),false);
      for (ConstRowIterator row_iter = A.begin1(); row_iter != A.end1(); ++row_iter)
        for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
          op (col_iter.index1(), col_iter.index2()) = *col_iter;
      
      // Permutation matrix has to be reinitialized with actual size. Do not clear() or resize()!
      Permutation = boost::numeric::ublas::permutation_matrix<> (op.size1());
      boost::numeric::ublas::lu_factorize(op,Permutation);
    }
    
    template <typename MatrixType>
    class amg_precond
    {
      typedef typename MatrixType::value_type ScalarType;
      typedef boost::numeric::ublas::vector<ScalarType> VectorType;
      typedef detail::amg::amg_sparsematrix<ScalarType> SparseMatrixType;
      typedef detail::amg::amg_pointvector PointVectorType;
      
      typedef typename SparseMatrixType::const_iterator1 InternalConstRowIterator;
      typedef typename SparseMatrixType::const_iterator2 InternalConstColIterator;
      typedef typename SparseMatrixType::iterator1 InternalRowIterator;
      typedef typename SparseMatrixType::iterator2 InternalColIterator;
      
      boost::numeric::ublas::vector <SparseMatrixType> A_setup;
      boost::numeric::ublas::vector <SparseMatrixType> P_setup;
      boost::numeric::ublas::vector <MatrixType> A;
      boost::numeric::ublas::vector <MatrixType> P;
      boost::numeric::ublas::vector <MatrixType> R;
      boost::numeric::ublas::vector <PointVectorType> Pointvector;
       
      mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
      mutable boost::numeric::ublas::permutation_matrix<> Permutation;  
      
      mutable boost::numeric::ublas::vector <VectorType> result;
      mutable boost::numeric::ublas::vector <VectorType> rhs;
      mutable boost::numeric::ublas::vector <VectorType> residual;
      
      mutable bool done_init_apply;
          
      amg_tag tag_;
    public:
    
      amg_precond(): Permutation(0) {}
      amg_precond(MatrixType const & mat, amg_tag const & tag): Permutation(0)
      {
        tag_ = tag;
        // Initialize data structures.
        amg_init (mat,A_setup,P_setup,Pointvector,tag_);
        
        done_init_apply = false;
      }
      
      void setup()
      {
        // Start setup phase.
        amg_setup(A_setup,P_setup,Pointvector,tag_);
        // Transform to CPU-Matrixtype for precondition phase.
        amg_transform_cpu(A,P,R,A_setup,P_setup,tag_);
        
        done_init_apply = false;
      }
      
      void init_apply() const
      {
        // Setup precondition phase (Data structures).
        amg_setup_apply(result,rhs,residual,A_setup,tag_);
        // Do LU factorization for direct solve.
        amg_lu(op,Permutation,A_setup[tag_.get_coarselevels()]);
        
        done_init_apply = true;
      }
      
      template <typename VectorType>
      ScalarType calc_complexity(VectorType & avgstencil)
      {
        avgstencil = VectorType (tag_.get_coarselevels()+1);
        unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
          
        for (unsigned int level=0; level < tag_.get_coarselevels()+1; ++level)
        {
          level_coefficients = 0;
          for (InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
          { 
            for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
            {
              if (level == 0)
                systemmat_nonzero++;
              nonzero++;
              level_coefficients++;
            }
          }
          avgstencil[level] = level_coefficients/static_cast<ScalarType>(A_setup[level].size1());
        }
        return nonzero/static_cast<ScalarType>(systemmat_nonzero);
      }

      template <typename VectorType>
      void apply(VectorType & vec) const
      {   
        // Build data structures and do lu factorization before first iteration step.
        if (!done_init_apply)
          init_apply();
        
        int level;
        
        // Precondition operation (Yang, p.3)
        rhs[0] = vec;
        for (level=0; level <static_cast<int>(tag_.get_coarselevels()); level++)
        {    
          result[level].clear();
          
          // Apply Smoother presmooth_ times.
          smooth_jacobi (level, tag_.get_presmooth(), result[level], rhs[level]);    
          
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "After presmooth:" << std::endl;
          printvector(result[level]);
          #endif

          // Compute residual.
          residual[level] = rhs[level] - boost::numeric::ublas::prod (A[level],result[level]);
          
          #ifdef VIENNACL_AMG_DEBUG          
          std::cout << "Residual:" << std::endl;
          printvector(residual[level]);
          #endif
          
          // Restrict to coarse level. Restricted residual is RHS of coarse level.
          rhs[level+1] = boost::numeric::ublas::prod (R[level],residual[level]);
          
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Restricted Residual: " << std::endl;
          printvector(rhs[level+1]);
          #endif
        }
          
        // On highest level use direct solve to solve equation.
        result[level] = rhs[level];
        boost::numeric::ublas::lu_substitute(op,Permutation,result[level]);

        #ifdef VIENNACL_AMG_DEBUG
        std::cout << "After direct solve: " << std::endl;
        printvector (result[level]);
        #endif
          
        for (level=tag_.get_coarselevels()-1; level >= 0; level--)
        {       
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Coarse Error: " << std::endl;
          printvector(result[level+1]);
          #endif
          
          // Interpolate error to fine level. Correct solution by adding error.
          result[level] += boost::numeric::ublas::prod (P[level], result[level+1]);
              
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Corrected Result: " << std::endl;
          printvector (result[level]);
          #endif
          
          // Apply Smoother postsmooth_ times.
          smooth_jacobi (level, tag_.get_postsmooth(), result[level], rhs[level]);
          
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "After postsmooth: " << std::endl;
          printvector (result[level]);
          #endif
        }    
        vec = result[0];
      }
      
      template <typename VectorType>
      void smooth_jacobi(int level, int const iterations, VectorType & x, VectorType const & rhs) const
      {
        VectorType old_result (x.size());
        unsigned int index;
        ScalarType sum = 0, diag = 1;
        
        for (int i=0; i<iterations; ++i)
        {
          old_result = x;
          x.clear();
#ifdef VIENNACL_WITH_OPENMP
          #pragma omp parallel for private (sum,diag) shared (rhs,x)
#endif          
          for (index=0; index<A_setup[level].size1(); ++index)  
          {
            InternalConstRowIterator row_iter = A_setup[level].begin1();
            row_iter += index; 
            sum = 0;
            diag = 1;
            for (InternalConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
            {
              if (col_iter.index1() == col_iter.index2())
                diag = *col_iter;
              else
                sum += *col_iter * old_result[col_iter.index2()];
            }
            x[index]= static_cast<ScalarType>(tag_.get_jacobiweight()) * (rhs[index] - sum) / diag + (1-static_cast<ScalarType>(tag_.get_jacobiweight())) * old_result[index];
          }
        }
      }
      
      amg_tag & tag() { return tag_; }
    };
    
    template <typename ScalarType, unsigned int MAT_ALIGNMENT>
    class amg_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
    {
      typedef viennacl::compressed_matrix<ScalarType, MAT_ALIGNMENT> MatrixType;
      typedef viennacl::vector<ScalarType> VectorType;
      typedef detail::amg::amg_sparsematrix<ScalarType> SparseMatrixType;
      typedef detail::amg::amg_pointvector PointVectorType;
      
      typedef typename SparseMatrixType::const_iterator1 InternalConstRowIterator;
      typedef typename SparseMatrixType::const_iterator2 InternalConstColIterator;
      typedef typename SparseMatrixType::iterator1 InternalRowIterator;
      typedef typename SparseMatrixType::iterator2 InternalColIterator;
      
      boost::numeric::ublas::vector <SparseMatrixType> A_setup;
      boost::numeric::ublas::vector <SparseMatrixType> P_setup;
      boost::numeric::ublas::vector <MatrixType> A;
      boost::numeric::ublas::vector <MatrixType> P;
      boost::numeric::ublas::vector <MatrixType> R;
      boost::numeric::ublas::vector <PointVectorType> Pointvector;
      
      mutable boost::numeric::ublas::compressed_matrix<ScalarType> op;
      mutable boost::numeric::ublas::permutation_matrix<> Permutation;  
      
      mutable boost::numeric::ublas::vector <VectorType> result;
      mutable boost::numeric::ublas::vector <VectorType> rhs;
      mutable boost::numeric::ublas::vector <VectorType> residual;
          
      mutable bool done_init_apply;
    
      amg_tag tag_;
      
    public:
      
      amg_precond(): Permutation(0) {}
      
      amg_precond(compressed_matrix<ScalarType, MAT_ALIGNMENT> const & mat, amg_tag const & tag): Permutation(0)
      {
        tag_ = tag;
        
        // Copy to CPU. Internal structure of sparse matrix is used for copy operation.
        std::vector<std::map<unsigned int, ScalarType> > mat2 = std::vector<std::map<unsigned int, ScalarType> >(mat.size1());
        viennacl::copy(mat, mat2);
        
        // Initialize data structures.
        amg_init (mat2,A_setup,P_setup,Pointvector,tag_);
          
        done_init_apply = false;
      }
      
      void setup()
      {
        // Start setup phase.
        amg_setup(A_setup,P_setup,Pointvector, tag_);  
        // Transform to GPU-Matrixtype for precondition phase.
        amg_transform_gpu(A,P,R,A_setup,P_setup, tag_);
        
        done_init_apply = false;
      }
      
      void init_apply() const
      {
        // Setup precondition phase (Data structures).
        amg_setup_apply(result,rhs,residual,A_setup,tag_);
        // Do LU factorization for direct solve.
        amg_lu(op,Permutation,A_setup[tag_.get_coarselevels()]);
        
        done_init_apply = true;
      }

      template <typename VectorType>
      ScalarType calc_complexity(VectorType & avgstencil)
      {
        avgstencil = VectorType (tag_.get_coarselevels()+1);
        unsigned int nonzero=0, systemmat_nonzero=0, level_coefficients=0;
        
        for (unsigned int level=0; level < tag_.get_coarselevels()+1; ++level)
        {
          level_coefficients = 0;
          for (InternalRowIterator row_iter = A_setup[level].begin1(); row_iter != A_setup[level].end1(); ++row_iter)
          { 
            for (InternalColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
            {
              if (level == 0)
                systemmat_nonzero++;
              nonzero++;
              level_coefficients++;
            }
          }
          avgstencil[level] = level_coefficients/(double)A[level].size1();
        }
        return nonzero/static_cast<double>(systemmat_nonzero);
      }

      template <typename VectorType>
      void apply(VectorType & vec) const
      {
        if (!done_init_apply)
          init_apply();  
        
        int level;
        
        // Precondition operation (Yang, p.3).
        rhs[0] = vec;
        for (level=0; level <static_cast<int>(tag_.get_coarselevels()); level++)
        {    
          result[level].clear();
          
          // Apply Smoother presmooth_ times.
          smooth_jacobi (level, tag_.get_presmooth(), result[level], rhs[level]);

          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "After presmooth: " << std::endl;
          printvector(result[level]);
          #endif
          
          // Compute residual.
          residual[level] = rhs[level] - viennacl::linalg::prod (A[level],result[level]);
          
          #ifdef VIENNACL_AMG_DEBUG             
          std::cout << "Residual: " << std::endl;
          printvector(residual[level]);
          #endif
          
          // Restrict to coarse level. Result is RHS of coarse level equation.
          //residual_coarse[level] = viennacl::linalg::prod(R[level],residual[level]);
          rhs[level+1] = viennacl::linalg::prod(R[level],residual[level]);
          
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Restricted Residual: " << std::endl;
          printvector(rhs[level+1]);
          #endif
        }
          
        // On highest level use direct solve to solve equation (on the CPU)
        //TODO: Use GPU direct solve!
        result[level] = rhs[level];
        boost::numeric::ublas::vector <ScalarType> result_cpu (result[level].size());

        copy (result[level],result_cpu);
        boost::numeric::ublas::lu_substitute(op,Permutation,result_cpu);
        copy (result_cpu, result[level]);
        
        #ifdef VIENNACL_AMG_DEBUG
        std::cout << "After direct solve: " << std::endl;
        printvector (result[level]);
        #endif
          
        for (level=tag_.get_coarselevels()-1; level >= 0; level--)
        {   
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Coarse Error: " << std::endl;
          printvector(result[level+1]);
          #endif
          
          // Interpolate error to fine level and correct solution.
          result[level] += viennacl::linalg::prod(P[level],result[level+1]);
              
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "Corrected Result: " << std::endl;
          printvector (result[level]);
          #endif
          
          // Apply Smoother postsmooth_ times.
          smooth_jacobi (level, tag_.get_postsmooth(), result[level], rhs[level]);
          
          #ifdef VIENNACL_AMG_DEBUG
          std::cout << "After postsmooth: " << std::endl;
          printvector (result[level]);
          #endif
        }    
        vec = result[0];
      }
      
      template <typename VectorType>
      void smooth_jacobi(int level, unsigned int iterations, VectorType & x, VectorType const & rhs) const
      {     
        VectorType old_result (x.size());
  
        viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::init();
        viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
                    "jacobi");
  
        for (unsigned int i=0; i<iterations; ++i)
        {
          old_result = x;    
          x.clear();
          viennacl::ocl::enqueue(k(A[level].handle1().opencl_handle(), A[level].handle2().opencl_handle(), A[level].handle().opencl_handle(),
                                  static_cast<ScalarType>(tag_.get_jacobiweight()), 
                                  viennacl::traits::opencl_handle(old_result),
                                  viennacl::traits::opencl_handle(x),
                                  viennacl::traits::opencl_handle(rhs),
                                  static_cast<cl_uint>(rhs.size()))); 
          
        }
      }
      
      amg_tag & tag() { return tag_; }
    };

  }
}



#endif