Precompute and save data for fast1 spherical harmonic transform.
More...
Defines |
|
#define | ENTEREXIT 0 |
Functions |
| REAL | fast1_tofile (int numpts, int band_limit, REAL *node, REAL *weight, REAL epsilon, FILE *fp) |
| void | init_fast1 (Pm_1D **out, Fast1_Save *saves_fordump, int halfpts, int band_limit, REAL *node, int *numlevels, REAL epsilon, REAL(*bell_choice)(int, int)) |
| void | dump_fast1 (int numpts, int band_limit, int num_levels, REAL *node, REAL *weight, Fast1_Save *saves_fordump, Pm_1D *pmn_matrix, FILE *fp) |
Detailed Description
Precompute and save data for fast1 spherical harmonic transform.
Summary:
- fast1_tofile drives the initialization and saving.
- init_fast1 does the actual initialization.
- dump_fast1 does the actual saving.
Function Documentation
This routine dumps the initialization constructed by init_fast1 and general parameters for the fast transform for spherical harmonics, version 1.
INPUTS:
- numpts -- the number of points on the full interval.
- band_limit -- the bound on the degree of the expansion. In \(P^m_n\) notation this bounds n.
- num_levels -- the number of dyadic levels used in the adaption
- node -- an array of length numpts with the nodes, usually Arcsin of the usual gaussian nodes.
- weights -- an array of length numpts with the quadrature weights
- saves_fordump -- a structure containing the information needed for the bells. Constructed by init_fast1
- pmn_matrix -- an array of band_limit (==order_limit) Pm_1D structures, each representing the compressed matrix for one m.
OUTPUTS: fp -- the objects above are written into the file pointed to by fp, in binary format. The file must be opened before calling this routine.
This routine initializes the fast transform for spherical harmonics (version 1) and dumps it to a file. This version is designed for the novice user, and so doesn't allow tuning.
INPUTS:
- numpts -- the number of points on the full interval. Must be divisible by 4. More factors of 2 will improve performance.
- band_limit -- the bound on the degree of the expansion. In \(P^m_n\) notation this bounds n.
- node -- an array of length numpts with the nodes, usually Arcsin of the usual gaussian nodes.
- weight -- an array of length numpts with the quadrature weights
- epsilon -- the level at which to truncate. The error in the transform will be slightly larger than this
OUTPUTS:
- fp -- the objects created are written onto the file pointed to by fp, in binary format. The file must be opened before calling this routine.
- returns the compression ratio, which is an upper bound on the speed up factor of using this algorithm
NOTES:
- node and weight are input so the user can select.
- we hard-wire the largest number of levels legal. This will slow initialization but is the best thing for novice users
- The safety factor on the bell may not be quite right
- we leak a little bit of memory from the trig plans
This routine initializes for the fast transform for spherical harmonic, version 1.
Version 1 compresses in 1D (each Pmn separately) and chooses a single partition for each m and parity. By choosing a single partition we minimize the trig transform cost, but sacrifice precision. This is appropriate for small problems, near the break-even point, since it reduces overhead.
INPUTS:
- halfpts -- half the number of points on the full interval. This implies the full number of points is even, and we also assume the nodes are symmetric. Currrently halfpts must also be even for our trig transforms to work.
- band_limit -- the bound on the degree of the expansion. In \(P^m_n\) notation this bounds n. For gaussian nodes this must be < 2*halfpts and for other nodes it should be even smaller. We confirn band_limit < 2*halfpts but leave the rest to the user.
- node -- an array with the nodes, usually Arcsin of the usual gaussian nodes. We use halfpts of these, representing half of the interval.
- bell_choice -- a pointer to the bell function to use. This should be tuned to the desired precision. The choice of bell is subtle, and handled elsewhere.
- numlevels -- the number of dyadic levels to use when trying to adapt. halfpts must be divisible by 2^numlevels. If numlevels is input too large or negative, this routine will choose the maximum allowed.
- epsilon -- the target precision. The fast algorithm is an approximate method, tunable to the desired precision. The algorithm is faster the less precision (larger epsilon but < 1) you request. On the other hand, doubling the number of digits slows the algorithm by less than a factor of 2, so asking for a little extra is safe. If epsilon is chosen too close to the round-off error, the algorithm will degenerate to direct matrix applications.
OUTPUTS:
- out -- memory is allocated for band_limit (==order_limit) Pm_1D structures, each representing the compressed matrix for one m. Various bits of memory are allocated for these, and values filled in.
- saves_fordump -- bell values and total length are filled in here so we can dump them to file efficiently.
- numlevels -- is modified to the maximum allowed number if the input number was illegal.
NOTES:
- halfpts must be even for this routine to function. Performance (speed) of the algorithm will improve if halfpts has more factors of 2. For best performance, the non-power-of-two factor in halfpts should be small (10 or so).
- output ? total coefficients,bands,rows,intervals ? modified num_levels ?some holder for the trig plans
- leaks some memory
- right FFTW flags? import and export wisdom?
- min_gap hard wired
- bdmat_raw may not have enough bands
1.8.0