Allocate df with similar, partially fixes GPU support #123

charleskawczynski · 2020-06-06T21:27:40Z

This PR allocates df using similar instead of eltype to support CuArrays.

Script: test/runtests_gpu.jl

using NLsolve
using LinearAlgebra
using CuArrays
# CuArrays.allowscalar(false)  # toggling this

function f!(F, x)
    mul!(F,A,x)
    F .-= R
end

A = cu(rand(2,2))
R = cu(rand(2))
x = cu([ 0.1; 1.2])
nlsolve(f!, x)

Before (without CuArrays.allowscalar(false)):

julia> include("test\\runtests_gpu.jl")
┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:43
ERROR: LoadError: ArgumentError: cannot take the CPU address of a CuArray{Float32,1,Nothing}
Stacktrace:
 [1] dogleg!(::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::Array{Float32,2}, ::Float32) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:61
 [2] trust_region_(::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Float32, ::Bool, ::NLsolve.NewtonTrustRegionCache{CuArray{Float32,1,Nothing}}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:167
 [3] trust_region at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:229 [inlined] (repeats 2 times)
 [4] #nlsolve#25(::Symbol, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Static, ::NLsolve.var"#27#29", ::Float32, ::Bool, ::Int64, ::Int64, ::Int64, ::Float32, ::typeof(nlsolve), ::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:26
 [5] #nlsolve at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:0 [inlined]
 [6] #nlsolve#30(::Symbol, ::Symbol, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(nlsolve), ::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:52
 [7] nlsolve(::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:46
 [8] top-level scope at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14
 [9] include at .\boot.jl:328 [inlined]
 [10] include_relative(::Module, ::String) at .\loading.jl:1105
 [11] include(::Module, ::String) at .\Base.jl:31
 [12] include(::String) at .\client.jl:424
 [13] top-level scope at REPL[1]:1
in expression starting at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14
caused by [exception 1]
ArgumentError: cannot take the CPU address of a CuArray{Float32,1,Nothing}
Stacktrace:
 [1] unsafe_convert(::Type{Ptr{Float32}}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\.julia\packages\CuArrays\l0gXB\src\array.jl:226
 [2] getrs!(::Char, ::Array{Float32,2}, ::Array{Int64,1}, ::CuArray{Float32,1,Nothing}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\lapack.jl:1017
 [3] ldiv!(::LU{Float32,Array{Float32,2}}, ::CuArray{Float32,1,Nothing}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\lu.jl:391
 [4] \(::LU{Float32,Array{Float32,2}}, ::CuArray{Float32,1,Nothing}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\factorization.jl:99
 [5] \(::Array{Float32,2}, ::CuArray{Float32,1,Nothing}) at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.3\LinearAlgebra\src\generic.jl:1050
 [6] dogleg!(::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::CuArray{Float32,1,Nothing}, ::Array{Float32,2}, ::Float32) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:51
 [7] trust_region_(::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Float32, ::Bool, ::NLsolve.NewtonTrustRegionCache{CuArray{Float32,1,Nothing}}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:167
 [8] trust_region at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:229 [inlined] (repeats 2 times)
 [9] #nlsolve#25(::Symbol, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Static, ::NLsolve.var"#27#29", ::Float32, ::Bool, ::Int64, ::Int64, ::Int64, ::Float32, ::typeof(nlsolve), ::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:26
 [10] #nlsolve at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:0 [inlined]
 [11] #nlsolve#30(::Symbol, ::Symbol, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(nlsolve), ::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:52
 [12] nlsolve(::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:46
 [13] top-level scope at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14
 [14] include at .\boot.jl:328 [inlined]
 [15] include_relative(::Module, ::String) at .\loading.jl:1105
 [16] include(::Module, ::String) at .\Base.jl:31
 [17] include(::String) at .\client.jl:424
 [18] top-level scope at REPL[1]:1

After (without CuArrays.allowscalar(false)):

julia> include("test\\runtests_gpu.jl")
[ Info: Precompiling NLsolve [2774e3e8-f4cf-5e23-947b-6d7e65073b56]
┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:43
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: Float32[0.1, 1.2]
 * Zero: Float32[1.8909574, -1.4453843]
 * Inf-norm of residuals: 0.000000
 * Iterations: 1000
 * Convergence: false
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: false
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 4

Before (with CuArrays.allowscalar(false)):

julia> include("test\\runtests_gpu.jl")
ERROR: LoadError: scalar getindex is disallowed
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] assertscalar(::String) at C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:41
 [3] getindex(::CuArray{Float32,1,Nothing}, ::Int64) at C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:96
 [4] _broadcast_getindex at .\broadcast.jl:596 [inlined]
 [5] _getindex at .\broadcast.jl:626 [inlined]
 [6] _broadcast_getindex at .\broadcast.jl:602 [inlined]
 [7] _getindex at .\broadcast.jl:626 [inlined]
 [8] _broadcast_getindex at .\broadcast.jl:602 [inlined]
 [9] getindex at .\broadcast.jl:563 [inlined]
 [10] macro expansion at .\broadcast.jl:909 [inlined]
 [11] macro expansion at .\simdloop.jl:77 [inlined]
 [12] copyto! at .\broadcast.jl:908 [inlined]
 [13] copyto! at .\broadcast.jl:863 [inlined]
 [14] materialize! at .\broadcast.jl:822 [inlined]
 [15] #finite_difference_jacobian!#40(::Float32, ::Float32, ::UnitRange{Int64}, ::Nothing, ::Bool, ::typeof(FiniteDiff.finite_difference_jacobian!), ::Array{Float32,2}, ::typeof(f!), ::CuArray{Float32,1,Nothing}, ::FiniteDiff.JacobianCache{CuArray{Float32,1,Nothing},CuArray{Float32,1,Nothing},CuArray{Float32,1,Nothing},UnitRange{Int64},Nothing,Val{:central},Float32}, ::Nothing) at C:\Users\kawcz\.julia\packages\FiniteDiff\hcMgi\src\jacobians.jl:371
 [16] finite_difference_jacobian! at C:\Users\kawcz\.julia\packages\FiniteDiff\hcMgi\src\jacobians.jl:295 [inlined] (repeats 2 times)
 [17] (::NLSolversBase.var"#fj_finitediff!#21"{typeof(f!),FiniteDiff.JacobianCache{CuArray{Float32,1,Nothing},CuArray{Float32,1,Nothing},CuArray{Float32,1,Nothing},UnitRange{Int64},Nothing,Val{:central},Float32}})(::CuArray{Float32,1,Nothing}, ::Array{Float32,2}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLSolversBase.jl\src\objective_types\oncedifferentiable.jl:139
 [18] value_jacobian!!(::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}, ::Array{Float32,2}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLSolversBase.jl\src\interface.jl:124
 [19] value_jacobian!! at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLSolversBase.jl\src\interface.jl:122 [inlined]
 [20] trust_region_(::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Float32, ::Bool, ::NLsolve.NewtonTrustRegionCache{CuArray{Float32,1,Nothing}}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:119
 [21] trust_region at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:229 [inlined] (repeats 2 times)
 [22] #nlsolve#25(::Symbol, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Static, ::NLsolve.var"#27#29", ::Float32, ::Bool, ::Int64, ::Int64, ::Int64, ::Float32, ::typeof(nlsolve), ::OnceDifferentiable{CuArray{Float32,1,Nothing},Array{Float32,2},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:26
 [23] #nlsolve at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:0 [inlined]
 [24] #nlsolve#30(::Symbol, ::Symbol, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(nlsolve), ::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:52
 [25] nlsolve(::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:46
 [26] top-level scope at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14
 [27] include at .\boot.jl:328 [inlined]
 [28] include_relative(::Module, ::String) at .\loading.jl:1105
 [29] include(::Module, ::String) at .\Base.jl:31
 [30] include(::String) at .\client.jl:424
 [31] top-level scope at REPL[1]:1
in expression starting at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14

After (with CuArrays.allowscalar(false)):

julia> include("test\\runtests_gpu.jl")
ERROR: LoadError: scalar setindex! is disallowed
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] assertscalar(::String) at C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:41
 [3] setindex!(::CuArray{Float32,1,Nothing}, ::Float32, ::Int64) at C:\Users\kawcz\.julia\packages\GPUArrays\JqOUg\src\host\indexing.jl:103
 [4] trust_region_(::OnceDifferentiable{CuArray{Float32,1,Nothing},CuArray{Float32,2,Nothing},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Float32, ::Bool, ::NLsolve.NewtonTrustRegionCache{CuArray{Float32,1,Nothing}}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:147
 [5] trust_region at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\solvers\trust_region.jl:229 [inlined] (repeats 2 times)
 [6] #nlsolve#25(::Symbol, ::Float32, ::Float32, ::Int64, ::Bool, ::Bool, ::Bool, ::Static, ::NLsolve.var"#27#29", ::Float32, ::Bool, ::Int64, ::Int64, ::Int64, ::Float32, ::typeof(nlsolve), ::OnceDifferentiable{CuArray{Float32,1,Nothing},CuArray{Float32,2,Nothing},CuArray{Float32,1,Nothing}}, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:26
 [7] #nlsolve at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:0 [inlined]
 [8] #nlsolve#30(::Symbol, ::Symbol, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(nlsolve), ::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:52
 [9] nlsolve(::Function, ::CuArray{Float32,1,Nothing}) at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\src\nlsolve\nlsolve.jl:46
 [10] top-level scope at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14
 [11] include at .\boot.jl:328 [inlined]
 [12] include_relative(::Module, ::String) at .\loading.jl:1105
 [13] include(::Module, ::String) at .\Base.jl:31
 [14] include(::String) at .\client.jl:424
 [15] top-level scope at REPL[1]:1
in expression starting at C:\Users\kawcz\Dropbox\Caltech\work\dev\clones\NLsolve.jl\test\runtests_gpu.jl:14

longemen3000

Are there any other places where this can be added? It seems like a good solution

pkofod · 2020-06-10T05:14:04Z

Thanks, I seem to remember that I had similar and went away from it for some very particular reason, but I cannot remember why. If I can’t think of why, I’ll merge :) but I’ll give optim and nlsolve a spin with it first.

pkofod · 2020-06-10T09:29:49Z

Yup, I think I remember. It's the ArrayPartition test in Optim... I have to look into this. Can you check if this doesn't fail for you with your version of NLSolversBase and latest Optim?

using RecursiveArrayTools
@testset "arraypartition input" begin

    function polynomial(x)
            return (10.0 - x[1])^2 + (7.0 - x[2])^4 + (108.0 - x[3])^4
        end

    function polynomial_gradient!(storage, x)
            storage[1] = -2.0 * (10.0 - x[1])
            storage[2] = -4.0 * (7.0 - x[2])^3
            storage[3] = -4.0 * (108.0 - x[3])^3
        end

    function polynomial_hessian!(storage, x)
            storage[1, 1] = 2.0
            storage[1, 2] = 0.0
            storage[1, 3] = 0.0
            storage[2, 1] = 0.0
            storage[2, 2] = 12.0 * (7.0 - x[2])^2
            storage[2, 3] = 0.0
            storage[3, 1] = 0.0
            storage[3, 2] = 0.0
            storage[3, 3] = 12.0 * (108.0 - x[3])^2
        end

    ap = ArrayPartition(rand(1), rand(2))

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, NelderMead())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, ParticleSwarm())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, SimulatedAnnealing())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, GradientDescent())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, AcceleratedGradientDescent())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, MomentumGradientDescent())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, ConjugateGradient())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, BFGS())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, LBFGS())

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, Newton())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, NewtonTrustRegion())
end

