Refactor AMG

[karlrupp]

Is fairly important for users. Need to be partly rewritten for higher efficiency and for supporting multiple compute backends.

[ddemidov]

You could in principle use amgcl; it supports viennacl types out of the box: https://github.com/ddemidov/amgcl/blob/master/examples/viennacl.cpp. But it relies on Boost, and I know you are reluctant to introduce any extra dependencies into viennacl.

[karlrupp]

The main reason for refactoring is actually multi-backend capabilities, i.e. OpenMP, CUDA and OpenCL. A user can already now use amgcl right away with ViennaCL types, so I don't want to unnecessarily duplicate that code ;-)

[karlrupp]

Copying a message from Yaron Keren to viennacl-devel: I am trying to use ViennaCL with Ceemple and encountered a possible bug.

In amg_base.hpp:1045, iter2 is decremented (code below). However, nothing (at least locally in this function) prevents iter2 from being pointlist.begin() and then the result is undefined, see

http://stackoverflow.com/questions/18225651/is-stdvectorbegin-1-undefined

Also, why are both iter and iter2 required? it seems iter alone could do the job.

Thanks, Yaron

Code:

        void add_influence(amg_point* point, unsigned int add)
        {
          ListType::iterator iter = pointlist.find(point);
          // If point is not in the list then stop.
          if (iter == pointlist.end()) return;

          // Save iterator and decrement
          ListType::iterator iter2 = iter;
          iter2--;

          // Point has to be erased first as changing the value does not re-order the std::set
          pointlist.erase(iter);
          point->add_influence(add);

          // Insert point back into the list. Using the iterator improves performance. The new position has to be at the same position or to the right of the old.
          pointlist.insert(iter2,point);
        }

[yrnkrn]

It happens with the amg.cpp example but not always reproducable as pointlist is a std::set and there are no guarantees on iteration order.

If you add an assert before it will trip:

          assert(iter2 != pointlist.begin());
          iter2--;

[karlrupp]

@yrnkrn Thanks, I could reproduce the problem and pushed a fix here: https://github.com/viennacl/viennacl-dev/commit/db903fb76bf20cab86fa6baa59fa3bffe204555e

[yrnkrn]

Thanks for the quick patch. I also tried

          if (iter != pointlist.begin())
            --iter;

which produces different output than the patch code. The first discrepancy appears in Operator complexity, below. I don't know if this indicates a problem.

with --

453139 lines read. -- CG solver (CPU, no preconditioner) --

Relative residual: 0.00811647 Iterations: 300 Relative deviation from result: 0.0101884 -- CG solver (GPU, no preconditioner) -- Relative residual: 0.00542515 Iterations: 300 Relative deviation from result: 0.0103569 -- CG with AMG preconditioner, RS COARSENING, DIRECT INTERPOLATION --

• Setup phase (ublas types)...
  • * Operator complexity: 1.00004*

without --

453139 lines read. -- CG solver (CPU, no preconditioner) --

Relative residual: 0.00811647 Iterations: 300 Relative deviation from result: 0.0101884 -- CG solver (GPU, no preconditioner) -- Relative residual: 0.00542515 Iterations: 300 Relative deviation from result: 0.0103569 -- CG with AMG preconditioner, RS COARSENING, DIRECT INTERPOLATION --

• Setup phase (ublas types)...
  • * Operator complexity: 1.06404*

2014-08-12 19:51 GMT+03:00 Karl Rupp notifications@github.com:

@yrnkrn https://github.com/yrnkrn Thanks, I could reproduce the problem and pushed a fix here: db903fb https://github.com/viennacl/viennacl-dev/commit/db903fb76bf20cab86fa6baa59fa3bffe204555e

Reply to this email directly or view it on GitHub https://github.com/viennacl/viennacl-dev/issues/16#issuecomment-51942405 .

[karlrupp]

This is most likely due to undefined behavior, because STL iterators to erased elements become invalid...