Skip to content

fast methods for sparse ± UniformScaling #631

@stevengj

Description

@stevengj

As described in this discourse thread, if A is a sparse matrix then A + x*I is hitting a specialized method, whereas A' + x*I hits a generic method that mutates one element of the diagonal at a time, which is very slow for sparse matrices.

Should be easily fixable by adding some specialized methods, e.g.

import Base: +, -
using LinearAlgebra: AdjOrTrans, wrapperop
(+)(A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}, J::UniformScaling{<:Number}) =
    wrapperop(A)(parent(A) + wrapperop(A)(J))
(-)(A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}, J::UniformScaling{<:Number}) =
    wrapperop(A)(parent(A) - wrapperop(A)(J))
(+)(J::UniformScaling{<:Number}, A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}) =
    wrapperop(A)(wrapperop(A)(J) + parent(A))
(-)(J::UniformScaling{<:Number}, A::AdjOrTrans{<:Any, <:AbstractSparseMatrix}) =
    wrapperop(A)(wrapperop(A)(J) - parent(A))

Should be an easy PR for someone?

Should be fixed by JuliaLang/LinearAlgebra.jl#1366

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions