OK starting from Greg's comment, I think I now see where I had erred. The
main problem was that, while blocks matrices are conceptually matrices,
they are not subtypes of AbstractArray{T,2}. That is to say, AbstractArrays
are not sufficiently abstract for my purposes, since I need to allow for
the possibility that Block may themselves contain blocks. To achieve this I
am now writing (foregoing my concerns about enforcing the eltype, since in
practice it will be Float64):
abstract AbstractBlockMatrix
typealias VeryAbstractMatrix Union(AbstractBlockMatrix,AbstractMatrix{
Float64})
typealias VAT VeryAbstractMatrix
type BlockMatrix{TA<:VAT,TB<:VAT,TC<:VAT,TD<:VAT} <: AbstractBlockMatrix
A::TA
B::TB
C::TC
D::TD
end
One can then construct various block matrices, with type information
attached
M1 = BlockMatrix([rand(2,2) for _=1:4]...)
M2 = BlockMatrix([Diagonal(rand(2)) for _=1:4]...)
M3 = BlockMatrix(M2,M1,M1,M2)
for M in [M1,M2,M3]
println(typeof(M),"\n\n")
end
yields the correct types. The fun can then be had
#not very efficient, but gets it done on first pass
unblocked{T<:Number}(x::T) = x
unblocked{T<:AbstractMatrix}(M::T) = M
unblocked(M::BlockMatrix) = map(unblocked,[M.A M.B; M.C M.D])
#one and one specialized method
det(M::BlockMatrix{AbstractMatrix{Float64},AbstractMatrix{Float64},
AbstractMatrix{Float64},AbstractMatrix{Float64}}) = unblocked(M)
function det(M::BlockMatrix{Diagonal{Float64},Diagonal{Float64},Diagonal{
Float64},Diagonal{Float64}})
return det(M.A*M.D-M.B*M.C)
end
abs(det(M2)-det(unblocked(M2)))
>>8.673617379884035e-19
Hope that renders correctly. Thanks for the help.
On Friday, March 13, 2015 at 3:20:45 AM UTC+1, Greg Plowman wrote:
>
> I don't really understand how this works, but this might point someone in
> the right direction.
> It seems Julia can't fully infer types, in particular the element type S.
> So we get further if we give a hint:
>
> type BlockMatrix{S,TA<:AbstractMatrix{S},TB<:AbstractMatrix{S},TC<:
> AbstractMatrix{S},TD<:AbstractMatrix{S}} <: AbstractMatrix{S}
> A::TA
> B::TB
> C::TC
> D::TD
> end
>
> typealias Block{S} BlockMatrix{S,AbstractMatrix{S},AbstractMatrix{S},
> AbstractMatrix{S},AbstractMatrix{S}}
>
> # not really sure what size() should be, but need to define for output
> Base.size(x::BlockMatrix) = (size(x.A,1) + size(x.C,1), size(x.A,2) + size
> (x.B,2))
>
> N = Block{Float64}(A,A,A,B)
>
>
> julia> N.A
> 4x4 Array{Float64,2}:
> 0.805914 0.473687 0.721984 0.464178
> 0.306 0.728015 0.148804 0.776728
> 0.439048 0.566558 0.72709 0.524761
> 0.255731 0.16528 0.331941 0.167353
> julia> N.B
> 4x4 Array{Float64,2}:
> 0.805914 0.473687 0.721984 0.464178
> 0.306 0.728015 0.148804 0.776728
> 0.439048 0.566558 0.72709 0.524761
> 0.255731 0.16528 0.331941 0.167353
> julia> N.C
> 4x4 Array{Float64,2}:
> 0.805914 0.473687 0.721984 0.464178
> 0.306 0.728015 0.148804 0.776728
> 0.439048 0.566558 0.72709 0.524761
> 0.255731 0.16528 0.331941 0.167353
> julia> N.D
> 4x4 Diagonal{Float64}:
> 1.0 0.0 0.0 0.0
> 0.0 2.0 0.0 0.0
> 0.0 0.0 3.0 0.0
> 0.0 0.0 0.0 4.0
>
>
>
>
>
> On Friday, March 13, 2015 at 8:49:41 AM UTC+11, Gabriel Mitchell wrote:
>
>> @ggggg Sorry, I guess I didn't state my intent that clearly. While your
>> example does enforce the Matrix/eltype constraint that is only part of what
>> I am after. Having a type parameter for each block is a main thing that I
>> am interested in. The reason is that I can write methods that dispatch on
>> those types. An example of such a method with an explicit 4-ary structure
>> would be
>>
>> #generic fallback
>> det(A::Matrix,B::Matrix,C::Matrix,D::Matrix) = det([A B; C D])
>>
>> #specialized method, should actually, check for A invertible, but you get
>> the idea
>> det(A::Diagonal,B::Diagonal,C::Diagonal,D::Diagonal) =
>> det(A)*det(D-C*inv(A)*B)
>>
>> #etc...
>>
>> In my applications there are at least a dozen situations where certain
>> block structure allow for significantly more efficient implementations than
>> the generic fallback. One would like to make the calls to these methods
>> (det,inv,trace, and so on) with the normal 1-ary argument, that matrix M
>> itself. This would be possible if the type information of the blocks could
>> be read out of the type of M. I hope this clear up my motivation for the
>> above question.
>>
>>
>>
>> On Thursday, March 12, 2015 at 8:49:05 PM UTC+1, ggggg wrote:
>>>
>>> BlockMatrix only needs one type parameter to fully specify the type, so
>>> you should probably only use one type parameter. Like so:
>>>
>>> *type BlockMatrix{S} <: AbstractMatrix{S}*
>>>>
>>>> * A::AbstractMatrix{S}*
>>>>
>>>> * B::AbstractMatrix{S}*
>>>>
>>>> * C::AbstractMatrix{S}*
>>>>
>>>> * D::AbstractMatrix{S}*
>>>>
>>>> *end*
>>>>
>>>>
>>> I'm sure someone else can explain in more detail why yours didn't work.
>>>
>>
