`lu!`: `DivideError: integer division error` for Arrays of `ForwardDiff.Dual`s with large chunkSize

[dafred94]

I am experiencing an error DivideError: integer division error when trying to apply lu! to a Matrix of ForwardDiff.Duals with a large chunk size.

After trying to reduce the problem as far as I could, I believe that the error is somewhere in this Package (probably in nsplit). See MWE below.

Note that when both matWidth and chunkSize are set to 25, the error does not happen. However a different error is then thrown (SingularException(9)) but I believe that this is due to the nature of f_test always returning a 0.

using ForwardDiff
using RecursiveFactorization
using LinearAlgebra

# DivideError
matWidth = 256
chunkSize = 256

# (different error with this parametrization)
# matWidth = 25
# chunkSize = 25

dualBuffer = 0


function f_test(x)
    global dualBuffer
    dualBuffer = copy(x)
    return zero(eltype(x))
end 


x = reshape(convert.(Float64, 1:matWidth^2), (matWidth, matWidth))

cfg = ForwardDiff.GradientConfig(f_test, x, ForwardDiff.Chunk{chunkSize}());
ForwardDiff.gradient(f_test, x, cfg)

RecursiveFactorization.lu!(dualBuffer)

Am I doing something illegal here or do you agree that the DivideError should not pop up in the first place?

[chriselrod]

If you look at the definition of nsplit https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl/blob/ae6debe2d733a91966f11c6dfddaf6d6f4317f1d/src/lu.jl#L103-L107 You'll see that a very large sizeof(T) (i.e., > 512) can make k = 0, resulting in an integer division error with (n + k_2) ÷ k.

We should probably instead define

     k = max(2, 512 ÷ (isbitstype(T) ? sizeof(T) : 8))

Also, I should probably split out the ForwardDiff.Dual multiplication functions from https://github.com/PumasAI/SimpleChains.jl/blob/main/src/forwarddiff_matmul.jl into their own repo, and then we can depend on it here for better/faster ForwardDiff.Dual support.

That code will need modification to support C -= A*B, but that should be straightforward.

EDIT: I guess Octavian.jl already is that library, but it is a bit hefty as a dependency...

[chriselrod]

Also, differentiating LU should be totally unnecessary. We can/will get around to adding dual overloads to linear solve.

[dafred94]

Also, differentiating LU should be totally unnecessary. We can/will get around to adding dual overloads to linear solve.

We like this idea! Any news on this? Or should we help somehow?