I'm in the middle of something else right now, but maybe we can specifically fix this version by branching out if all(xi -> length(xi) == 1, x) or something similar.

pkofod · 2020-06-10T17:38:30Z

Or maybe this is to be considered a bug in RecursiveArrayTools...

charleskawczynski · 2020-06-10T17:41:25Z

Are there any other places where this can be added? It seems like a good solution

Potentially the other alloc_DF, but I'd have to test it.

Yup, I think I remember. It's the ArrayPartition test in Optim... I have to look into this. Can you check if this doesn't fail for you with your version of NLSolversBase and latest Optim?

The test you shared does fail in Optim.jl with these changes. My guess is that the existing test relies on the increased precision for the gradient for Float32 tests. Looking into this now...

pkofod · 2020-06-10T17:43:37Z

The test you shared does fail in Optim.jl with these changes. My guess is that the existing test relies on the increased precision for the gradient for Float32 tests. Looking into this now...

What do you mean here? Can you share the error?

charleskawczynski · 2020-06-10T17:47:36Z

Or maybe this is to be considered a bug in arrayapartitions...

Maybe this is it

What do you mean here?

I thought I was getting NaNs in gradients, potentially from reduced precision but I think I was mistaken (I see the NaNs in the existing passing tests). I think the issue is somewhere else.

