Parallel ILU factorization with a reference implementation

[thoasm]

This PR adds a parallel incomplete LU factorization to Ginkgo. I did not completely follow the example of #27 because:

• now, the creation process is the same like the solvers (::build().with_iterations(X).on(exec))
• it is still usable as a Composition (see the test ApplyMethodDenseSmall, where apply is called)
• With this class, we can show that other classes depend on a LU decomposition (like a preconditioner), so they get the proper format.

For the last point it is possible for the user to change the content of L and U by calling copy_from(), so it is not 100% guaranteed that the matrices are actually triangular matrices (or that there are two of them).

closes #134, closes #27

Additionally, as part of this PR, a sort_by_column_index() method was added to the Csr matrix because the current 'ParIlu' implementation only works if the given system matrix is sorted. A reference implementation is provided which works with the custom made iterator IteratorFactory::Iterator, so no array has to be copied to work with std::sort.
It should be possible to also use this iterator to implement the OpenMP version (sort chunks in parallel, followed by a merge step), although a custom sort implementation might be faster.

TODO before merging:

• Rebase

[hartwiganzt]

Is the CI failing because it excesses the time limit?

[thoasm]

No, it fails because the test GenerateForReverseCooSmall and GenerateForReverseCsrSmall fail, which test what happens if the system_matrix is in COO format and in reverse order (stores the matrix from bottom right to upper left), or in CSR format and in reverse order (each row from right to left). All the other tests work.
I added these tests yesterday to see if it would work, and it does not (which was expected). Now, we either have to add a sort function to the CSR matrix (which would require us to implement that on all executors), or, we try to do it inside the ParIlu.

I think adding a sort function to the CSR matrix might be a good idea, including a test to see if it is sorted, or not. I expect that other solver / preconditioner / operations might also benefit from a sort function, and in CSR, it should also be possible to implement an efficient solution for the GPU.

[tcojean]

@thoasm for some reason there are quite a bit of conflicts of the CSR files (with the Hybrid conversions?) so you will have to rebase.

[thoasm]

@tcojean I know, I am not done yet with the implementation of the Csr sort function. I will do the rebase as soon as I am done there, which will be in the next couple of hours. I just did not want to do multiple rebases if multiple PRs were merged before I am done.
Currently, I am writing the tests to ensure that everything is working, then, I will push everything and remove the WIP label.

[thoasm]

This PR can finally be reviewed. All tests should succeed.

In addition to the ParIlu, a sort function was added for the Csr matrix (including a reference implementation). This was necessary because the current LU computation requires each row of the Csr matrix to be sorted by column index. There are two test cases which cover this behaviour: GenerateForReverseCsrSmall and GenerateForReverseCooSmall.

[thoasm]

Unfortunately, we get a note from g++ from the added Iterator when compiling (output see below). For clang, there is no warning. We could disable this note with #pragma GCC diagnostic and push -Wno-psabi for this function, which removes all warnings / notes for ABI changes for this function (see this stackoverflow answer for details). Alternatively, we could also just leave it as it is.

Edit: I tried the in SO suggested pragma, but it did not work. Looks like we will keep this note since I think it is a very bad idea to put -Wno-psabi (which removes this note, I tested it) in our default compiler flags since it might also remove warnings / notes for actual problems in our code.

