Hi,
I'm trying to define my own custom operator type that will allow me to
implement my own * and '* operations for use inside eigenvalue or singular
value routines like eigs and svds. But I'm having trouble making this work.
Here's a simple example reimplementing matrix multiplication, with
questions sprinkled as comments in the code.
-----------
import Base: *, Ac_mul_b, size
# Ac_mul_b doesn't seem to live in base; where does it live?
type MyOperator
A::Array{Float64}
end
o = MyOperator(randn(5,4))
*(o::MyOperator, v::AbstractVector) = (println("using custom * method");
o.A*v)
o*rand(4) # this works
Ac_mul_b(o::MyOperator, v::AbstractVector) = (println("using custom '*
method"); o.A'*b)
o'*rand(5) # this doesn't work; instead, it calls ctranspose(o)*v, which is
not what I want
size(o::MyOperator) = size(o.A)
size(o::MyOperator, i::Int) = size(o)[i]
svds(o, nsv=1)
# this doesn't work; svds requires a `zero` method for the parametrized
type of my operator. If I know the type will always be Float64, what's the
easiest way to tell this to svds?
-----------
Questions, collected:
1. Ac_mul_b doesn't seem to live in base; where does it live? Or should I
be extending a different function? This causes o'*v to call
ctranspose(o)*v, which is not what I want. (The default ctranspose for new
types seems to be the identity.)
2. svds requires a `zero` method for the parametrized type of my operator.
If I know the type will always be Float64, what's the easiest way to tell
this to svds? Here's the eigs/svds
<https://github.com/JuliaLang/julia/blob/master/base/linalg/arnoldi.jl>
code, but I haven't been able to track down enough information on
AbstractMatrices and parametrized types to make this work.
3. Any other methods I should implement for my operator?
Thanks!
Madeleine
