amg_coarse.hpp

Go to the documentation of this file.

#ifndef VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP
#define VIENNACL_LINALG_DETAIL_AMG_AMG_COARSE_HPP

/* =========================================================================
   Copyright (c) 2010-2014, 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 <cmath>
#include "viennacl/linalg/detail/amg/amg_base.hpp"

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

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

namespace viennacl
{
  namespace linalg
  {
    namespace detail
    {
      namespace amg
      {

    template <typename InternalType1, typename InternalType2, typename InternalType3>
    void amg_coarse(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, InternalType3 & Slicing, amg_tag & tag)
    {
      switch (tag.get_coarse())
      {
        case VIENNACL_AMG_COARSE_RS: amg_coarse_classic (level, A, Pointvector, tag); break;
        case VIENNACL_AMG_COARSE_ONEPASS: amg_coarse_classic_onepass (level, A, Pointvector, tag); break;
        case VIENNACL_AMG_COARSE_RS0: amg_coarse_rs0 (level, A, Pointvector, Slicing, tag); break;
        case VIENNACL_AMG_COARSE_RS3: amg_coarse_rs3 (level, A, Pointvector, Slicing, tag); break;
        case VIENNACL_AMG_COARSE_AG:   amg_coarse_ag (level, A, Pointvector, tag); break;
      }
    }

    template <typename InternalType1, typename InternalType2>
    void amg_influence(unsigned int level, InternalType1 const & A, 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::value_type ScalarType;
      typedef typename SparseMatrixType::const_iterator1 ConstRowIterator;
      typedef typename SparseMatrixType::const_iterator2 ConstColIterator;

      ScalarType max;
      int diag_sign;
      //unsigned int i;

#ifdef VIENNACL_WITH_OPENMP
      #pragma omp parallel for private (max,diag_sign)
#endif
      for (long i=0; i<static_cast<long>(A[level].size1()); ++i)
      {
        diag_sign = 1;
        if (A[level](i,i) < 0)
          diag_sign = -1;

        ConstRowIterator row_iter = A[level].begin1();
        row_iter += i;
        // Find greatest non-diagonal negative value (positive if diagonal is negative) in row
        max = 0;
        for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
        {
            if (i == (unsigned int) col_iter.index2()) continue;
            if (diag_sign == 1)
              if (max > *col_iter)  max = *col_iter;
            if (diag_sign == -1)
              if (max < *col_iter)  max = *col_iter;
        }

        // If maximum is 0 then the row is independent of the others
        if (max == 0)
          continue;

        // Find all points that strongly influence current point (Yang, p.5)
        for (ConstColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
        {
          unsigned int j = static_cast<unsigned int>(col_iter.index2());
          if (i == j) continue;
          if (diag_sign * (-*col_iter) >= tag.get_threshold() * (diag_sign * (-max)))
          {
            // Strong influence from j to i found, save information
            Pointvector[level][i]->add_influencing_point(Pointvector[level][j]);
          }
        }
      }

      #ifdef VIENNACL_AMG_DEBUG
      std::cout << "Influence Matrix: " << std::endl;
      boost::numeric::ublas::matrix<bool> mat;
      Pointvector[level].get_influence_matrix(mat);
      printmatrix (mat);
      #endif

      // Save influenced points
      for (typename PointVectorType::iterator iter = Pointvector[level].begin(); iter != Pointvector[level].end(); ++iter)
      {
        for (typename amg_point::iterator iter2 = (*iter)->begin_influencing(); iter2 != (*iter)->end_influencing(); ++iter2)
        {
          (*iter2)->add_influenced_point(*iter);
        }
      }

      #ifdef VIENNACL_AMG_DEBUG
      std::cout << "Influence Measures: " << std::endl;
      boost::numeric::ublas::vector<unsigned int> temp;
      Pointvector[level].get_influence(temp);
      printvector (temp);
      std::cout << "Point Sorting: " << std::endl;
      Pointvector[level].get_sorting(temp);
      printvector (temp);
      #endif
    }

    template <typename InternalType1, typename InternalType2>
    void amg_coarse_classic_onepass(unsigned int level, InternalType1 & A, InternalType2 & Pointvector, amg_tag & tag)
    {
      amg_point* c_point, *point1, *point2;

      // Check and save all strong influences
      amg_influence (level, A, Pointvector, tag);

      // Traverse through points and calculate initial influence measure
      long i;
#ifdef VIENNACL_WITH_OPENMP
      #pragma omp parallel for private (i)
#endif
      for (i=0; i<static_cast<long>(Pointvector[level].size()); ++i)
  Pointvector[level][i]->calc_influence();

       // Do initial sorting
      Pointvector[level].sort();

      // Get undecided point with highest influence measure
      while ((c_point = Pointvector[level].get_nextpoint()) != NULL)
      {
        // Make this point C point
        Pointvector[level].make_cpoint(c_point);

        // All strongly influenced points become F points