Can you share the error?

    11     1.000000e+00              NaN
 * time: 0.045770883560180664
 * x: [0.0, 0.0]
  0.138849 seconds (795.50 k allocations: 19.416 MiB)
./multivariate/array.jl
arraypartition input: Error During Test at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/multivariate/array.jl:57
  Got exception outside of a @test
  BoundsError: attempt to access 1-element Array{Float64,1} at index [2]
  Stacktrace:
   [1] throw_boundserror(::Array{Float64,1}, ::Tuple{Int64}) at ./abstractarray.jl:538
   [2] checkbounds at ./abstractarray.jl:503 [inlined]
   [3] setindex! at /Users/charliekawczynski/.julia/packages/RecursiveArrayTools/GrX6g/src/array_partition.jl:193 [inlined]
   [4] (::var"#polynomial_hessian!#315")(::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/multivariate/array.jl:71
   [5] hessian!!(::TwiceDifferentiable{Float64,ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}}, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/forks/NLSolversBase.jl/src/interface.jl:95
   [6] initial_state(::Newton{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}}}, ::Optim.Options{Float64,Nothing}, ::TwiceDifferentiable{Float64,ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}}, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/src/multivariate/solvers/second_order/newton.jl:46
   [7] optimize(::TwiceDifferentiable{Float64,ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}},ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}}, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}, ::Newton{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}}}, ::Optim.Options{Float64,Nothing}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/src/multivariate/optimize/optimize.jl:33
   [8] #optimize#96(::Bool, ::Symbol, ::typeof(optimize), ::Function, ::Function, ::Function, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}, ::Newton{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}}}, ::Optim.Options{Float64,Nothing}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/src/multivariate/optimize/interface.jl:136
   [9] optimize(::Function, ::Function, ::Function, ::ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}, ::Newton{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}}}, ::Optim.Options{Float64,Nothing}) at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/src/multivariate/optimize/interface.jl:134 (repeats 2 times)
   [10] top-level scope at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/multivariate/array.jl:96
   [11] top-level scope at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
   [12] top-level scope at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/multivariate/array.jl:59
   [13] include at ./boot.jl:328 [inlined]
   [14] include_relative(::Module, ::String) at ./loading.jl:1105
   [15] include(::Module, ::String) at ./Base.jl:31
   [16] include(::String) at ./client.jl:424
   [17] top-level scope at util.jl:155
   [18] top-level scope at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/runtests.jl:251
   [19] top-level scope at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/Test/src/Test.jl:1107
   [20] top-level scope at /Users/charliekawczynski/Dropbox/Caltech/work/dev/clones/Optim.jl/test/runtests.jl:249
   [21] include at ./boot.jl:328 [inlined]
   [22] include_relative(::Module, ::String) at ./loading.jl:1105
   [23] include(::Module, ::String) at ./Base.jl:31
   [24] include(::String) at ./client.jl:424
   [25] top-level scope at none:6
   [26] eval(::Module, ::Any) at ./boot.jl:330
   [27] exec_options(::Base.JLOptions) at ./client.jl:263
   [28] _start() at ./client.jl:460
  
 20.424594 seconds (52.30 M allocations: 2.449 GiB, 4.76% gc time)
