Add a new permutation matrix format

[pratikvn]

This PR adds a new permutation matrix format, in preparation for the Reordering methods.

Additionally, it also adds a Permutable interface, which allows for LinOps to implement the row_permute, column_permute ,inverse_row_permute and inverse_column_permute functions. The apply method of the permutation matrix format, should automatically call these functions and there should be no need for the user to call these functions.

The Permutation format, holds two permutation arrays, the row_permutation and the col_permutation. It implements methods from EnableLinOp and hence will provide the apply methods. The apply methods will permute the given LinOp by calling the respective LinOp's Xrow_permute/Xcol_permute functions if the LinOp has a Permutable interface.

Along with the Permutable interface, permuting a matrix should look something like this:

using Perm = gko::matrix::Permutation<>;
using Csr = gko::matrix::Csr<>;
using Dense = gko::matrix::Dense<>;

// See the constructors of Permutation for available options
auto p = Perm::create(exec, perm_array);

auto csr_mat=Csr::create(exec, <row_ptrs,  vals ...>);
auto dense_mat=Dense::create(exec, <vals ...>);

auto perm_csr = Csr::create(exec, csr_mat->get_size(), csr_mat->get_num_stored_elements());
auto perm_dense = Dense::create(exec, dense->get_size());

// Permutes the CSR matrix.
p->apply(csr_mat, perm_csr);

// Permutes the dense matrix
p->apply(dense_mat, perm_dense);

TODO

• Merge Sparsity pattern matrix format #349.
• Merge Sparsity format Reference and OpenMP kernels #350.

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (develop@45a800a). Click here to learn what that means.
The diff coverage is 96.17%.

@@            Coverage Diff             @@
##             develop     #352   +/-   ##
==========================================
  Coverage           ?   89.45%           
==========================================
  Files              ?      289           
  Lines              ?    22588           
  Branches           ?        0           
==========================================
  Hits               ?    20205           
  Misses             ?     2383           
  Partials           ?        0

Impacted Files	Coverage Δ
include/ginkgo/core/matrix/dense.hpp	100% <ø> (ø)
include/ginkgo/core/base/exception_helpers.hpp	100% <ø> (ø)
core/matrix/sparsity_csr.cpp	97.18% <ø> (ø)
include/ginkgo/core/matrix/csr.hpp	85.48% <ø> (ø)
reference/matrix/csr_kernels.cpp	100% <100%> (ø)
include/ginkgo/core/base/lin_op.hpp	91.48% <100%> (ø)
core/matrix/csr.cpp	99.43% <100%> (ø)
reference/test/matrix/dense_kernels.cpp	100% <100%> (ø)
core/test/matrix/permutation.cpp	100% <100%> (ø)
omp/matrix/csr_kernels.cpp	97.05% <100%> (ø)
... and 8 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 45a800a...85c98d8. Read the comment docs.

[tcojean]

@pratikvn @thoasm I updated the CUDA docker images to include libthrust-dev, which is required for HIP builds. Could that have messed up the ParILU tests? I don't think the test is complex related but who knows. That's the only thing I can think of right now from my part.

Another thing is, since this code will anyway go away with #353, we may as well not care about it until #353 is merged and everything is rebased on top of it. Debugging dead code is not very useful.

[thoasm]

@tcojean I think it is possible that the problem is similar to what @yhmtsai had on Windows: In the old ILU preconditioner test, I did not set any initial value for the X vector or the cache vector, so it is possible that it reuses parts of previously allocated memory. That could result in wrong computations.
That part is fixed as part of the PR #353, so if that was the problem, rebasing should resolve it (after it is merged in master ofc).

[pratikvn]

Yes, I will leave it alone until #353 is merged and then rebase.

[tcojean]

@pratikvn Nice. As suggestion improvement, you could consider whether it is worth making it similar to precision_reduction , i.e. to turn this into a class eventually with functionalities for easy testing such as is_row_permutation and similar.

ginkgo/include/ginkgo/core/base/types.hpp

