Refactor algorithm selection in type domain

[lkdvos]

This alters the algorithm selection and default algorithms to be implemented in the type domain instead.
While doing this, I also realized that the new default_algorithm(f, ...) and default_f_alg are tautologous, so I removed the latter.
Given that these were not public, and the implementations in the "value domain" still dispatch to the type domain, I think this is non-breaking, so should be v0.2.1.

Fixes #29.

[codecov]

Codecov Report

Attention: Patch coverage is 73.68421% with 20 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/algorithms.jl	90.62%	3 Missing ⚠️
src/interface/eigh.jl	62.50%	3 Missing ⚠️
src/interface/lq.jl	50.00%	3 Missing ⚠️
src/interface/polar.jl	57.14%	3 Missing ⚠️
src/interface/qr.jl	50.00%	3 Missing ⚠️
src/interface/svd.jl	62.50%	3 Missing ⚠️
src/interface/eig.jl	71.42%	2 Missing ⚠️

Files with missing lines	Coverage Δ
src/interface/schur.jl	33.33% <100.00%> (-16.67%)	⬇️
src/yalapack.jl	67.69% <ø> (ø)
src/interface/eig.jl	53.84% <71.42%> (-19.49%)	⬇️
src/algorithms.jl	82.92% <90.62%> (+0.49%)	⬆️
src/interface/eigh.jl	50.00% <62.50%> (-23.34%)	⬇️
src/interface/lq.jl	57.14% <50.00%> (-42.86%)	⬇️
src/interface/polar.jl	41.66% <57.14%> (-21.97%)	⬇️
src/interface/qr.jl	57.14% <50.00%> (-42.86%)	⬇️
src/interface/svd.jl	50.00% <62.50%> (-23.34%)	⬇️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mtfishman]

Looks good to me, this will be helpful.

[Jutho]

This alters the algorithm selection and default algorithms to be implemented in the type domain instead.

While I am not opposed to this change, I was wondering if this helps with anything? Semantically, I guess it makes sense that the default algorithm only depends on the type of A and not on the values, but that was what you anyway do with multiple dispatch if you define new methods for a given combination of types, in this case mainly the type of A.

Furthermore, there might actually be cases where the default algorithm selection does want to depend on properties of A that are not encoded in the type, e.g. isposdef or ishermitian (I don't like the Hermitian wrapper from LinearAlgebra.jl very much, and we don't have it for TensorMap objects.)

While doing this, I also realized that the new default_algorithm(f, ...) and default_f_alg are tautologous, so I removed the latter.

Not really, as the default_f_algorithm was another level of where to specify default algorithms in a way that affected all methods related to the decomposition in f, so e.g. all of qr_full(!), qr_compact(!) and qr_null(!). While it was perhaps not yet declared as a public and documented method, I did find it a useful entry point, for example in the (WIP) CUDA extension. It is an easy enough change to the new style of this PR, but it looked natural that one would often want to specify a new default for all methods related to the same mathematical decomposition.

[lkdvos]

The main point where this appeared was in trying to select a default algorithm for either TensorMap or BlockSparseArrays: I wanted to reuse the default_algorithm that is defined for the individual blocks, but these might not be present (an annoying edge-case is for example an empty block-sparse array), or not so easy to obtain.

I somehow assumed this solution was the easiest, since both TensorMap and BlockSparseArrays can quite easily tell you their blocktype, but an alternative solution is to simply define a DefaultAlgorithm dummy type to effectively delay choosing the algorithm type, similar to how TensorOperations handles this. I'm definitely open to either (or both?).
In principle currently the code should still work in the value domain, so if a type really wants to do runtime checks in order to alter the default algorithm selection that could still be possible.

I'll reintroduce the other default_... functions, I see what you mean.

[mtfishman]

EDIT: Sorry for repeating your points you Lukas, I posted this before noticing your post.

This alters the algorithm selection and default algorithms to be implemented in the type domain instead.

While I am not opposed to this change, I was wondering if this helps with anything? Semantically, I guess it makes sense that the default algorithm only depends on the type of A and not on the values, but that was what you anyway do with multiple dispatch if you define new methods for a given combination of types, in this case mainly the type of A.

The use case we have in mind is a block diagonal matrix, where you might want to extract the type of the block to determine the default algorithm for the blocks, i.e.:

default_algorithm(f::typeof(svd_compact!), A::BlockDiagonalMatrix) = BlockDiagonalAlgorithm(default_algorithm(f, blocktype(A))

It's a bit easier to reason about working in the type domain, rather than grabbing instances of blocks.

Furthermore, there might actually be cases where the default algorithm selection does want to depend on properties of A that are not encoded in the type, e.g. isposdef or ishermitian (I don't like the Hermitian wrapper from LinearAlgebra.jl very much, and we don't have it for TensorMap objects.)

That's true, but from what I can tell this PR still allows customizing based on instances if you want to.

While doing this, I also realized that the new default_algorithm(f, ...) and default_f_alg are tautologous, so I removed the latter.

Not really, as the default_f_algorithm was another level of where to specify default algorithms in a way that affected all methods related to the decomposition in f, so e.g. all of qr_full(!), qr_compact(!) and qr_null(!). While it was perhaps not yet declared as a public and documented method, I did find it a useful entry point, for example in the (WIP) CUDA extension. It is an easy enough change to the new style of this PR, but it looked natural that one would often want to specify a new default for all methods related to the same mathematical decomposition.

I agree with this, I think it is a bit easier to reason about code like:

default_svd_algorithm(A) = LAPACK_DivideAndConquer()
default_algorithm(::typeof(svd_compact!), A) = default_svd_algorithm(A)
default_algorithm(::typeof(svd_full!), A) = default_svd_algorithm(A)

rather than:

default_algorithm(::typeof(svd_compact!), A) = LAPACK_DivideAndConquer()
default_algorithm(::typeof(svd_full!), A) = default_algorithm(svd_compact!, A)

since you can more easily remember that the actual definition is default_svd_algorithm, rather than having a convention that svd_compact! is the "core" definition, which could be documented but is a bit arbitrary.

[mtfishman]

I somehow assumed this solution was the easiest, since both TensorMap and BlockSparseArrays can quite easily tell you their blocktype, but an alternative solution is to simply define a DefaultAlgorithm dummy type to effectively delay choosing the algorithm type, similar to how TensorOperations handles this. I'm definitely open to either (or both?). In principle currently the code should still work in the value domain, so if a type really wants to do runtime checks in order to alter the default algorithm selection that could still be possible.

A DefaultAlgorithm dummy type may be best anyway, for example that could handle the case where there are different types of blocks. In that case, the block type wouldn't be that helpful, since it would be some abstract supertype of the block types.

[Jutho]

I see the argument about the type domain, that is indeed a very valid remark. The lack of certain Base functionality in the type domain (e.g. similar) has indeed been very annoying to work around, so I am all in with making default_algorithm operate in the type domain by default.

No strong opinions on the default_f_algorithm convention, I just wanted to point out the potential usefulness, but I am happy to adapt to either choice.

[lkdvos]

I've reinstated the desired functions and opened an issue to remind myself to add the default algorithm, which I'll implement in a separate PR.

[mtfishman]

@Jutho do you have any more comments on this? I think it would be nice to merge this, I have some ongoing work that would benefit from this.

[Jutho]

Yes looks good to me, feel free to merge.