src/LinOps.jl

########################################################
# Abstract Hilbert space linear operators. Does not set
# any limits on the input and output data types unlike
# the package LinearOperators.
########################################################

__precompile__()

module LinOps

export LinOp, IdOp, AdjointOp, AdjointableOp, MatrixOp,
       calc_adjoint, calc_adjoint!, inplace!,
       opnorm_estimate

# General type

abstract type LinOp{X, Y} end

function inplace!(y::Y, op::LinOp{X,Y}, x::X) where {X, Y}
    y .= op(x)
end

# Adjoint operator for operators that define calc_adjoint and calc_adjoint!

abstract type AdjointableOp{X, Y} <: LinOp{X,Y} end

struct AdjointOp{X,Y,T<:AdjointableOp{X,Y}} <: LinOp{Y,X}
    op :: T
end

function Base.adjoint(op::T) where {X, Y, T<:AdjointableOp{X,Y}}
    return AdjointOp{X,Y,T}(op)
end

function Base.adjoint(adj::AdjointOp{X,Y,T}) where {X, Y, T<:AdjointableOp{X,Y}}
    return adj.op
end

function (adj::AdjointOp{X,Y,T})(y::Y) where {X, Y, T<:AdjointableOp{X,Y}}
    calc_adjoint(adj, y)
end

function inplace!(x::X, adj::AdjointOp{X,Y,T}, y::Y) where {X, Y, T<:AdjointableOp{X,Y}}
    calc_adjoint!(x, adj.op, y)
end

function calc_adjoint(op::T, y::Y) where {X, Y, T<:AdjointableOp{X,Y}}
    @error "unimplemented"
end

function calc_adjoint!(x::X, adj::T, y::Y) where {X, Y, T<:AdjointableOp{X,Y}}
    x .= calc_adjoint(adj, y)
end

function opnorm_estimate(adj::AdjointOp{X,Y,T}) where {X, Y, T<:AdjointableOp{X,Y}}
    return opnorm_estimate(adj.op)
end

# Identity operator

struct IdOp{X} <: LinOp{X,X}
end

function (op::IdOp{X})(x::X) where X
    return x
end

function Base.adjoint(op::IdOp{X}) where X
    return op
end

function opnorm_estimate(op::IdOp{X}) where X
    return 1
end

# Matrix operator

struct MatrixOp{T<:Real} <: LinOp{Vector{T}, Vector{T}}
    m :: AbstractMatrix{T}
end

function (op::MatrixOp{T})(v::Vector{T}) where T
    return op.m*v
end

function Base.adjoint(op::MatrixOp{T}) where T
    return MatrixOp(op.m')
end

function opnorm_estimate(op::MatrixOp{T}) where T
    return opnorm(op.m, 2)
end

function Base.:+(a::MatrixOp{T}, b::MatrixOp{T}) where T
    return MatrixOp(a.m+b.m)
end

function Base.:-(a::MatrixOp{T}, b::MatrixOp{T}) where T
    return MatrixOp(a.m-b.m)
end

function Base.:*(a::MatrixOp{T}, b::MatrixOp{T}) where T
    return MatrixOp(a.m*b.m)
end

end