bmark_v1.jl

include("symmetric_matvec_v1.jl")
include("test.jl")

bmark_v1.log

BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     22.350 ms (0.00% GC)
  median time:      31.364 ms (0.00% GC)
  mean time:        30.879 ms (0.68% GC)
  maximum time:     39.376 ms (0.00% GC)
  --------------
  samples:          162
  evals/sample:     1
BenchmarkTools.Trial:
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     22.205 ms (0.00% GC)
  median time:      30.661 ms (0.00% GC)
  mean time:        30.149 ms (1.43% GC)
  maximum time:     37.602 ms (4.31% GC)
  --------------
  samples:          166
  evals/sample:     1
BenchmarkTools.Trial:
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     21.137 ms (0.00% GC)
  median time:      29.954 ms (0.00% GC)
  mean time:        29.176 ms (1.25% GC)
  maximum time:     35.663 ms (5.17% GC)
  --------------
  samples:          172
  evals/sample:     1

bmark_v2.jl

include("symmetric_matvec_v2.jl")
include("test.jl")

bmark_v2.log

BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     22.572 ms (0.00% GC)
  median time:      31.799 ms (0.00% GC)
  mean time:        31.534 ms (1.19% GC)
  maximum time:     43.668 ms (0.00% GC)
  --------------
  samples:          159
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     26.667 ms (0.00% GC)
  median time:      35.275 ms (0.00% GC)
  mean time:        35.609 ms (1.36% GC)
  maximum time:     47.691 ms (0.00% GC)
  --------------
  samples:          141
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     44.258 ms (0.00% GC)
  median time:      50.940 ms (0.00% GC)
  mean time:        52.069 ms (1.42% GC)
  maximum time:     66.515 ms (7.20% GC)
  --------------
  samples:          97
  evals/sample:     1

bmark_v3.jl

include("symmetric_matvec_v3.jl")
include("test.jl")

bmark_v3.log

BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     21.828 ms (0.00% GC)
  median time:      30.849 ms (0.00% GC)
  mean time:        30.339 ms (0.87% GC)
  maximum time:     38.438 ms (0.00% GC)
  --------------
  samples:          165
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     27.358 ms (0.00% GC)
  median time:      36.745 ms (0.00% GC)
  mean time:        37.298 ms (1.58% GC)
  maximum time:     52.893 ms (0.00% GC)
  --------------
  samples:          135
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     25.483 ms (0.00% GC)
  median time:      36.345 ms (0.00% GC)
  mean time:        36.801 ms (1.87% GC)
  maximum time:     50.449 ms (17.34% GC)
  --------------
  samples:          136
  evals/sample:     1

bmark_v3_joint.log

BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     23.519 ms (0.00% GC)
  median time:      33.706 ms (0.00% GC)
  mean time:        34.171 ms (0.73% GC)
  maximum time:     54.202 ms (0.00% GC)
  --------------
  samples:          147
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     26.395 ms (0.00% GC)
  median time:      35.423 ms (0.00% GC)
  mean time:        36.131 ms (1.37% GC)
  maximum time:     51.034 ms (0.00% GC)
  --------------
  samples:          139
  evals/sample:     1
BenchmarkTools.Trial: 
  memory estimate:  3.85 MiB
  allocs estimate:  2
  --------------
  minimum time:     27.737 ms (0.00% GC)
  median time:      36.531 ms (0.00% GC)
  mean time:        36.936 ms (1.51% GC)
  maximum time:     48.634 ms (0.00% GC)
  --------------
  samples:          136
  evals/sample:     1

symmetric_matvec_v1.jl

import Base: Symmetric, A_mul_B!, *

# from Base but not exported
function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

# from Base but not exported
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

function Symmetric(A::SparseMatrixCSC, uplo::Symbol=:U)
  checksquare(A)
  Symmetric{eltype(A), typeof(A)}(uplo == :U ? triu(A) : tril(A), char_uplo(uplo))
end

function (*){T <: AbstractFloat, Ti <: Integer}(A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                b::AbstractVector{T})

  x = Array{T}(size(A, 1))
  A_mul_B!(x, A, b)
end

function A_mul_B!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                     A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                     b::AbstractVector{T})

  # assume only one triangle of A is stored
  data = A.data
  length(x) ≥ data.n || throw(DimensionMismatch())
  length(b) == data.m || throw(DimensionMismatch())

  fill!(x, T(0))
  @inbounds for j = 1 : data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = data.colptr[j] : (data.colptr[j + 1] - 1)
      i = data.rowval[k]
      aij = data.nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

symmetric_matvec_v2.jl

import Base: Symmetric, A_mul_B!, *

# from Base but not exported
function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

# from Base but not exported
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

function Symmetric(A::SparseMatrixCSC, uplo::Symbol=:U)
  checksquare(A)
  Symmetric{eltype(A), typeof(A)}(A, char_uplo(uplo))  # preserve A
