add butterfly constructor for orthogonal matrices

Author: MikaelSlevinsky

src/SphericalHarmonics/Butterfly.jl

@@ -101,6 +101,79 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64)
     Butterfly(columns, factors, permutations, indices, zeros(T, kk), zeros(T, kk), zeros(T, kk))
 end
 
+function orthogonalButterfly{T}(A::AbstractMatrix{T}, L::Int64)
+    m, n = size(A)
+    tL = 2^L
+
+    LRAOpts = LRAOptions(T; rtol = eps(real(T))*max(m, n))
+
+    columns = Vector{Matrix{T}}(tL)
+    factors = Vector{Vector{IDPackedV{T}}}(L+1)
+    permutations = Vector{Vector{ColumnPermutation}}(L+1)
+    indices = Vector{Vector{Int64}}(L+1)
+    cs = Vector{Vector{Vector{Int64}}}(L+1)
+
+    factors[1] = Vector{IDPackedV{T}}(tL)
+    permutations[1] = Vector{ColumnPermutation}(tL)
+    indices[1] = Vector{Int64}(tL+1)
+    cs[1] = Vector{Vector{Int64}}(tL)
+
+    nd = n÷tL
+    nl = 1
+    nu = nd
+    indices[1][1] = 1
+    for j = 1:tL
+        factors[1][j] = IDPackedV{T}(collect(1:nd),Int64[],Array{T}(nd,0))
+        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)
+    for l = 2:L+1
+        factors[l] = Vector{IDPackedV{T}}(tL)
+        permutations[l] = Vector{ColumnPermutation}(tL)
+        indices[l] = Vector{Int64}(tL+1)
+        cs[l] = Vector{Vector{Int64}}(tL)
+
+        ctr = 0
+        md = m÷ii
+        ml = 1
+        mu = md
+        indices[l][1] = 1
+        for i = 1:ii
+            shft = 2jj*div(ctr,2jj)
+            for j = 1:jj
+                cols = vcat(cs[l-1][2j-1+shft],cs[l-1][2j+shft])
+                factors[l][j+ctr] = idfact(view(A, ml:mu, cols), LRAOpts)
+                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
+    for i = 1:tL
+        columns[i] = A[ml:mu,cs[L+1][i]]
+        ml += md
+        mu += md
+    end
+
+    kk = sumkmax(indices)
+
+    Butterfly(columns, factors, permutations, indices, zeros(T, kk), zeros(T, kk), zeros(T, kk))
+end
+
 function sumkmax(indices::Vector{Vector{Int64}})
     ret = 0
     @inbounds for j = 1:length(indices)

src/SphericalHarmonics/fastplan.jl

@@ -23,13 +23,13 @@ function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int)
     BF = Vector{Butterfly{T}}(n-2)
     for j = 1:2:n-2
         A_mul_B!(Ce, RP.layers[j])
-        BF[j] = Butterfly(Ce, L)
-        println("Level: ",j)
+        BF[j] = orthogonalButterfly(Ce, L)
+        #println("Level: ",j)
     end
     for j = 2:2:n-2
         A_mul_B!(Co, RP.layers[j])
-        BF[j] = Butterfly(Co, L)
-        println("Level: ",j)
+        BF[j] = orthogonalButterfly(Co, L)
+        #println("Level: ",j)
     end
     FastSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
 end

test/butterflytests.jl

@@ -77,3 +77,26 @@ for n in 7:N
     @time A_mul_B!(u, B, b)
     println(norm(u-uf)/2^n)
 end
+
+
+N = 14
+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)
+    @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)
+    uf = zero(b)
+    u = zero(b)
+    @time for k = 1:100
+        A_mul_B!(uf, A[n], b)
+    end
+    @time for k = 1:100
+        A_mul_B!(u, B[n], b)
+    end
+end

test/test_fastplan.jl

@@ -59,13 +59,25 @@ function maxcolnorm(A::AbstractMatrix)
     norm(nrm, Inf)
 end
 
-n = 1023
+function zero_spurious_modes!(A::AbstractMatrix)
+    M, N = size(A)
+    n = N÷2
+    for j = 1:n
+        for i = M-j+1:M
+            A[i,2j] = 0
+            A[i,2j+1] = 0
+        end
+    end
+    A
+end
+
+n = 255
 
 A = sphrandn(Float64, n+1, n+1)
 normalizecolumns!(A)
 Ac = deepcopy(A)
 
-FP = FastSphericalHarmonicPlan(A, 5);
+FP = FastSphericalHarmonicPlan(A, 2);
 
 B = zero(A);
 
@@ -91,20 +103,8 @@ println("The difference between fast plan and original: ", maxcolnorm(A-D))
 println("The difference between slow plan and original: ", maxcolnorm(A-E))
 println("The difference between slow and fast plans: ", maxcolnorm(D-E))
 
-function zero_spurious_modes!(A::AbstractMatrix)
-    M, N = size(A)
-    n = N÷2
-    for j = 1:n
-        for i = M-j+1:M
-            A[i,2j] = 0
-            A[i,2j+1] = 0
-        end
-    end
-    A
-end
-
-zero_spurious_modes!(D)
-zero_spurious_modes!(E)
+zero_spurious_modes!(D);
+zero_spurious_modes!(E);
 
 println("The difference between fast plan and original: ", maxcolnorm(A-D))
 println("The difference between slow plan and original: ", maxcolnorm(A-E))