Improved Matrix organisation

[EmilyBourne]

Currently the super-class Matrix is designed with LAPACK in mind. As a result it defines some pure functions with prototypes that are poorly adapted to other classes:

ddc/include/ddc/kernels/splines/matrix.hpp

Lines 127 to 131 in 9ef820e

	protected:
	virtual int factorize_method() = 0;
	virtual int solve_inplace_method(double* b, char transpose, int n_equations) const = 0;
	};

I think we should reorganise this somewhat. However it is not clear what the correct implementation should be.

In any case the current Matrix class should be renamed to better reflect its purpose. E.g LAPACKMatrix. If it is still useful to have a Matrix abstract class then I think this should look something like:

class Matrix
{
    int m_n;

public:
    explicit Matrix(const int mat_size) : m_n(mat_size) {}

    virtual ~Matrix() = default;

    virtual void factorize() = 0;

    template <class... Args>
    virtual Kokkos::View<double*, Args...> solve_inplace(Kokkos::View<double*, Args...> const b) const = 0;

    template <class... Args>
    virtual Kokkos::View<double*, Args...> solve_transpose_inplace(Kokkos::View<double*, Args...> const b) const = 0;

    template <class... Args>
    virtual Kokkos::View<double**, Args...> solve_multiple_inplace(Kokkos::View<double**, Args...> const bx) const = 0;

    template <class... Args>
    virtual Kokkos::View<double**, Args...> solve_multiple_transpose_inplace(Kokkos::View<double**, Args...> const bx) const = 0;

    template <class... Args>
    Kokkos::View<double**, Args...> solve_batch_inplace(Kokkos::View<double**, Args...> const bx) const = 0;

    int get_size() const
    {
        return m_n;
    }
};

however I believe this code will not compile (templates + virtual). Maybe there is another solution to enforce the implementation of a minimum subset of functions?

One option could be to make Matrix templated by Args... I think this makes sense for MemSpace but I'm not sure it is logical for the layout.

@tpadioleau @jbigot what do you think?

@blegouix any additional suggestions?

[EmilyBourne]

I'm not sure if factorize should be part of the interface. It is specific to certain LAPACK functions but if we want to be able to switch between Matrix implementations then maybe we should ensure that it is callable from all as it will need to be called for LAPACK.

Alternatively we could break up the data and the solver and call factorize at the initialisation of the LAPACK solver. But this would probably involve unnecessary data copies unless the LAPACK constructor consumes the unfactorized matrix

[jbigot]

Could factorize be an implementation detail called from the constructor of child classes? Is there a good reason to separate it?

[blegouix]

Yes we need to fill the matrix first with set_element

[jbigot]

There is no set_element in the API. I expect this to be done in an implementation-specific way, so that factorize could also be implementation-specific. Either

• you provide the content of the matrix in the constructor and it's done behind the scene,
• or it's done in multiple implementation-specific steps and the last one is factorize, but this is just one more implementation-specific step

[blegouix]

The first solution has the drawback to be storage-specific. Atm using Ell storage in place of Csr is one line of code changed in matrix_sparse.hpp. I guess we could pass the prefilled Ginkgo matrix object to the MatrixSparse constructor and we would keep the flexibility though.

However to pre-fill the arrays passed to a MatrixCornerBlock constructor (each submatrix in their own storage) it may be cumbersome. Doable and probably more in the spirit of GPU programming but I would at least prefer to do it after all the matrices have been merged first with the current set/get_element API in order to avoid mistakes.

[tpadioleau]

I agree we should use Kokkos::View. We cannot mix virtual functions and templates so i would template the base class by the execution space and use a fixed layout to right or left. I am not convinced by the layout stride in the end, we kind of need to know the child class to know what layout is expected.

Factorize could be renamed fill_complete like in trilinos ? This way we could keep the same base for sparse and lapack based methods

[EmilyBourne]

I don't see where the stride is necessary in your link. You use a stride but it is extracted from b so if it is disallowed for b then won't u and v also have a stride of 1?

[blegouix]

The left stride of u is 1 indeed (contiguous in left dimension) but the right stride is different from the left extent because you "skip" the values belonging to v

uuuuuuuuuu
uuuuuuuuuu
uuuuuuuuuu
vvvvvvvvvvvv
vvvvvvvvvvvv