end

function (*){T <: AbstractFloat, Ti <: Integer}(A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                b::AbstractVector{T})

  x = Array{T}(size(A, 1))
  A_mul_B!(x, A, b)
end

function A_mul_B!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                     A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                     b::AbstractVector{T})

  # assume only one triangle of A is stored
  length(x) ≥ A.data.n || throw(DimensionMismatch())
  length(b) == A.data.m || throw(DimensionMismatch())

  fill!(x, T(0))
  return A.uplo == 'U' ? A_mul_B_U_kernel!(x, A, b) : A_mul_B_L_kernel!(x, A, b)
end

function A_mul_B_U_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = colptr[j] : (colptr[j + 1] - 1)
      i = rowval[k]
      i > j && break  # assume indices are sorted
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

function A_mul_B_L_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = (colptr[j+1] - 1) : -1 : colptr[j]
      i = rowval[k]
      i < j && break  # assume indices are sorted
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

symmetric_matvec_v3.jl

import Base: Symmetric, size, A_mul_B!, *

# from Base but not exported
function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

# from Base but not exported
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

type Symmetric{T,I<:Integer,S<:SparseMatrixCSC{T,I}} <: AbstractMatrix{T}
    data::S
    uplo::Char
    rowidx::Vector{I}
end

# annoying: must redefine all this stuff
size(A::Symmetric{T,I,S}, d) where {T,I,S} = size(A.data, d)
size(A::Symmetric{T,I,S}) where {T,I,S} = size(A.data)

function Symmetric(A::SparseMatrixCSC, uplo::Symbol=:U)
  checksquare(A)
  colptr = A.colptr
  rowval = A.rowval
  # uplo == :U : index in rowval of the last nonzero above the diagonal in each column
  # uplo == :L : index in rowval of the first nonzero below the diagonal in each column
  findfn(j) = uplo == :U ?
    findlast(i -> i ≤ j, view(rowval, colptr[j] : (colptr[j + 1] - 1))) :
    findfirst(i -> i ≥ j, view(rowval, colptr[j] : (colptr[j + 1] - 1)))

  rowidx = Array{eltype(A.colptr)}(A.n)
  for j = 1 : A.n
    i = findfn(j)  # i = 0 means nothing above/below diagonal in current column
    if uplo == :U
      rowidx[j] = i == 0 ? (colptr[j] - 1)  : (colptr[j] + i - 1)
    else
      rowidx[j] = i == 0 ? colptr[j+1]  : (colptr[j] + i - 1)
    end
  end
  Symmetric{eltype(A), eltype(A.colptr), typeof(A)}(A, char_uplo(uplo), rowidx)
end

function (*){T <: AbstractFloat, Ti <: Integer}(A::Symmetric{T, Ti, SparseMatrixCSC{T, Ti}},
                                                b::AbstractVector{T})

  x = Array{T}(size(A, 1))
  A_mul_B!(x, A, b)
end

function A_mul_B!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                     A::Symmetric{T, Ti, SparseMatrixCSC{T, Ti}},
                                                     b::AbstractVector{T})

  length(x) ≥ A.data.n || throw(DimensionMismatch())
  length(b) == A.data.m || throw(DimensionMismatch())

  fill!(x, T(0))
  # return A.uplo == 'U' ? A_mul_B_U_kernel!(x, A, b) : A_mul_B_L_kernel!(x, A, b)
  A_mul_B_UL_kernel!(x, A, b, Val{A.uplo})
end

function A_mul_B_U_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, Ti, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  rowidx = A.rowidx
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = colptr[j] : rowidx[j]
      i = rowval[k]
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

function A_mul_B_L_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, Ti, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  rowidx = A.rowidx
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = rowidx[j] : (colptr[j + 1] - 1)
      i = rowval[k]
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

function A_mul_B_UL_kernel!(x::AbstractVector{T},
                            A::Symmetric{T, Ti, SparseMatrixCSC{T, Ti}},
                            b::AbstractVector{T},
                            ::Type{Val{uplo}}) where {T,Ti,uplo}

  colptr = A.data.colptr
  rowval = A.data.rowval
  rowidx = A.rowidx
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = (uplo == :U ? (colptr[j] : rowidx[j]) :
                                    (rowidx[j] : (colptr[j + 1] - 1)))
      i = rowval[k]
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

symmetric_matvec_v4.jl

import Base: Symmetric, size, A_mul_B!, *

# from Base but not exported
function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

# from Base but not exported
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

type Symmetric{T,S<:AbstractSparseMatrix} <: AbstractMatrix{T}
    data::S
    uplo::Char
    is_triangular::Bool
end

# annoying: must redefine all this stuff
size(A::Symmetric{T,S}, d) where {T,S} = size(A.data, d)
size(A::Symmetric{T,S}) where {T,S} = size(A.data)