./multivariate/extrapolate.jl

pkofod · 2020-06-10T18:25:28Z

Yeah, it appears that similar doesn't correctly produce the approriate arraypartition vector. Maybe this is simply a bug in RecursiveArrayTools.jl.

charleskawczynski · 2020-06-10T18:32:25Z

The only tests failing are

    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, Newton())
    optimize(polynomial, polynomial_gradient!, polynomial_hessian!, ap, NewtonTrustRegion())

Others pass fine. So, I think there's an issue in how TwiceDifferentiable is initialized

pkofod · 2020-06-10T19:12:49Z

Interesting, let me see.

charleskawczynski · 2020-06-10T19:16:14Z

Printing info when there are differences before and after the changes is telling:

function alloc_DF(x, F)
  a = (Base.OneTo(length(F)), Base.OneTo(length(x)))
  df = similar(F, a)
  fill!(df, NaN)
  df_old = fill(eltype(F)(NaN), length(F), length(x))
  if !(size(df) == size(df_old))
    println("------------------- ")
    @show F
    @show x
    @show typeof(F)
    @show typeof(x)
    @show eltype(F)
    @show eltype(x)
    @show df
    @show df_old
    @show typeof(df)
    @show typeof(df_old)
    @show eltype(df)
    @show eltype(df_old)
    println("-------------------")
  end
  return df
