You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: README.md
+59-59Lines changed: 59 additions & 59 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -4,71 +4,100 @@
4
4
5
5
`FastTransforms.jl` allows the user to conveniently work with orthogonal polynomials with degrees well into the millions.
6
6
7
-
Transforms include conversion between Jacobi polynomial expansions, with Chebyshev, Legendre, and ultraspherical polynomial transforms as special cases. For the signal processor, all three types of nonuniform fast Fourier transforms available. As well, spherical harmonic transforms and transforms between orthogonal polynomials on the triangle allow for the efficient simulation of partial differential equations of evolution.
7
+
This package provides a Julia wrapper for the [C library](https://github.com/MikaelSlevinsky/FastTransforms) of the same name. Additionally, all three types of nonuniform fast Fourier transforms available, as well as the Padua transform.
8
8
9
-
Algorithms include methods based on asymptotic formulae to relate the transforms to a small number of fast Fourier transforms, matrix factorizations based on the Hadamard product, hierarchical matrix decompositions à la Fast Multipole Method, and the butterfly algorithm.
9
+
## Installation
10
10
11
-
## The Chebyshev—Legendre Transform
12
-
13
-
The Chebyshev—Legendre transform allows the fast conversion of Chebyshev expansion coefficients to Legendre expansion coefficients and back.
11
+
The build script, which works on macOS and Linux systems with x86_64 processors, downloads precompiled binaries of the latest version of [FastTransforms](https://github.com/MikaelSlevinsky/FastTransforms). This library depends on `FFTW`, `MPFR`, and `OpenBLAS` (on Linux). Therefore, installation may be as straightforward as:
The 26 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:25)`. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:
The implementation separates pre-computation into a type of plan. This type is constructed with either `plan_leg2cheb` or `plan_cheb2leg`. Let's see how much faster it is if we pre-compute.
44
+
The implementation separates pre-computation into an `FTPlan`. This type is constructed with either `plan_leg2cheb` or `plan_cheb2leg`. Let's see how much faster it is if we pre-compute.
31
45
32
46
```julia
33
47
julia> p1 =plan_leg2cheb(c);
34
48
35
49
julia> p2 =plan_cheb2leg(c);
36
50
37
51
julia>@timeleg2cheb(c);
38
-
0.082615 seconds (11.94 k allocations:31.214 MiB, 6.75% gc time)
52
+
0.433938 seconds (9allocations:64.641 KiB)
39
53
40
54
julia>@time p1*c;
41
-
0.004297 seconds (6 allocations:78.422 KiB)
55
+
0.007927 seconds (79 allocations:68.563 KiB)
42
56
43
57
julia>@timecheb2leg(c);
44
-
0.110388 seconds (11.94 k allocations:31.214 MiB, 8.16% gc time)
58
+
0.423865 seconds (9allocations:64.641 KiB)
45
59
46
60
julia>@time p2*c;
47
-
0.004500 seconds (6 allocations:78.422 KiB)
61
+
0.009164 seconds (89 allocations:69.672 KiB)
62
+
63
+
```
64
+
65
+
Furthermore, for orthogonal polynomial connection problems that are degree-preserving, we should expect to be able to apply the transforms in-place:
66
+
67
+
```julia
68
+
julia>lmul!(p1, c);
69
+
70
+
julia>lmul!(p2, c);
71
+
72
+
julia>ldiv!(p1, c);
73
+
74
+
julia>ldiv!(p2, c);
75
+
48
76
```
49
77
50
-
## The Chebyshev—Jacobi Transform
78
+
###The spherical harmonic transform
51
79
52
-
The Chebyshev—Jacobi transform allows the fast conversion of Chebyshev expansion coefficients to Jacobi expansion coefficients and back.
80
+
Let `F` be an array of spherical harmonic expansion coefficients with columns arranged by increasing order in absolute value, alternating between negative and positive orders. Then `sph2fourier` converts the representation into a bivariate Fourier series, and `fourier2sph` converts it back. Once in a bivariate Fourier series on the sphere, `plan_sph_synthesis` converts the coefficients to function samples on an equiangular grid that does not sample the poles, and`plan_sph_analysis` converts them back.
53
81
54
82
```julia
55
-
julia> c =rand(10001);
83
+
julia> F =sphrandn(Float64, 1024, 2047); # convenience method
Composition of transforms allows the Jacobi—Jacobi transform, computed via `jjt`. The remainder in Hahn's asymptotic expansion is valid for the half-open square `(α,β) ∈ (-1/2,1/2]^2`. Therefore, the fast transform works best when the parameters are inside. If the parameters `(α,β)` are not exceptionally beyond the square, then increment/decrement operators are used with linear complexity (and linear conditioning) in the degree.
100
+
Due to the structure of the spherical harmonic connection problem, these transforms may also be performed in-place with `lmul!`and `ldiv!`.
72
101
73
102
## Nonuniform fast Fourier transforms
74
103
@@ -124,41 +153,12 @@ julia> N = div((n+1)*(n+2), 2);
124
153
125
154
julia> v =rand(N); # The length of v is the number of Padua points
Let `F` be a matrix of spherical harmonic expansion coefficients with columns arranged by increasing order in absolute value, alternating between negative and positive orders. Then `sph2fourier` converts the representation into a bivariate Fourier series, and `fourier2sph` converts it back.
As with other fast transforms, `plan_sph2fourier` saves effort by caching the pre-computation. Be warned that for dimensions larger than `1,000`, this is no small feat!
161
-
162
162
# References:
163
163
164
164
[1] B. Alpert and V. Rokhlin. <ahref="http://dx.doi.org/10.1137/0912009">A fast algorithm for the evaluation of Legendre expansions</a>, *SIAM J. Sci. Stat. Comput.*, **12**:158—179, 1991.
@@ -171,7 +171,7 @@ As with other fast transforms, `plan_sph2fourier` saves effort by caching the pr
171
171
172
172
[5] R. M. Slevinsky. <ahref="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, *IMA J. Numer. Anal.*, **38**:102—124, 2018.
173
173
174
-
[6] R. M. Slevinsky. <ahref="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, in press at *Appl. Comput. Harmon. Anal.*, 2017.
174
+
[6] R. M. Slevinsky. <ahref="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, *Appl. Comput. Harmon. Anal.*, **47**:585—606, 2019.
175
175
176
176
[7] R. M. Slevinsky, <ahref="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.
0 commit comments