function Symmetric(A::SparseMatrixCSC, uplo::Symbol=:U)
  checksquare(A)
  colptr = A.colptr
  rowval = A.rowval
  # set is_triangular to true if the matrix is indeed "`uplo` triangular"
  findfn(j) = uplo == :U ?
    findfirst(i -> i > j, view(rowval, colptr[j] : (colptr[j + 1] - 1))) :
    findlast(i -> i < j, view(rowval, colptr[j] : (colptr[j + 1] - 1)))
  is_triangular = true
  for j = 1:size(A,2)
    if findfn(j) > 0
      is_triangular = false
      break
    end
  end
  Symmetric{eltype(A), typeof(A)}(A, char_uplo(uplo), is_triangular)
end

function (*){T <: AbstractFloat, Ti <: Integer}(A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                b::AbstractVector{T})

  x = Array{T}(size(A, 1))
  A_mul_B!(x, A, b)
end

function A_mul_B!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                     A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                     b::AbstractVector{T})

  length(x) ≥ A.data.n || throw(DimensionMismatch())
  length(b) == A.data.m || throw(DimensionMismatch())

  fill!(x, T(0))
  if A.is_triangular
    return A_mul_B_tri_kernel!(x, A, b)  # :U or :L doesn't matter
  else
    return A.uplo == 'U' ? A_mul_B_U_kernel!(x, A, b) : A_mul_B_L_kernel!(x, A, b)
  end
end

function A_mul_B_tri_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                                A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                                b::AbstractVector{T})

  # assume only one triangle of A is stored
  data = A.data
  length(x) ≥ data.n || throw(DimensionMismatch())
  length(b) == data.m || throw(DimensionMismatch())

  fill!(x, T(0))
  @inbounds for j = 1 : data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = data.colptr[j] : (data.colptr[j + 1] - 1)
      i = data.rowval[k]
      aij = data.nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

function A_mul_B_U_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = colptr[j] : (colptr[j + 1] - 1)
      i = rowval[k]
      i > j && break  # assume indices are sorted
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

function A_mul_B_L_kernel!{T <: AbstractFloat, Ti <: Integer}(x::AbstractVector{T},
                                                              A::Symmetric{T, SparseMatrixCSC{T, Ti}},
                                                              b::AbstractVector{T})

  colptr = A.data.colptr
  rowval = A.data.rowval
  nzval = A.data.nzval
  @inbounds for j = 1 : A.data.n
    bj = b[j]
    tmp = T(0)
    @inbounds for k = (colptr[j+1] - 1) : -1 : colptr[j]
      i = rowval[k]
      i < j && break  # assume indices are sorted
      aij = nzval[k]
      x[i] += aij * bj
      i == j || (tmp += aij * b[i])
    end
    x[j] += tmp
  end
  x
end

test.jl

using MatrixMarket
# curl -O https://www.cise.ufl.edu/research/sparse/MM/Schenk_AFE/af_shell8.tar.gz && tar zxf af_shell8.tar.gz
A = MatrixMarket.mmread("af_shell8/af_shell8.mtx")
B = Symmetric(A, :U)
b = ones(A.n)
Ab = A * b
Ab_norm = norm(Ab)
ϵ = 100 * eps(eltype(A))
@assert norm(Ab - B * b) ≤ ϵ * Ab_norm

C = Symmetric(A, :L)
@assert norm(Ab - C * b) ≤ ϵ * Ab_norm

using BenchmarkTools
Astat = @benchmark A * b
show(STDOUT, MIME"text/plain"(), Astat); println()

Bstat = @benchmark B * b
show(STDOUT, MIME"text/plain"(), Bstat); println()

Cstat = @benchmark C * b
show(STDOUT, MIME"text/plain"(), Cstat); println()

test_v4.jl

n = 10
A = sprandn(n, n, 0.5) |> t -> t + t'
b = rand(n)
Ab = A * b
Ab_norm = norm(Ab)
ϵ = 100 * eps(eltype(A))

# simulate a symmetric matrix from the upper triangle of A
B = Symmetric(A, :U)
@assert B.is_triangular == false
@assert norm(Ab - B * b) ≤ ϵ * Ab_norm

# simulate a symmetric matrix from the lower triangle of A
C = Symmetric(A, :L)
@assert C.is_triangular == false
@assert norm(Ab - C * b) ≤ ϵ * Ab_norm

# simulate a symmetric matrix from an upper triangular matrix
U = triu(A)
B = Symmetric(U, :U)
@assert B.is_triangular == true
@assert norm(Ab - B * b) ≤ ϵ * Ab_norm

# simulate a symmetric matrix from a lower triangular matrix
L = tril(A)
C = Symmetric(L, :L)
@assert C.is_triangular == true
@assert norm(Ab - C * b) ≤ ϵ * Ab_norm