modify butterfly to remove power-of-two restriction

see issue #16

Author: MikaelSlevinsky

src/SphericalHarmonics/Butterfly.jl

@@ -46,11 +46,12 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64; isorthogonal::Bool = false
     indices[1] = Vector{Int64}(tL+1)
     cs[1] = Vector{Vector{Int64}}(tL)
 
-    nd = n÷tL
-    nl = 1
-    nu = nd
+    ninds = linspace(1, n+1, tL+1)
     indices[1][1] = 1
     for j = 1:tL
+        nl = round(Int, ninds[j])
+        nu = round(Int, ninds[j+1]) - 1
+        nd = nu-nl+1
         if isorthogonal
             factors[1][j] = IDPackedV{T}(collect(1:nd),Int64[],Array{T}(nd,0))
         else
@@ -59,8 +60,6 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64; isorthogonal::Bool = false
         permutations[1][j] = factors[1][j][:P]
         indices[1][j+1] = indices[1][j] + size(factors[1][j], 1)
         cs[1][j] = factors[1][j].sk+nl-1
-        nl += nd
-        nu += nd
     end
 
     ii, jj = 2, (tL>>1)
@@ -71,41 +70,37 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64; isorthogonal::Bool = false
         cs[l] = Vector{Vector{Int64}}(tL)
 
         ctr = 0
-        md = m÷ii
-        ml = 1
-        mu = md
+        minds = linspace(1, m+1, ii+1)
         indices[l][1] = 1
         for i = 1:ii
             shft = 2jj*div(ctr,2jj)
+            ml = round(Int, minds[i])
+            mu = round(Int, minds[i+1]) - 1
             for j = 1:jj
                 cols = vcat(cs[l-1][2j-1+shft],cs[l-1][2j+shft])
                 lc = length(cols)
                 Av = A[ml:mu,cols]
                 if maximum(abs, Av) < realmin(real(T))/eps(real(T))
                     factors[l][j+ctr] = IDPackedV{T}(Int64[], collect(1:lc), Array{T}(0,lc))
                 else
-                    LRAOpts.rtol = eps(real(T))*max(length(ml:mu),lc)
+                    LRAOpts.rtol = eps(real(T))*max(mu-ml+1, lc)
                     factors[l][j+ctr] = idfact!(Av, LRAOpts)
                 end
                 permutations[l][j+ctr] = factors[l][j+ctr][:P]
                 indices[l][j+ctr+1] = indices[l][j+ctr] + size(factors[l][j+ctr], 1)
                 cs[l][j+ctr] = cols[factors[l][j+ctr].sk]
             end
-            ml += md
-            mu += md
             ctr += jj
         end
         ii <<= 1
         jj >>= 1
     end
 
-    md = m÷tL
-    ml = 1
-    mu = md
+    minds = linspace(1, m+1, tL+1)
     for i = 1:tL
-        columns[i] = A[ml:mu,cs[L+1][i]]
-        ml += md
-        mu += md
+        ml = round(Int, minds[i])
+        mu = round(Int, minds[i+1]) - 1
+        columns[i] = A[ml:mu, cs[L+1][i]]
     end
 
     kk = sumkmax(indices)

src/SphericalHarmonics/fastplan.jl

@@ -10,7 +10,6 @@ end
 
 function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
-    @assert ispow2(M)
     @assert N == 2M-1
     n = (N+1)÷2
     RP = RotationPlan(T, n-1)
@@ -36,7 +35,7 @@ function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     FastSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
 end
 
-FastSphericalHarmonicPlan(A::Matrix; opts...) = FastSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
+FastSphericalHarmonicPlan(A::Matrix; opts...) = FastSphericalHarmonicPlan(A, floor(Int, log2(size(A, 1)+1)-6); opts...)
 
 function Base.A_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1, p2, B = FP.RP, FP.BF, FP.p1, FP.p2, FP.B

src/SphericalHarmonics/thinplan.jl

@@ -14,7 +14,6 @@ end
 
 function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
-    @assert ispow2(M)
     @assert N == 2M-1
     n = (N+1)÷2
     RP = RotationPlan(T, n-1)
@@ -40,7 +39,7 @@ function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     ThinSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
 end
 
-ThinSphericalHarmonicPlan(A::Matrix; opts...) = ThinSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
+ThinSphericalHarmonicPlan(A::Matrix; opts...) = ThinSphericalHarmonicPlan(A, floor(Int, log2(size(A, 1)+1)-6); opts...)
 
 function Base.A_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1, p2, B = TP.RP, TP.BF, TP.p1, TP.p2, TP.B

test/butterflytests.jl

@@ -31,7 +31,8 @@ end
 for n in 7:N
     println("N = ", n)
     @time B = Butterfly(A[n], n-5)
-    b = rand(Complex{Float64},2^n)./(1:2^n)
+    nb = size(A[n], 2)
+    b = rand(Complex{Float64}, nb)./(1:nb)
     u = zero(b)
     @time uf = A[n]*b
     @time A_mul_B!(u, B, b)
@@ -46,20 +47,21 @@ end
 N = 10
 A = Vector{Matrix{Float64}}(N)
 for n in 1:N
-    A[n] = Float64[1/(i+j-1) for i = 1:2^n,j=1:2^n]
+    A[n] = Float64[1/(i+j-1) for i = 1:2^n+50,j=1:2^n+50]
     println(n)
 end
 
 for n in 7:N
     println("N = ", n)
     @time B = Butterfly(A[n], n-5)
-    b = rand(Float64,2^n)./(1:2^n)
+    nb = size(A[n], 2)
+    b = rand(Float64, nb)./(1:nb)
     u = zero(b)
     @time uf = A[n]*b
     @time A_mul_B!(u, B, b)
     w = zero(b)
     @time At_mul_B!(w, B, b)
-    println(norm(u-uf)/2^n)
+    println(norm(u-uf)/nb)
     println(norm(w-A[n]'b))
     println(norm(u-w))
 end
@@ -68,14 +70,15 @@ N = 10
 A = Vector{Matrix{Complex{Float64}}}(N)
 B = Vector{Butterfly{Complex{Float64}}}(N)
 for n in 7:N
-    A[n] = randnfft(2^n,2^n,0.1)
+    A[n] = randnfft(2^n+50,2^n+50,0.1)
     @time B[n] = Butterfly(A[n], n-6)
     println(n)
 end
 
 for n in 7:N
     println("N = ", n)
-    b = rand(Complex{Float64},2^n)./(1:2^n)
+    nb = size(A[n], 2)
+    b = rand(Complex{Float64}, nb)./(1:nb)
     uf = zero(b)
     u = zero(b)
     @time A_mul_B!(uf, A[n], b)

test/test_fastplan.jl

@@ -1,6 +1,6 @@
 import FastTransforms: allranks, normalizecolumns!, maxcolnorm
 
-n = 255
+n = 362
 
 A = sphrandn(Float64, n+1, n+1);
 normalizecolumns!(A);

test/test_thinplan.jl

@@ -1,6 +1,6 @@
 import FastTransforms: allranks, normalizecolumns!, maxcolnorm
 
-n = VERSION < v"0.6.0-" ? 255 : 511
+n = VERSION < v"0.6.0-" ? 362 : 724
 
 A = sphrandn(Float64, n+1, n+1);
 normalizecolumns!(A);