Lines 233 to 344 in 3ba4a47

	class precision_reduction {
	public:
	/**
	* The underlying datatype used to store the encoding.
	*/
	using storage_type = uint8;
	private:
	static constexpr auto nonpreserving_bits = 4u;
	static constexpr auto preserving_bits =
	byte_size * sizeof(storage_type) - nonpreserving_bits;
	static constexpr auto nonpreserving_mask =
	storage_type{(0x1 << nonpreserving_bits) - 1};
	static constexpr auto preserving_mask =
	storage_type{(0x1 << preserving_bits) - 1} << nonpreserving_bits;
	public:
	/**
	* Creates a default precision_reduction encoding.
	*
	* This encoding represents the case where no conversions are performed.
	*/
	GKO_ATTRIBUTES constexpr precision_reduction() noexcept : data_{0x0} {}
	/**
	* Creates a precision_reduction encoding with the specified number of
	* conversions.
	*
	* @param preserving the number of range preserving conversion
	* @param nonpreserving the number of range non-preserving conversions
	*/
	GKO_ATTRIBUTES constexpr precision_reduction(
	storage_type preserving, storage_type nonpreserving) noexcept
	: data_((GKO_ASSERT(preserving < (0x1 << preserving_bits) - 1),
	GKO_ASSERT(nonpreserving < (0x1 << nonpreserving_bits) - 1),
	(preserving << nonpreserving_bits) | nonpreserving))
	{}
	/**
	* Extracts the raw data of the encoding.
	*
	* @return the raw data of the encoding
	*/
	GKO_ATTRIBUTES constexpr operator storage_type() const noexcept
	{
	return data_;
	}
	/**
	* Returns the number of preserving conversions in the encoding.
	*
	* @return the number of preserving conversions in the encoding.
	*/
	GKO_ATTRIBUTES constexpr storage_type get_preserving() const noexcept
	{
	return (data_ & preserving_mask) >> nonpreserving_bits;
	}
	/**
	* Returns the number of non-preserving conversions in the encoding.
	*
	* @return the number of non-preserving conversions in the encoding.
	*/
	GKO_ATTRIBUTES constexpr storage_type get_nonpreserving() const noexcept
	{
	return data_ & nonpreserving_mask;
	}
	/**
	* Returns a special encoding which instructs the algorithm to automatically
	* detect the best precision.
	*
	* @return a special encoding instructing the algorithm to automatically
	* detect the best precision.
	*/
	GKO_ATTRIBUTES constexpr static precision_reduction autodetect() noexcept
	{
	return precision_reduction{preserving_mask | nonpreserving_mask};
	}
	/**
	* Returns the common encoding of input encodings.
	*
	* The common encoding is defined as the encoding that does not have more
	* preserving, nor non-preserving conversions than the input encodings.
	*
	* @param x an encoding
	* @param y an encoding
	*
	* @return the common encoding of `x` and `y`
	*/
	GKO_ATTRIBUTES constexpr static precision_reduction common(
	precision_reduction x, precision_reduction y) noexcept
	{
	return precision_reduction(
	min(x.data_ & preserving_mask, y.data_ & preserving_mask) |
	min(x.data_ & nonpreserving_mask, y.data_ & nonpreserving_mask));
	}
	private:
	GKO_ATTRIBUTES constexpr precision_reduction(storage_type data)
	: data_{data}
	{}
	GKO_ATTRIBUTES constexpr static storage_type min(storage_type x,
	storage_type y) noexcept
	{
	return x < y ? x : y;
	}
	storage_type data_;
	};

[pratikvn]

The design looks mostly good to me, a few comments:

1. For column_permute and inverse_row_permute in CSR, the implementation currently uses a linear scan to effectively determine an entry of the inverse permutation for each row/non-zero. We could speed this up considerably by computing the inverse permutation explicitly for the calls or even storing it in the permutation class. This would not make a difference for Dense and only a slight difference for COO (i.e., sorting after permutation), but benefit CSR etc. This would also remove a bit of code duplication from the CSR permute code.

2. Concerning Terry's comment w.r.t the interface, if we look at it from the view of a matrix, we could limit the choices to two options: `inverse_permute` determines if we multiply with `P` or `P*` (transposed), while `row_permute` and `column_permute` describe a multiplication from the left or right. Can you clarify why we need two separate flags for these operations?

3. Just out of curiosity: Do we plan to be able to combine (i.e., multiply) permutation matrices?

@upsj ,

1. Thanks, that is indeed better.

2. row_permute can be seen as P A whereas the column_permute can be seen as A P, so I added two different options.

3. I guess it would be nice to have this as well. But I think as the Permutation class here does not strictly represent a matrix, I think it might be a bit difficult to group two permutations into one.

[upsj]

LGTM! You could remove some duplicate code as the permute and inverse_permute now only differ in whether they use perm[...] or inv_pindx[...].