{"action":"closed","issue":{"url":"https://api.github.com/repos/JuliaLang/julia/issues/15630","repository_url":"https://api.github.com/repos/JuliaLang/julia","labels_url":"https://api.github.com/repos/JuliaLang/julia/issues/15630/labels{/name}","comments_url":"https://api.github.com/repos/JuliaLang/julia/issues/15630/comments","events_url":"https://api.github.com/repos/JuliaLang/julia/issues/15630/events","html_url":"https://github.com/JuliaLang/julia/issues/15630","id":143637130,"number":15630,"title":"RFC: Facilitate not paying double lookup cost when mutating a value in a LinearSlow AbstractArray","user":{"login":"KristofferC","id":1282691,"avatar_url":"https://avatars.githubusercontent.com/u/1282691?v=3","gravatar_id":"","url":"https://api.github.com/users/KristofferC","html_url":"https://github.com/KristofferC","followers_url":"https://api.github.com/users/KristofferC/followers","following_url":"https://api.github.com/users/KristofferC/following{/other_user}","gists_url":"https://api.github.com/users/KristofferC/gists{/gist_id}","starred_url":"https://api.github.com/users/KristofferC/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/KristofferC/subscriptions","organizations_url":"https://api.github.com/users/KristofferC/orgs","repos_url":"https://api.github.com/users/KristofferC/repos","events_url":"https://api.github.com/users/KristofferC/events{/privacy}","received_events_url":"https://api.github.com/users/KristofferC/received_events","type":"User","site_admin":false},"labels":[],"state":"closed","locked":false,"assignee":null,"milestone":null,"comments":5,"created_at":"2016-03-26T00:33:55Z","updated_at":"2016-03-26T17:57:25Z","closed_at":"2016-03-26T17:57:25Z","body":"There are many cases where one wants to modify a value in an array like `A[i, j] = f(A[i, j])`, If `A` is a `LinearSlow` array this would currently need two transformation between the Cartesian index and the linear index, one in `getindex` and one in `setindex!`.\r\n\r\nA typical real world example is in Finite Element Analysis where one wants to \"assemble\" a local small dense matrix `Ke` into a large sparse matrix `K`\r\n\r\n```jl\r\nfunction assemble!(K, Ke, dofs)\r\n    for (loc1, glob1) in enumerate(dofs)\r\n        for (loc2, glob2) in enumerate(dofs)\r\n            K[glob1, glob2] += Ke[loc1, loc2]\r\n        end\r\n    end\r\nend\r\n```\r\n\r\nSince `K[i,j] += Ke[a,b]` lowers to `K[i,j] = K[i,j] + K[a,b]` we will do the Cartesian to linear index searches twice. Typically in Finite Element codes to avoid this, one usually creates a `add(K, v, i::Integer, j::Integer)` method one `sub(K, v, i::Integer, j::Integer)` etc, see for example https://github.com/dealii/dealii/blob/0cb18a1aba58177c59bfe3755d13e67cc6099244/examples/step-5/step-5.cc#L391.\r\n\r\nIt could perhaps be useful for all `AbstractArrays` to optionally support a fast version of `A[i, j] = f(A[i, j])`, lets call it `mutateindex!`\r\n\r\n`mutateindex!` would have the signature `mutateindex(A, f, i...)` where `f` is a function `(v) -> f(v)` where `v` is the value of `A[i...]`\r\n\r\nA default fall back for `mutateindex` would be simple:\r\n```jl\r\nfunction mutateindex!(A, f, i...)\r\n    A[i...] = f(A[i...])\r\nend\r\n```\r\n\r\nOptionally an `AbstractArray` could implement a fast version of `mutateindex`, like for a SparseMatrix:\r\n\r\n```jl\r\nimport Base.Order.Forward\r\n\r\nfunction mutateindex!{T,Ti}(A::SparseMatrixCSC{T,Ti}, f, i0::Integer, i1::Integer)\r\n    i0 = convert(Ti, i0)\r\n    i1 = convert(Ti, i1)\r\n    if !(1 <= i0 <= A.m && 1 <= i1 <= A.n); throw(BoundsError()); end\r\n    r1 = Int(A.colptr[i1])\r\n    r2 = Int(A.colptr[i1+1]-1)\r\n    i = (r1 > r2) ? r1 : searchsortedfirst(A.rowval, i0, r1, r2, Forward)\r\n\r\n    if (i <= r2) && (A.rowval[i] == i0)\r\n        A.nzval[i] = f(A.nzval[i])\r\n    else\r\n        insert!(A.rowval, i, i0)\r\n        insert!(A.nzval, i, f(zero(T)))\r\n        @simd for j = (i1+1):(A.n+1)\r\n            @inbounds A.colptr[j] += 1\r\n        end\r\n    end\r\n    return A\r\nend\r\n```\r\n\r\nWe could now write our assembly function again as:\r\n```jl\r\nadd!(A::AbstractMatrix, v, i0::Integer, i1::Integer) = mutateindex!(A, (i) -> i + v, i0, i1)\r\n\r\nfunction assemble2!(K, Ke, dofs)\r\n    for (loc1, glob1) in enumerate(dofs)\r\n        for (loc2, glob2) in enumerate(dofs)\r\n            add!(K, Ke[loc1, loc2], glob1, glob2) \r\n        end\r\n    end\r\nend\r\n```\r\n\r\nand by doing some benchmarks:\r\n\r\n```jl\r\nfunction benchmark(K, Ke, dofs)\r\n    @time for i in 1:10^4\r\n        assemble!(K, Ke, dofs)\r\n    end\r\n\r\n    @time for i in 1:10^4\r\n        assemble2!(K, Ke, dofs)\r\n    end\r\nend\r\n\r\nconst Ke = rand(8,8)\r\nconst dofs = [20,21, 500, 501, 1001, 1002, 1086, 1087]\r\nconst K = sprand(10^5, 10^5, 0.001)\r\nassemble!(K, Ke, dofs)\r\n\r\njulia> benchmark(K, Ke, dofs)\r\n  0.033032 seconds\r\n  0.017603 seconds\r\n```\r\n\r\nwe see that we are now about twice as fast due to halving the dominating lookup cost.\r\n\r\nIt would of course be nice to have some good syntax for this, which ties on to https://github.com/JuliaLang/julia/issues/249.\r\n\r\nWe could also just have `setindex!(A, v, i, j)` be defined as `mutateindex!(A, (_) -> v, i...)  but that is probably too disruptive.\r\n\r\nMaybe our `AbstractArray` Super Saiyans @timholy and @mbauman has any comments on the feasibility of something like this.\r\n"}}