comparison: src/LinOps.jl
src/LinOps.jl
- changeset 20
- 8b80aa64adec
- child 23
- 4399bf266660
equal
deleted
inserted
replaced
| 19:63849571a046 | 20:8b80aa64adec |
|---|---|
| 1 ######################################################## | |
| 2 # Abstract Hilbert space linear operators. Does not set | |
| 3 # any limits on the input and output data types unlike | |
| 4 # the package LinearOperators. | |
| 5 ######################################################## | |
| 6 | |
| 7 __precompile__() | |
| 8 | |
| 9 module LinOps | |
| 10 | |
| 11 export LinOp, IdOp, AdjointOp, AdjointableOp, | |
| 12 calc_adjoint, calc_adjoint!, inplace!, | |
| 13 opnorm_estimate | |
| 14 | |
| 15 # General type | |
| 16 | |
| 17 abstract type LinOp{X, Y} end | |
| 18 | |
| 19 function inplace!(y::Y, op::LinOp{X,Y}, x::X) where {X, Y} | |
| 20 y .= op(x) | |
| 21 end | |
| 22 | |
| 23 # Adjoint operator for operators that define calc_adjoint and calc_adjoint! | |
| 24 | |
| 25 abstract type AdjointableOp{X, Y} <: LinOp{X,Y} end | |
| 26 | |
| 27 struct AdjointOp{X,Y,T<:AdjointableOp{X,Y}} <: LinOp{Y,X} | |
| 28 op :: T | |
| 29 end | |
| 30 | |
| 31 function Base.adjoint(op::T) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 32 return AdjointOp{X,Y,T}(op) | |
| 33 end | |
| 34 | |
| 35 function Base.adjoint(adj::AdjointOp{X,Y,T}) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 36 return adj.op | |
| 37 end | |
| 38 | |
| 39 function (adj::AdjointOp{X,Y,T})(y::Y) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 40 calc_adjoint(adj, y) | |
| 41 end | |
| 42 | |
| 43 function inplace!(x::X, adj::AdjointOp{X,Y,T}, y::Y) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 44 calc_adjoint!(x, adj.op, y) | |
| 45 end | |
| 46 | |
| 47 function calc_adjoint(op::T, y::Y) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 48 @error "unimplemented" | |
| 49 end | |
| 50 | |
| 51 function calc_adjoint!(x::X, adj::T, y::Y) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 52 x .= calc_adjoint(adj, y) | |
| 53 end | |
| 54 | |
| 55 function opnorm_estimate(adj::AdjointOp{X,Y,T}) where {X, Y, T<:AdjointableOp{X,Y}} | |
| 56 return opnorm_estimate(adj.op) | |
| 57 end | |
| 58 | |
| 59 # Identity operator | |
| 60 | |
| 61 struct IdOp{X} <: LinOp{X,X} | |
| 62 end | |
| 63 | |
| 64 function (op::IdOp{X})(x::X) where X | |
| 65 return x | |
| 66 end | |
| 67 | |
| 68 function Base.adjoint(op::IdOp{X}) where X | |
| 69 return op | |
| 70 end | |
| 71 | |
| 72 function opnorm_estimate(op::IdOp{X}) where X | |
| 73 return 1 | |
| 74 end | |
| 75 | |
| 76 # Matrix operator | |
| 77 | |
| 78 struct MatrixOp{T<:Real} | |
| 79 m :: AbstractArray{T, 2} | |
| 80 end | |
| 81 | |
| 82 function (op::MatrixOp{T})(v::Vector{T}) where T | |
| 83 return op.m*v | |
| 84 end | |
| 85 | |
| 86 function Base.adjoint(op::MatrixOp{T}) where T | |
| 87 return MatrixOp(op.m') | |
| 88 end | |
| 89 | |
| 90 function opnorm_estimate(op::MatrixOp{T}) where T | |
| 91 return opnorm(op.m, 2) | |
| 92 end | |
| 93 | |
| 94 end |
