Chapel: Distributed LinearAlgebra.diag() design

[rahulghangas]

I will add questions as they come along

• [ ] Should the result of diag of a distributed matrix be distributed?

[rahulghangas]

@ben-albrecht need label: gsoc: linearalgebra

[ben-albrecht]

At a high level, I don't think we need to over-optimize this operation, as it's not particularly computationally demanding.

Should the result of diag of a distributed matrix be distributed?

I think we can support both a distributed and localized return vector. That opens a few more questions:

• What should the distribution of the return vector be if distributed?
• How should a user specify the return vector distribution?
• What should the default behavior be?

[bradcray]

For a distributed diagonal matrix, I think we might want a custom domain map that was optimized for that behavior. We could also use a more general one (like a sparse subdomain of the full square matrix), and that might be the best way to proceed in the short-term, but it's likely to be higher-overhead / more general than necessary. Then again, since diagonal operations tend to be asymptotically cheap-ish by nature, maybe investing more effort optimizing for it isn't worthwhile... (or should be put off until the sparse subdomain turns out to have too much overhead).

[ben-albrecht]

What should the distribution of the return vector be if distributed?

For a distributed diagonal matrix, I think we might want a custom domain map that was optimized for that behavior.

Agreed.

We could also use a more general one (like a sparse subdomain of the full square matrix), and that might be the best way to proceed in the short-term, but it's likely to be higher-overhead / more general than necessary

That's an interesting idea.

It seems like there are a lot of reasonable answers to this question, which may depend on the user's needs. I think the interface should support specifying different distributions for the result, but we can keep the initially supported distributions simple for now.

How should a user specify the return vector distribution?

An enum would be clean and minimal, while passing in a distribution object would be the most flexible. We could also potentially support both. I don't have a strong opinion one way or another yet.

@rahulghangas - Maybe writing out some examples would help us decide.

What should the default behavior be?

The LinearAlgebra library always tries to inherit the domain/distribution from the given arrays. Following that philosophy, it could make sense to default to the distribution of the input array (with a dimension removed, of course).

Alternatively, I'd guess that that local vector is the most common use-case, so that would also be a reasonable choice.

[e-kayrakli]

An enum would be clean and minimal, while passing in a distribution object would be the most flexible. We could also potentially support both. I don't have a strong opinion one way or another yet.

I feel like distribution objects are a bit ~second-class citizen~ under-used in the language in general. Maybe it is not because any language limitations and it is just needed that much in real-life, maybe because we don't have a distribution type, maybe those two are related. Anyways, that makes me feel a bit hesitant to pass distribution objects and make users create them etc. I think they can create an array just as simply and pass it as the output array.

[bradcray]

Like @ben-albrecht, my intuition is also that the diagonal matrix should be distributed / follow the original matrix's distribution by default, similar to how a sparse subdomain of a domain D inherits the distribution of D.

[rahulghangas]

@rahulghangas - Maybe writing out some examples would help us decide.

@ben-albrecht I'll try to open a draft PR.

Meanwhile, what about creating a diagonal matrix from a vector? How should the output look like if the input vetor was distributed?

[ben-albrecht]

Meanwhile, what about creating a diagonal matrix from a vector? How should the output look like if the input vetor was distributed?

Following the same principles from before, I think we'll want to inherit from the vector's distribution by default (so if it's a 1D block-distributed vector, we return a 2D block-distributed array). If the user wants to change this distribution they'll have to set an argument. The options discussed so far are:

1. an array, e.g. proc diag(A, ret=B); (proposed by @e-kayrakli)
2. an enum, e.g. proc diag(A, ret=Distributions.BlockDist)
3. a distribution, e.g. proc diag(A, dist=new dmap(new BlockDist(...)))

I've listed these in the order of my preference at the moment. I think passing in the return array is an intuitive way to specify the distribution, and lets users optimize away temporary allocations if they are doing repeated operations. I also agree with @e-kayrakli about distribution objects being under-used. They are also a bit verbose to declare, as demonstrated above.