end

yields

------------------- 
F = [0.874056113176056][0.3378228594999251, 0.4733472739579154]
x = [0.874056113176056][0.3378228594999251, 0.4733472739579154]
typeof(F) = ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}
typeof(x) = ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}
eltype(F) = Float64
eltype(x) = Float64
df = [NaN][NaN, NaN]
df_old = [NaN NaN NaN; NaN NaN NaN; NaN NaN NaN]
typeof(df) = ArrayPartition{Float64,Tuple{Array{Float64,1},Array{Float64,1}}}
typeof(df_old) = Array{Float64,2}
eltype(df) = Float64
eltype(df_old) = Float64
-------------------

It's not clear to me how df is used in NLSolversBase.jl, but it seems that similar is working correctly for ArrayPartition. Is this a correct assumption that df does not necessarily match the same shape as the input? I think the two broken cases are the only places where this difference pops up. It seems we can easily handle this edge case and still use similar.

pkofod · 2020-06-10T19:24:59Z

It's not clear to me how df is used in NLSolversBase.jl, but it seems that similar is working correctly for ArrayPartition.

Hm, your a is oneto3 and oneto2 but size(df_old) is (3,) and only holds 3 elements.

I have an idea to make this work

charleskawczynski · 2020-06-10T19:29:04Z

Yeah, IIUC how ArrayPartitions works, the rank changes and this complicates using similar.

pkofod · 2020-06-10T19:32:13Z

Hessian rarely plays well with these types of arrays, so I'm just going to change the hessian cache allocator. Let me run my tests and then I'll conclude :) Thanks for the help so far btw.

Change H allocator.

charleskawczynski · 2020-06-10T20:34:38Z

Alright, all tests are passing locally with the new H allocator. @pkofod feel free to merge, or let me know if you want me to squash first.

Allocate df with similar, fixes GPU support

811d4b7

longemen3000 reviewed Jun 10, 2020

View reviewed changes

pkofod and others added 2 commits June 10, 2020 21:38

Change H allocator.

2260e7b

Merge pull request #1 from JuliaNLSolvers/GPU

9aba553

Change H allocator.

pkofod merged commit e7f78c9 into JuliaNLSolvers:master Jun 10, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Allocate df with similar, partially fixes GPU support #123

Allocate df with similar, partially fixes GPU support #123

charleskawczynski commented Jun 6, 2020

longemen3000 left a comment

pkofod commented Jun 10, 2020

pkofod commented Jun 10, 2020

pkofod commented Jun 10, 2020 •

edited

Loading

charleskawczynski commented Jun 10, 2020 •

edited

Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 •

edited

Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 •

edited

Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 •

edited

Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020

Allocate df with similar, partially fixes GPU support #123

Allocate df with similar, partially fixes GPU support #123

Conversation

charleskawczynski commented Jun 6, 2020

longemen3000 left a comment

Choose a reason for hiding this comment

pkofod commented Jun 10, 2020

pkofod commented Jun 10, 2020

pkofod commented Jun 10, 2020 • edited Loading

charleskawczynski commented Jun 10, 2020 • edited Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 • edited Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 • edited Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020 • edited Loading

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020

pkofod commented Jun 10, 2020

charleskawczynski commented Jun 10, 2020

pkofod commented Jun 10, 2020 •

edited

Loading

charleskawczynski commented Jun 10, 2020 •

edited

Loading

charleskawczynski commented Jun 10, 2020 •

edited

Loading

charleskawczynski commented Jun 10, 2020 •

edited

Loading

charleskawczynski commented Jun 10, 2020 •

edited

Loading