In file included from /home/thomas/projects/ginkgo_github/reference/matrix/csr_kernels.cpp:51:
/home/thomas/projects/ginkgo_github/core/base/iterator_factory.hpp: In instantiation of ‘gko::detail::IteratorFactory<ToSortType, SecondaryType>::Reference::operator gko::detail::IteratorFactory<ToSortType, SecondaryType>::element() const [with ToSortType = int; SecondaryType = std::complex<float>]’:
/usr/include/c++/8.3.0/bits/stl_algo.h:1850:3:   required from ‘void std::__insertion_sort(_RandomAccessIterator, _RandomAccessIterator, _Compare) [with _RandomAccessIterator = gko::detail::IteratorFactory<int, std::complex<float> >::Iterator; _Compare = __gnu_cxx::__ops::_Iter_less_iter]’
/usr/include/c++/8.3.0/bits/stl_algo.h:1885:25:   required from ‘void std::__final_insertion_sort(_RandomAccessIterator, _RandomAccessIterator, _Compare) [with _RandomAccessIterator = gko::detail::IteratorFactory<int, std::complex<float> >::Iterator; _Compare = __gnu_cxx::__ops::_Iter_less_iter]’
/usr/include/c++/8.3.0/bits/stl_algo.h:1971:31:   required from ‘void std::__sort(_RandomAccessIterator, _RandomAccessIterator, _Compare) [with _RandomAccessIterator = gko::detail::IteratorFactory<int, std::complex<float> >::Iterator; _Compare = __gnu_cxx::__ops::_Iter_less_iter]’
/usr/include/c++/8.3.0/bits/stl_algo.h:4834:18:   required from ‘void std::sort(_RAIter, _RAIter) [with _RAIter = gko::detail::IteratorFactory<int, std::complex<float> >::Iterator]’
/home/thomas/projects/ginkgo_github/reference/matrix/csr_kernels.cpp:486:18:   required from ‘void gko::kernels::reference::csr::sort_by_column_index(std::shared_ptr<const gko::ReferenceExecutor>, gko::matrix::Csr<ValueType, IndexType>*) [with ValueType = std::complex<float>; IndexType = int]’
/home/thomas/projects/ginkgo_github/reference/matrix/csr_kernels.cpp:490:1:   required from here
/home/thomas/projects/ginkgo_github/core/base/iterator_factory.hpp:140:9: note: the ABI of passing structure with complex float member has changed in GCC 4.4
         operator element() const { return element{dominant(), secondary()}; }
         ^~~~~~~~

[codecov]

Codecov Report

Merging #305 into develop will decrease coverage by <.01%.
The diff coverage is 98%.

@@             Coverage Diff             @@
##           develop     #305      +/-   ##
===========================================
- Coverage    98.17%   98.17%   -0.01%     
===========================================
  Files          215      223       +8     
  Lines        16529    17180     +651     
===========================================
+ Hits         16227    16866     +639     
- Misses         302      314      +12

Impacted Files	Coverage Δ
cuda/base/device_guard.hpp	100% <ø> (ø)	⬆️
include/ginkgo/core/matrix/csr.hpp	85.48% <ø> (ø)	⬆️
omp/factorization/par_ilu_kernels.cpp	0% <0%> (ø)
omp/matrix/csr_kernels.cpp	91.66% <0%> (-2.62%)	⬇️
core/base/iterator_factory.hpp	100% <100%> (ø)
core/matrix/csr.cpp	100% <100%> (ø)	⬆️
include/ginkgo/core/factorization/par_ilu.hpp	100% <100%> (ø)
reference/matrix/csr_kernels.cpp	100% <100%> (ø)	⬆️
core/factorization/par_ilu.cpp	100% <100%> (ø)
core/test/base/iterator_factory.cpp	100% <100%> (ø)
... and 14 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2f2f09e...dbabfdb. Read the comment docs.

[thoasm]

@tcojean I am currently failing SonarCloud because I have Duplicated Lines on New Code 5.3% (see here for details). This duplication is only in {omp,cuda}/factorization/par_ilu_kernels.cpp, which are the stubs for the new kernels, which are not implemented.
Is there a way to ignore or get around that?

Also, the code coverage is still not perfect apparently, but this is due to the exact same reason: I am not testing the stubs. Additionally, it detects the two apply_impl in my DummyLinOp as also not tested, even though they are inside the test file. They are not required for the test, but are required by the class, so it is not an abstract class.
I could probably use an already existing class that is not convertible to Csr (the important property for the test), but that would create a dependency that I would like to avoid (if we use a matrix, we would need to change the test as soon as we actually support a conversion to Csr).

[tcojean]

LGTM. Some minor documentation and include issues and other improvement suggestions.

I think you should also add some files for the documentation, namely a file in doc/headers which defines a group factor (see linop.hpp as example), and in addition, an update to the modules.dot at the same time (in fact is this used @pratikvn, I can't remember?).

[pratikvn]

Mostly minor and nits. One possibly major: Is the factorization of big_mtx in the reference test correct ? I am getting a different result in MATLAB, but it is possible that I entered the matrix incorrectly.

[thoasm]

All comments should now be dealt with. I also added doc/headers/factor.hpp and inserted factor to doc/headers/modules.dot (a very short description, but I am not sure what else I should have wrote there).

[pratikvn]

LGTM!

[pratikvn]

Some minor things

[tcojean]

For the CI system, we can add to .before_script_template in .gitlab-ci.yml the following:

  - export CUDA_VISIBLE_DEVICES=0

This should allow only the usage of the first GPU for all job. But we have a test such as sending data from a CUDA executor to another, in which case that will happen on the same GPU instead of on the two distinct GPUs.