Skip to content

Commit

Permalink
Added option to sort Csr Matrices
Browse files Browse the repository at this point in the history 
Sorting of Csr matrices is required for ParIlu to work with any
kind of matrix input.
Currently, only the reference implementation is done. It is performed
with `std::sort` with a custom Iterator. That helps to prevent both
unnecessary copies and the need to implement a custom sorting algorithm.
  • Loading branch information
Thomas Grützmacher committed May 23, 2019
1 parent 7339437 commit 2b27fc7
Show file tree
Hide file tree
Showing 11 changed files with 459 additions and 34 deletions.
12 changes: 12 additions & 0 deletions core/device_hooks/common_kernels.inc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,6 +360,18 @@ GKO_NOT_COMPILED(GKO_HOOK_MODULE);
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_CALCULATE_MAX_NNZ_PER_ROW_KERNEL);

template <typename ValueType, typename IndexType>
GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX(ValueType, IndexType)
GKO_NOT_COMPILED(GKO_HOOK_MODULE);
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX);

template <typename ValueType, typename IndexType>
GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX(ValueType, IndexType)
GKO_NOT_COMPILED(GKO_HOOK_MODULE);
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX);


} // namespace csr

Expand Down
6 changes: 5 additions & 1 deletion core/factorization/par_ilu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,13 @@ GKO_REGISTER_OPERATION(csr_transpose, csr::transpose);
template <typename ValueType, typename IndexType>
std::unique_ptr<Composition<ValueType>>
ParIlu<ValueType, IndexType>::generate_l_u(
const std::shared_ptr<const LinOp> &system_matrix) const
const std::shared_ptr<const LinOp> &system_matrix, bool skip_sorting) const
{
using CsrMatrix = matrix::Csr<ValueType, IndexType>;
using CooMatrix = matrix::Coo<ValueType, IndexType>;

// TODO include the sorting!

const auto exec = this->get_executor();
// Only copies the matrix if it is not on the same executor or was not in
// the right format. Throws an exception if it is not convertable.
Expand Down Expand Up @@ -105,6 +107,8 @@ ParIlu<ValueType, IndexType>::generate_l_u(

// At first, test if the given system_matrix was already a Coo matrix,
// so no conversion would be necessary.
// TODO if the matrix is unsorted, it maybe needs be converted from the
// sorted CSR matrix (TEST it)
std::unique_ptr<CooMatrix> coo_system_matrix_unique_ptr{nullptr};
auto coo_system_matrix_ptr =
dynamic_cast<const CooMatrix *>(system_matrix.get());
Expand Down
25 changes: 23 additions & 2 deletions core/matrix/csr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,10 +58,13 @@ GKO_REGISTER_OPERATION(convert_to_dense, csr::convert_to_dense);
GKO_REGISTER_OPERATION(convert_to_sellp, csr::convert_to_sellp);
GKO_REGISTER_OPERATION(calculate_total_cols, csr::calculate_total_cols);
GKO_REGISTER_OPERATION(convert_to_ell, csr::convert_to_ell);
GKO_REGISTER_OPERATION(transpose, csr::transpose);
GKO_REGISTER_OPERATION(conj_transpose, csr::conj_transpose);
GKO_REGISTER_OPERATION(calculate_max_nnz_per_row,
csr::calculate_max_nnz_per_row);
GKO_REGISTER_OPERATION(transpose, csr::transpose);
GKO_REGISTER_OPERATION(conj_transpose, csr::conj_transpose);
GKO_REGISTER_OPERATION(sort_by_column_index, csr::sort_by_column_index);
GKO_REGISTER_OPERATION(is_sorted_by_column_index,
csr::is_sorted_by_column_index);


} // namespace csr
Expand Down Expand Up @@ -257,6 +260,24 @@ std::unique_ptr<LinOp> Csr<ValueType, IndexType>::conj_transpose() const
}


template <typename ValueType, typename IndexType>
void Csr<ValueType, IndexType>::sort_by_column_idx()
{
auto exec = this->get_executor();
exec->run(csr::make_sort_by_column_index(this));
}


template <typename ValueType, typename IndexType>
bool Csr<ValueType, IndexType>::is_sorted_by_column_idx() const
{
auto exec = this->get_executor();
bool is_sorted;
exec->run(csr::make_is_sorted_by_column_index(this, &is_sorted));
return is_sorted;
}


#define GKO_DECLARE_CSR_MATRIX(ValueType, IndexType) \
class Csr<ValueType, IndexType>
GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_CSR_MATRIX);
Expand Down
57 changes: 35 additions & 22 deletions core/matrix/csr_kernels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <ginkgo/core/matrix/coo.hpp>
#include <ginkgo/core/matrix/csr.hpp>
#include <ginkgo/core/matrix/dense.hpp>
#include <ginkgo/core/matrix/sellp.hpp>
#include <ginkgo/core/matrix/ell.hpp>
#include <ginkgo/core/matrix/sellp.hpp>


namespace gko {
Expand Down Expand Up @@ -100,27 +100,40 @@ namespace kernels {
std::shared_ptr<const DefaultExecutor> exec, \
const matrix::Csr<ValueType, IndexType> *source, size_type *result)

#define GKO_DECLARE_ALL_AS_TEMPLATES \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_SPMV_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_ADVANCED_SPMV_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_DENSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_COO_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_SELLP_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_ELL_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CALCULATE_TOTAL_COLS_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_TRANSPOSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONJ_TRANSPOSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CALCULATE_MAX_NNZ_PER_ROW_KERNEL(ValueType, IndexType)
#define GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX(ValueType, IndexType) \
void sort_by_column_index(std::shared_ptr<const DefaultExecutor> exec, \
matrix::Csr<ValueType, IndexType> *to_sort)

#define GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX(ValueType, IndexType) \
void is_sorted_by_column_index( \
std::shared_ptr<const DefaultExecutor> exec, \
const matrix::Csr<ValueType, IndexType> *to_check, bool *is_sorted)

#define GKO_DECLARE_ALL_AS_TEMPLATES \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_SPMV_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_ADVANCED_SPMV_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_DENSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_COO_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_SELLP_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONVERT_TO_ELL_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CALCULATE_TOTAL_COLS_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_TRANSPOSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CONJ_TRANSPOSE_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_CALCULATE_MAX_NNZ_PER_ROW_KERNEL(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX(ValueType, IndexType); \
template <typename ValueType, typename IndexType> \
GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX(ValueType, IndexType)


namespace omp {
Expand Down
23 changes: 23 additions & 0 deletions cuda/matrix/csr_kernels.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1418,6 +1418,29 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_CALCULATE_MAX_NNZ_PER_ROW_KERNEL);


template <typename ValueType, typename IndexType>
void sort_by_column_index(std::shared_ptr<const CudaExecutor> exec,
matrix::Csr<ValueType, IndexType> *to_sort)
{
GKO_NOT_IMPLEMENTED;
}

GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX);


template <typename ValueType, typename IndexType>
void is_sorted_by_column_index(
std::shared_ptr<const CudaExecutor> exec,
const matrix::Csr<ValueType, IndexType> *to_check, bool *is_sorted)
{
GKO_NOT_IMPLEMENTED;
}

GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX);


} // namespace csr
} // namespace cuda
} // namespace kernels
Expand Down
29 changes: 26 additions & 3 deletions include/ginkgo/core/factorization/par_ilu.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,25 @@ class ParIlu : public Composition<ValueType> {
* the factorization.
*/
unsigned int GKO_FACTORY_PARAMETER(iterations, 1);

/**
* @brief `true` means it is known that the matrix given to this
* factory will be sorted first by row, then by column index,
* `false` means it is unknown or not sorted, so an additional
* sorting step will be performed during the factorization
* (it will not change the matrix given).
* The matrix must be sorted for this factorization to work.
*
* The `system_matrix`, which will be given to this factory, must be
* sorted (first by row, then by column) in order for the algorithm
* to work. If it is known that the matrix will be sorted, this
* parameter can be set to `true` to skip the sorting (therefore,
* shortening the runtime).
* However, if it is unknown or if the matrix is known to be not sorted,
* it must remain `false`, otherwise, the factorization might be
* incorrect.
*/
bool GKO_FACTORY_PARAMETER(skip_sorting, false);
};
GKO_ENABLE_LIN_OP_FACTORY(ParIlu, parameters, Factory);
GKO_ENABLE_BUILD_METHOD(Factory);
Expand All @@ -114,7 +133,7 @@ class ParIlu : public Composition<ValueType> {
: Composition<ValueType>(factory->get_executor()),
parameters_{factory->get_parameters()}
{
generate_l_u(system_matrix)->move_to(this);
generate_l_u(system_matrix, parameters_.skip_sorting)->move_to(this);
}

/**
Expand All @@ -125,12 +144,16 @@ class ParIlu : public Composition<ValueType> {
*
* @param system_matrix the source matrix used to generate the factors.
* @note: system_matrix must be convertable to a Csr
Matrix, otherwise, an exception is thrown.
* Matrix, otherwise, an exception is thrown.
* @param skip_sorting if set to `true`, the sorting will be skipped.
* @note: If the matrix was not sorted, the
* factorization might be wrong.
* @return A Composition, containing the incomplete LU factors for the
* given system_matrix (first element is L, then U)
*/
std::unique_ptr<Composition<ValueType>> generate_l_u(
const std::shared_ptr<const LinOp> &system_matrix) const;
const std::shared_ptr<const LinOp> &system_matrix,
bool skip_sorting) const;
};


Expand Down
13 changes: 13 additions & 0 deletions include/ginkgo/core/matrix/csr.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,6 +300,19 @@ class Csr : public EnableLinOp<Csr<ValueType, IndexType>>,

std::unique_ptr<LinOp> conj_transpose() const override;

/**
* Sorts all (value, col_idx) pairs in each row by column index
*/
void sort_by_column_idx();

/*
* Tests if all row entry pairs (value, col_idx) are sorted by column index
*
* @returns True if all row entry pairs (value, col_idx) are sorted by
* column index
*/
bool is_sorted_by_column_idx() const;

/**
* Returns the values of the matrix.
*
Expand Down
23 changes: 23 additions & 0 deletions omp/matrix/csr_kernels.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,29 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_CALCULATE_MAX_NNZ_PER_ROW_KERNEL);


template <typename ValueType, typename IndexType>
void sort_by_column_index(std::shared_ptr<const OmpExecutor> exec,
matrix::Csr<ValueType, IndexType> *to_sort)
{
GKO_NOT_IMPLEMENTED;
}

GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_SORT_BY_COLUMN_INDEX);


template <typename ValueType, typename IndexType>
void is_sorted_by_column_index(
std::shared_ptr<const OmpExecutor> exec,
const matrix::Csr<ValueType, IndexType> *to_check, bool *is_sorted)
{
GKO_NOT_IMPLEMENTED;
}

GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
GKO_DECLARE_CSR_IS_SORTED_BY_COLUMN_INDEX);


} // namespace csr
} // namespace omp
} // namespace kernels
Expand Down
Loading

0 comments on commit 2b27fc7

Please sign in to comment.