# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1610083602 18000
# Node ID 8b80aa64adec9dc9b6b90597b5ea9a76182e9497
# Parent  63849571a046b03e894829ae85cbf8835236aa30
Add abstract LinOps

diff -r 63849571a046 -r 8b80aa64adec README.md
--- a/README.md	Fri May 08 14:47:33 2020 -0500
+++ b/README.md	Fri Jan 08 00:26:42 2021 -0500
@@ -12,6 +12,7 @@
   * Calculation of norms, dot products, projections
   * Conditional thread execution macros
   * Template for conveniently writing iterative algorithms using `do`-block anonymous functions as algorithm steps, with visualisation/verbosity implemented generally in iterator.
+  * Abstract linear operators with arbitrary domains and codomains; `LinOp{X,Y}`.
 
   The code is used by [ImageTools][] and both, for example, by <http://dx.doi.org/10.5281/zenodo.3659180>.
 
diff -r 63849571a046 -r 8b80aa64adec src/AlgTools.jl
--- a/src/AlgTools.jl	Fri May 08 14:47:33 2020 -0500
+++ b/src/AlgTools.jl	Fri Jan 08 00:26:42 2021 -0500
@@ -8,5 +8,6 @@
 include("Iterate.jl")
 include("Util.jl")
 include("Comms.jl")
+include("LinOps.jl")
 
 end
diff -r 63849571a046 -r 8b80aa64adec src/LinOps.jl
--- /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