Skip to content

Update syntax (no more At_mul, etc.) #74

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Sep 3, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
docs/build/
docs/site/
.DS_Store
2 changes: 1 addition & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@ os:
- linux
- osx
julia:
- 0.7
- 1.0
- 1.1
- 1.2
- 1.3
- nightly
matrix:
allow_failures:
Expand Down
32 changes: 20 additions & 12 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
name = "FastTransforms"
uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
version = "0.5"
version = "0.6"

[deps]
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
HierarchicalMatrices = "7c893195-952b-5c83-bb6e-be12f22eed51"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73"
Expand All @@ -15,13 +15,21 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"

[compat]
AbstractFFTs = "≥ 0.3.1"
Compat = "≥ 0.17.0"
DSP = "≥ 0.4.0"
FFTW = "≥ 0.0.4"
HierarchicalMatrices = "≥ 0.1.1"
LowRankApprox = "≥ 0.1.1"
ProgressMeter = "≥ 0.3.4"
SpecialFunctions = "≥ 0.3.4"
ToeplitzMatrices = "≥ 0.4.0"
julia = "0.7, 1"
AbstractFFTs = "0.4"
DSP = "0.6"
FFTW = "0.3"
HierarchicalMatrices = "0.2"
LowRankApprox = "0.2.3"
ProgressMeter = "1"
SpecialFunctions = "0.8"
ToeplitzMatrices = "0.6"
FastGaussQuadrature = "0.4"
julia = "1"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[targets]
test = ["Test","Random","Statistics"]
4 changes: 2 additions & 2 deletions appveyor.yml
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
environment:
matrix:
- julia_version: 0.7
- julia_version: 1.0
- julia_version: 1.1
- julia_version: 1.2
- julia_version: 1.2
- julia_version: 1.3
- julia_version: nightly

platform:
Expand Down
6 changes: 3 additions & 3 deletions examples/sphere.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,14 +56,14 @@ PA = FastTransforms.plan_analysis(F);

# Its spherical harmonic coefficients demonstrate that it is degree-3:
V = zero(F);
A_mul_B!(V, PA, F);
mul!(V, PA, F);
U3 = P\V

# Similarly, on the tensor product grid, the Legendre polynomial P₄(z⋅y) is:
F = [P4(z(θ,φ)⋅y) for θ in θ, φ in φ]

# Its spherical harmonic coefficients demonstrate that it is exact-degree-4:
A_mul_B!(V, PA, F);
mul!(V, PA, F);
U4 = P\V

nrm1 = vecnorm(U4);
Expand All @@ -73,7 +73,7 @@ F = [P4(z(θ,φ)⋅x) for θ in θ, φ in φ]

# It only has one nonnegligible spherical harmonic coefficient.
# Can you spot it?
A_mul_B!(V, PA, F);
mul!(V, PA, F);
U4 = P\V

# That nonnegligible coefficient should be approximately √(2π/(4+1/2)),
Expand Down
58 changes: 24 additions & 34 deletions src/FastTransforms.jl
Original file line number Diff line number Diff line change
@@ -1,45 +1,30 @@
__precompile__()

module FastTransforms

using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Compat,
AbstractFFTs, SpecialFunctions

if VERSION < v"0.7-"
using Base.FFTW
import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan, plan_r2r!,
REDFT00, REDFT01, REDFT10, REDFT11,
RODFT00, RODFT01, RODFT10, RODFT11
const LAmul! = Base.A_mul_B!
import Base: Factorization
rmul!(A::AbstractArray, c::Number) = scale!(A,c)
lmul!(c::Number, A::AbstractArray) = scale!(c,A)
lmul!(A::AbstractArray, B::AbstractArray) = mul!(A,B)
rmul!(A::AbstractArray, B::AbstractArray) = mul!(A,B)
const floatmin = realmin
else
using FFTW, LinearAlgebra, DSP
import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
REDFT00, REDFT01, REDFT10, REDFT11,
RODFT00, RODFT01, RODFT10, RODFT11
const LAmul! = LinearAlgebra.mul!
const libfftw = libfftw3
const libfftwf = libfftw3f
import LinearAlgebra: Factorization
flipdim(A,d) = reverse(A; dims=d)
end
using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter,
AbstractFFTs, SpecialFunctions, FastGaussQuadrature

using FFTW, LinearAlgebra, DSP
import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
REDFT00, REDFT01, REDFT10, REDFT11,
RODFT00, RODFT01, RODFT10, RODFT11

const libfftw = libfftw3
const libfftwf = libfftw3f
import LinearAlgebra: Factorization
flipdim(A,d) = reverse(A; dims=d)

import Base: *, \, inv, size, view
import Base: getindex, setindex!, length
import Compat.LinearAlgebra: BlasFloat, BlasInt
import Base: getindex, setindex!, length, axes
import LinearAlgebra: BlasFloat, BlasInt, transpose, adjoint, lmul!, rmul!
import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!
import HierarchicalMatrices: mul!, At_mul_B!, Ac_mul_B!
import HierarchicalMatrices: mul!
import HierarchicalMatrices: ThreadSafeVector, threadsafezeros
import LowRankApprox: ColPerm
import LowRankApprox: ColPerm, RowPerm
import AbstractFFTs: Plan
import Compat: range, transpose, adjoint, axes

import FastGaussQuadrature: unweightedgausshermite

export cjt, icjt, jjt, plan_cjt, plan_icjt
export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac, plan_jac2jac
Expand All @@ -64,6 +49,8 @@ export SlowTriangularHarmonicPlan
export tri2cheb, cheb2tri, plan_tri2cheb
export triones, trizeros, trirand, trirandn, trievaluate

export hermitepoints, weightedhermitetransform, iweightedhermitetransform

# Other module methods and constants:
#export ChebyshevJacobiPlan, jac2cheb, cheb2jac
#export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
Expand All @@ -73,6 +60,8 @@ export triones, trizeros, trirand, trirandn, trievaluate
#export plan_fejer2, fejernodes2, fejerweights2
#export RecurrencePlan, forward_recurrence!, backward_recurrence

lgamma(x) = logabsgamma(x)[1]

include("stepthreading.jl")
include("fftBigFloat.jl")
include("specialfunctions.jl")
Expand All @@ -97,6 +86,7 @@ include("inufft.jl")
include("cjt.jl")

include("toeplitzhankel.jl")
include("hermite.jl")

#leg2cheb(x...)=th_leg2cheb(x...)
#cheb2leg(x...)=th_cheb2leg(x...)
Expand Down
64 changes: 29 additions & 35 deletions src/SphericalHarmonics/Butterfly.jl
Original file line number Diff line number Diff line change
Expand Up @@ -110,10 +110,9 @@ function Butterfly(A::AbstractMatrix{T}, L::Int; isorthogonal::Bool = false, opt
Butterfly(columns, factors, permutations, indices, threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk))
end

if VERSION ≥ v"0.7-"
LinearAlgebra.adjoint(B::Butterfly) = Adjoint(B)
LinearAlgebra.transpose(B::Butterfly) = Transpose(B)
end
LinearAlgebra.adjoint(B::Butterfly) = Adjoint(B)
LinearAlgebra.transpose(B::Butterfly) = Transpose(B)


function sumkmax(indices::Vector{Vector{Int}})
ret = 0
Expand Down Expand Up @@ -173,66 +172,61 @@ function rowperm!(fwd::Bool, y::AbstractVector, x::AbstractVector, p::Vector{Int
end

## ColumnPermutation
mul!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(false, B, A.p, jstart)
At_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(true, B, A.p, jstart)
Ac_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = At_mul_B!(A, B, jstart)
lmul!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(false, B, A.p, jstart)
lmul!(At::RowPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(true, B, transpose(At).p, jstart)

mul!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(false, y, x, A.p, jstart)
At_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(true, y, x, A.p, jstart)
Ac_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = At_mul_B!(y, x, A, jstart)
mul!(y::AbstractVector, At::RowPerm, x::AbstractVector, jstart::Int) = rowperm!(true, y, x, transpose(At).p, jstart)

# Fast mul!, At_mul_B!, and Ac_mul_B! for an ID. These overwrite the output.
# Fast mul! for an ID. These overwrite the output.


function mul!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
k, n = size(A)
At_mul_B!(P, x, jstart)
lmul!(transpose(P), x, jstart)
copyto!(y, istart, x, jstart, k)
mul!(y, A.T, x, istart, jstart+k)
mul!(P, x, jstart)
lmul!(P, x, jstart)
y
end

function mul!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
k, n = size(A)
At_mul_B!(temp, P, x, jstart)
mul!(temp, transpose(P), x, jstart)
copyto!(y, istart, temp, jstart, k)
mul!(y, A.T, temp, istart, jstart+k)
y
end

### mul!, At_mul_B!, and Ac_mul_B! for a Butterfly factorization.
### mul! for a Butterfly factorization.
mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)

for f! in (:At_mul_B!, :Ac_mul_B!)
for Trans in (:Transpose, :Adjoint)
@eval begin
function $f!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
function mul!(y::AbstractVecOrMat{T}, Ac::$Trans{T,IDPackedV{T}}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
A = parent(Ac)
k, n = size(A)
copyto!(y, istart, x, jstart, k)
$f!(y, A.T, x, istart+k, jstart)
mul!(P, y, istart)
mul!(y, $Trans(A.T), x, istart+k, jstart)
lmul!(P, y, istart)
y
end

function $f!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
function mul!(y::AbstractVector{T}, Ac::$Trans{T,IDPackedV{T}}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
A = parent(Ac)
k, n = size(A)
copyto!(temp, istart, x, jstart, k)
$f!(temp, A.T, x, istart+k, jstart)
mul!(temp, $Trans(A.T), x, istart+k, jstart)
mul!(y, P, temp, istart)
y
end
end
end

if VERSION < v"0.7-"
Base.A_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
Base.At_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = At_mul_B_col_J!(u, B, b, 1)
Base.Ac_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, B, b, 1)
else
LinearAlgebra.mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
LinearAlgebra.mul!(u::Vector{T}, Bt::Transpose{T,Butterfly{T}}, b::Vector{T}) where T = At_mul_B_col_J!(u, parent(Bt), b, 1)
LinearAlgebra.mul!(u::Vector{T}, Bc::Adjoint{T,Butterfly{T}}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, parent(Bc), b, 1)
end
LinearAlgebra.mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
LinearAlgebra.mul!(u::Vector{T}, Bt::Transpose{T,Butterfly{T}}, b::Vector{T}) where T = mul_col_J!(u, Bt, b, 1)
LinearAlgebra.mul!(u::Vector{T}, Bc::Adjoint{T,Butterfly{T}}, b::Vector{T}) where T = mul_col_J!(u, Bc, b, 1)



function mul_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
Expand Down Expand Up @@ -290,10 +284,10 @@ function mul_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) whe
u
end

for f! in (:At_mul_B!,:Ac_mul_B!)
f_col_J! = Meta.parse(string(f!)[1:end-1]*"_col_J!")
for Trans in (:Transpose,:Adjoint)
@eval begin
function $f_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
function mul_col_J!(u::VecOrMat{T}, Bc::$Trans{T,Butterfly{T}}, b::VecOrMat{T}, J::Int) where T
B = parent(Bc)
L = length(B.factors) - 1
tL = 2^L

Expand All @@ -315,7 +309,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
for i = 1:tL
ml = mu+1
mu += size(columns[i], 1)
$f!(temp1, columns[i], b, inds[i], ml+COLSHIFT)
mul!(temp1, $Trans(columns[i]), b, inds[i], ml+COLSHIFT)
end

ii, jj = tL, 1
Expand All @@ -329,7 +323,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
shft = 2jj*div(ctr,2jj)
fill!(temp4, zero(T))
for j = 1:jj
$f!(temp3, factors[j+ctr], permutations[j+ctr], temp1, temp4, indsout[2j+shft-1], indsin[j+ctr])
mul!(temp3, $Trans(factors[j+ctr]), permutations[j+ctr], temp1, temp4, indsout[2j+shft-1], indsin[j+ctr])
addtemp3totemp2!(temp2, temp3, indsout[2j+shft-1], indsout[2j+shft+1]-1)
end
ctr += jj
Expand All @@ -346,7 +340,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
for j = 1:tL
nl = nu+1
nu += size(factors[j], 2)
$f!(u, factors[j], permutations[j], temp1, nl+COLSHIFT, inds[j])
mul!(u, $Trans(factors[j]), permutations[j], temp1, nl+COLSHIFT, inds[j])
end

u
Expand Down
11 changes: 3 additions & 8 deletions src/SphericalHarmonics/SphericalHarmonics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,9 @@ function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
mul!(zero(X), P, X)
end

if VERSION < v"0.7-"
function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
At_mul_B!(zero(X), P, X)
end
else
function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
mul!(zero(X), transpose(P), X)
end

function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
mul!(zero(X), transpose(P), X)
end

include("sphfunctions.jl")
Expand Down
Loading