la vie est belle sans SubArray

Author: MikaelSlevinsky

src/SphericalHarmonics/Butterfly.jl

@@ -28,11 +28,12 @@ end
 
 size(B::Butterfly) = size(B, 1), size(B, 2)
 
-function Butterfly{T}(A::AbstractMatrix{T}, L::Int64)
+function Butterfly{T}(A::AbstractMatrix{T}, L::Int64; opts...)
     m, n = size(A)
     tL = 2^L
 
-    LRAOpts = LRAOptions(T; rtol = eps(real(T))*max(m, n))
+    LRAOpts = LRAOptions(T; opts...)
+    LRAOpts.rtol = eps(real(T))*max(m, n)
 
     columns = Vector{Matrix{T}}(tL)
     factors = Vector{Vector{IDPackedV{T}}}(L+1)
@@ -50,7 +51,7 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64)
     nu = nd
     indices[1][1] = 1
     for j = 1:tL
-        factors[1][j] = idfact(view(A,:,nl:nu), LRAOpts)
+        factors[1][j] = idfact(A[:,nl:nu], LRAOpts)
         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
@@ -74,7 +75,14 @@ function Butterfly{T}(A::AbstractMatrix{T}, L::Int64)
             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)
+                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)
+                    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]
@@ -101,11 +109,12 @@ 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)
+function orthogonalButterfly{T}(A::AbstractMatrix{T}, L::Int64; opts...)
     m, n = size(A)
     tL = 2^L
 
-    LRAOpts = LRAOptions(T; rtol = eps(real(T))*max(m, n))
+    LRAOpts = LRAOptions(T; opts...)
+    LRAOpts.rtol = eps(real(T))*max(m, n)
 
     columns = Vector{Matrix{T}}(tL)
     factors = Vector{Vector{IDPackedV{T}}}(L+1)
@@ -147,7 +156,14 @@ function orthogonalButterfly{T}(A::AbstractMatrix{T}, L::Int64)
             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)
+                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)
+                    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]

src/SphericalHarmonics/fastplan.jl

@@ -8,7 +8,7 @@ immutable FastSphericalHarmonicPlan{T}
     B::Matrix{T}
 end
 
-function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int)
+function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
     @assert ispow2(M)
     @assert N == 2M-1
@@ -25,18 +25,18 @@ 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] = orthogonalButterfly(Ce, L)
+        BF[j] = orthogonalButterfly(Ce, L; opts...)
         println("Layer: ",j)
     end
     for j = 2:2:n-2
         A_mul_B!(Co, RP.layers[j])
-        BF[j] = orthogonalButterfly(Co, L)
+        BF[j] = orthogonalButterfly(Co, L; opts...)
         println("Layer: ",j)
     end
     FastSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
 end
 
-FastSphericalHarmonicPlan(A::Matrix) = FastSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6))
+FastSphericalHarmonicPlan(A::Matrix; opts...) = FastSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
 
 function 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

@@ -12,7 +12,7 @@ immutable ThinSphericalHarmonicPlan{T}
     B::Matrix{T}
 end
 
-function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int)
+function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
     M, N = size(A)
     @assert ispow2(M)
     @assert N == 2M-1
@@ -29,16 +29,16 @@ function ThinSphericalHarmonicPlan{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])
-        checklayer(j+1) && (BF[j] = orthogonalButterfly(Ce, L);println("Layer: ",j))
+        checklayer(j+1) && (BF[j] = orthogonalButterfly(Ce, L; opts...);println("Layer: ",j))
     end
     for j = 2:2:n-2
         A_mul_B!(Co, RP.layers[j])
-        checklayer(j) && (BF[j] = orthogonalButterfly(Co, L);println("Layer: ",j))
+        checklayer(j) && (BF[j] = orthogonalButterfly(Co, L; opts...);println("Layer: ",j))
     end
     ThinSphericalHarmonicPlan(RP, BF, p1, p2, p1inv, p2inv, B)
 end
 
-ThinSphericalHarmonicPlan(A::Matrix) = ThinSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6))
+ThinSphericalHarmonicPlan(A::Matrix; opts...) = ThinSphericalHarmonicPlan(A, round(Int, log2(size(A, 1)+1)-6); opts...)
 
 function A_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
     RP, BF, p1, p2, B = TP.RP, TP.BF, TP.p1, TP.p2, TP.B