[tpadioleau]

Could we get rid of solve_inplace_method that only takes a single contiguous right-hand side ?

[EmilyBourne]

Could we get rid of solve_inplace_method that only takes a single contiguous right-hand side ?

I agree that it seems overkill to have solve_inplace_method and solve_batch_inplace_method

[blegouix]

This is the case in #329, solve_batch_inplace_method is removed and solve_inplace_method takes m-rhs

[EmilyBourne]

Maybe we could write some kind of complicated is_matrix_v template function to enforce the interface

[tpadioleau]

OtherMatrixSolver would need to be the result something like get_adapted_matrix_t<BSplines, BCondL, BCondR>.

I do not see yet the problem with runtime choice ?

[EmilyBourne]

I do not see yet the problem with runtime choice ?

For runtime choice we need to have a superclass but I believe it is not possible to inherit from the class proposed in the discussion description because you cannot have templated virtual functions

[tpadioleau]

What do you think of the following

template <class ExecSpace>
class Matrix
{
    virtual void solve_multiple_inplace(Kokkos::View<double**, typename ExecSpace::memory_space, Kokkos::LayoutRight> const& bx) const = 0;
};

[EmilyBourne]

This looks good to me, but @blegouix will have to confirm if Kokkos::LayoutRight is sufficiently flexible

[blegouix]

I think it's not
#425 (reply in thread)

But LayoutRight or LayoutLeft Kokkos::View should be casted to LayoutStride one if we use:

virtual void solve_multiple_inplace(Kokkos::View<double**, typename ExecSpace::memory_space, Kokkos::LayoutStride> const& bx) const = 0;

right ?

[blegouix]

I think #329 partially solves this by using mdspan (which is equivalent to Kokkos::View but more standard and probably more in the spirit of the rest of DDC, especially Chunk class) . What remains to be done is the LapackMatrix class.

If you want this issue to be addressed before merging #329 I think we can git checkout lapack_matrices -- spline_builder.hpp matrix.hpp matrix_sparse.hpp in a new branch as a starting point.

[EmilyBourne]

I don't know. I imagine it could be useful for other matrix types. You can ask the user to take care of this, e.g. by specifying that the class is in a host exec space so they should always pass host memory, but it is always better to enforce it if possible.
Given that some of our matrices will run on GPU I think it would be nice to have it

[jbigot]

The long-term goal in DDC is to move-away from mdspan and rely more on Kokkos and especially rebase Chunk on Kokkos::view, so I would advise to use View instead of mdspan

[blegouix]

In #329 all Matrix classes except Matrix are templated with ExecSpace so we have the information, so we have the information in principle, But I do not know if you have in mind the problem you mentionned (and that I already encountered by the past) of mixing templates and virtuals.

[blegouix]

ok @jbigot

[EmilyBourne]

In #329 all Matrix classes except Matrix are templated with ExecSpace so we have the information

We have the information about what we want to do but if the user forgets and passes us data in the wrong memory space the code will crash. Using Views the code would simply not compile.

But I do not know if you have in mind the problem you mentionned (and that I already encountered by the past) of mixing templates and virtuals.

For the moment I am just discussing the ideal solution without worrying about the implementation details. We can work out the details when we know what we're aiming for. E.g. as suggested above we could write a is_matrix_v template function instead of using a superclass (e..g using https://en.cppreference.com/w/cpp/types/is_invocable)

[jbigot]

First, let me say I like the API, it is nice and clear.

Maybe more than Matrix, we could see this as a LinearProblemSolver operator, it happens to be representable as a matrix, but here we don't provide everything a matrix can do, only the linear problem resolution part.

Regarding inheritance, do we really need it? Couldn't we just converge on a templated API like the one you propose and have multiple independent implementations of it. It would then be up to the user to be general and template or to make a virtual-based wrapper with specific types of parameters

[jbigot]

And as you mention in a comment above, provide a C++20 Concept to enforce the interface

[tpadioleau]

Isn't DDC the only user ? I do not see the benefit of using concepts here, it seems quite rigid.

[EmilyBourne]

Isn't it the same rigidity as having an abstract superclass? In both cases it is still possible to add additional functions, it just enforces that the minimum requirements are met

[EmilyBourne]

Assembling some of the comments here to show a new proposal:

template <class ExecSpace>
class Matrix
{
    int m_n;

    using DView2D = Kokkos::View<double**, typename ExecSpace::memory_space, Kokkos::LayoutStride>;

public:
    explicit Matrix(const int mat_size) : m_n(mat_size) {}

    virtual ~Matrix() = default;

    virtual DView2D solve_inplace(DView2D const b) const = 0;

    virtual DView2D solve_transpose_inplace(DView2D const b) const = 0;

    virtual DView2D solve_multiple_inplace(DView2D const bx) const = 0;

    virtual DView2D solve_multiple_transpose_inplace(DView2D const bx) const = 0;

    int get_size() const
    {
        return m_n;
    }
};

[jbigot]

We discussed with Thomas. The goal is

1. to simplify to the maximum and fit the needs of the usage in DDC splines, no more,
2. target GPU (or more exactly DefaultExecutionSpace), anything CPU-only is a placeholder that we should not try to optimize and doing useless copies, transpose, etc., for these cases is completely fine.

Here's what we came up width:

1. make the constructor protected,
2. keep set_element and finished_filling (in replacement of factorize) in the base API to simplify factorization of the matrix filling between various implems,
3. Do not parameterize based on Execution or Memory space, offer a DefaultMemorySpace and DefaultExecutionSpace based API for now, only if more is required should it be provided,
4. use LayoutRight or LayoutLeft, see point above regarding providing more flexibility,
5. take a DefaultExecutionSpace instance as first parameter, fence according to it, so that if execution is done on a different space (e.g. CPU-only LAPACK), then the data is available for a following kernel on the provided space instance, but no need to fence if the execution can be done on this instance directly,
6. see if the non batched (multiple) versions of solve are really required,
7. see if the transpose should be specified through different functions or through a boolean parameter (here is your bikeshedding point, have fun ;p).

@tpadioleau let me know if I mispresented our discussion

[EmilyBourne]

This sounds good to me. But I am not sure I followed all the parametrisation suggestions. Can you share the suggested code? In particular I am not clear on how you want to handle DefaultExecutionSpace and whether this means that the code contains templated functions or if it retains virtual methods.

Regarding point 6. "non batched (multiple)" (solve_multiple_inplace) is effectively batching over the RHS. This is used in the splines. The alternative (solve_inplace) is batched over the matrices. In Gysela this is used for the collisions. I think both are useful. They could be combined into 1 method with a dimension for each type of batch but care would need to be taken to clearly document which batch dimension corresponds to which batch type. I think I've got this mixed up. I agree that there doesn't need to be 2 methods

For 7, I did wonder if the transpose should be restricted to the LAPACK implementations for now. I'm not sure how common it is to need this or how easy it is to implement for some of the other solvers

[jbigot]

This sounds good to me. But I am not sure I followed all the parametrisation suggestions

None :)

[jbigot]

Can you share the suggested code?

class LinearProblem
{
public:
    using MultiRHS = Kokkos::View<double**, Kokkos::DefaultExecutionSpace>;
    
private:
    size_t m_size;

protected:
    explicit LinearProblem(size_t const size) : m_size(size) {}

public:
    virtual ~LinearProblem() = default;

    virtual void set_element(size_t row, size_t column, double value) = 0;

    virtual void finished_filling() const = 0;

    virtual void solve(Kokkos::DefaultExecutionSpace const& exec_space, MultiRHS const& multi_rhs, bool transpose = false) const = 0;

    size_t size() const { return m_n; }

};

[tpadioleau]

@tpadioleau let me know if I mispresented our discussion

Well I would still go for the template version because the spline builder is template by the execution space

template <class ExecSpace>
class SplinesLinearProblem
{
public:
    using MultiRHS = Kokkos::View<double**, typename ExecSpace::memory_space, SomeLayout>;
    
private:
    std::size_t m_size;

protected:
    explicit SplinesLinearProblem(std::size_t const size) : m_size(size) {}

public:
    virtual ~SplinesLinearProblem() = default;

    virtual void set_element(std::size_t row, std::size_t column, double value) = 0;

    virtual void finished_filling() const = 0;

    virtual void solve(ExecSpace const& exec_space, MultiRHS const& multi_rhs, bool transpose = false) const = 0;

    std::size_t size() const { return m_n; }
};

Calling set_element after finished_filling is an error.