Skip to content

Conversation

oscardssmith
Copy link
Member

This is potentially much faster since it doesn't have to initialize values. For constructing a 1000x1000 Matrix{Float64} I benchmark this to be 10x faster which is very important because in cases like an autoswitching solver, you often want to preallocate a jacobian matrix, but might never use it if the equation isn't stiff.

This is potentially much faster since it doesn't have to initialize values. For constructing a 1000x1000 `Matrix{Float64}` I benchmark this to be 10x faster which is very important because in cases like an autoswitching solver, you often want to preallocate a jacobian matrix, but might never use it if the equation isn't stiff.
@codecov
Copy link

codecov bot commented Oct 10, 2022

Codecov Report

Merging #353 (2b2e40c) into master (aa1fcd2) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##           master     #353   +/-   ##
=======================================
  Coverage   90.27%   90.27%           
=======================================
  Files           9        9           
  Lines        1337     1337           
=======================================
  Hits         1207     1207           
  Misses        130      130           
Impacted Files Coverage Δ
src/ArrayInterface.jl 85.96% <ø> (ø)

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

oscardssmith pushed a commit to oscardssmith/OrdinaryDiffEq.jl that referenced this pull request Oct 10, 2022
JuliaArrays/ArrayInterface.jl#353 adds an `undefmatrix`which is significantly faster than `zeromatrix` since it doesn't have to initialize values. This matters for autostiff solvers on non-stiff equations since in these cases filling the jacobian matrix with zeros can be as expensive as the entire ODE solve.
@oscardssmith
Copy link
Member Author

Should we merge this (so i can make the PR that adds this to ODE solvers)?

@ChrisRackauckas
Copy link
Member

This doesn't get the type correct for static arrays, and it doesn't have tests.

@oscardssmith
Copy link
Member Author

What should I be testing here? Array, CuArray, SArray?

@oscardssmith
Copy link
Member Author

Ah the first version worked with static array and sparse (and the performance difference is only there for large values anyway).

@oscardssmith
Copy link
Member Author

I've added tests for base types. The StaticArrays issue is arguably their own fault. They haven't defined a form of similar which it really feels like is their job.

@chriselrod
Copy link
Collaborator

I've added tests for base types. The StaticArrays issue is arguably their own fault. They haven't defined a form of similar which it really feels like is their job.

The problem is probably the call to ArrayInterface.length, which either is overloaded to return StaticInts, or returns regular Ints.
If it is the latter, they're regular Ints and thus it'd risk type instabilities to have similar(::Type{<:StaticArray}, ::Int, ::Int) return a StaticArray.
If the former, yeah, it's probably StaticArray's or ArrayInterfaceStaticArray's job to define this overload.

@ChrisRackauckas
Copy link
Member

similar is defined in Base to return a mutable array.

help?> similar
search: similar

  similar(array, [element_type=eltype(array)], [dims=size(array)])

  Create an uninitialized mutable array with the given element type and size, based upon the given source array. The       
  second and third arguments are both optional, defaulting to the given array's eltype and size. The dimensions may be
  specified either as a single tuple argument or as a series of integer arguments.

  Custom AbstractArray subtypes may choose which specific array type is best-suited to return for the given element type   
  and dimensionality. If they do not specialize this method, the default is an Array{element_type}(undef, dims...).        

  For example, similar(1:10, 1, 4) returns an uninitialized Array{Int,2} since ranges are neither mutable nor support 2    
  dimensions:

@chriselrod
Copy link
Collaborator

By StaticArray I meant the abstract type for which MArray <: StaticArray, not the immutable SArray.
Creating an undef SArray is useless.

Creating an MArray would be type unstable given Int args, while StaticInt would require a package to implement it.
It'd be type piracy if ArrayInterfaceStaticArrays does it, but that's probably okay.

@ChrisRackauckas
Copy link
Member

Add it to the docs

@ChrisRackauckas
Copy link
Member

Another test failure.

@ChrisRackauckas
Copy link
Member

Bro.

@oscardssmith
Copy link
Member Author

fair. This time I've actually guaranteed that tests pass on my computer.

@ChrisRackauckas ChrisRackauckas merged commit a5e8275 into JuliaArrays:master Oct 13, 2022
@oscardssmith oscardssmith deleted the patch-1 branch October 13, 2022 00:23
Comment on lines +12 to +19
function ArrayInterface.undefmatrix(::MArray{S, T, N, L}) where {S, T, N, L}
return MMatrix{L, L, T, L*L}(undef)
end
# SArray doesn't have an undef constructor and is going to be small enough that this is fine.
function ArrayInterface.undefmatrix(s::SArray)
v = vec(s)
return v.*v'
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These could probably be staticarrayscore?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants