Mercurial > repos > AlgTools / file diff

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/LinOps.jl	Fri Jan 08 00:26:42 2021 -0500
@@ -0,0 +1,94 @@
+########################################################
+# 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,
+       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}
+    m :: AbstractArray{T, 2}
+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
+
+end