Proposal: moving fftpack under fortran-lang

We have discussed in the past the possibility of maintaining packages under fortran-lang at GitHub.

I would like to propose one such move: https://github.com/certik/fftpackhttps://github.com/fortran-lang/fftpack

@zoziha has done a great job modernizing it by adding fpm support, CI, modern interfaces, etc. So it would be at least he and I maintaining it. And we are hoping others would join also.

Let me know what you think.

2 Likes

Fast Fourir Transform (FFT) is a well-known and practical algorithm. fftpack is the native fft implementation of Fortran stored in netlib. I believe we have the need to maintain and use it. My original intention was to make fftpack easier-to-use when fortran-lang was rebuilding the Fortran community. We have some APIs we can refer to:

  1. keurfonluu/FFTPack: FFTPack aims to provide an easily usable package of functions using FFTPack library (Fortran 77). (github.com)
  2. Legacy discrete Fourier transforms (scipy.fftpack) — SciPy v1.7.1 Manual
  3. There is some more information here (some are outdated):
    To update the fftpack package · Issue #1 · certik/fftpack (github.com) ;
    add fftpack package by zoziha · Pull Request #40 · fortran-lang/fpm-registry (github.com).

@certik I still have 2 tasks waiting to be submitted. Please wait until I finish these 2 tasks today or tomorrow before we migrate certik/fftpack to fortran-lang/fftpack.

1 Like

Yes of course, we’ll wait.

I also have a modern Fortran reimplementation of FFTPACK here: hfsolver/fourier.f90 at master · certik/hfsolver · GitHub, but that should go into a separate package. For fortran-lang/fftpack I would like to preserve the original algorithm, for benchmarking and reference.

Yes please! I haven’t actually used the version at certik/fftpack yet, but I’ve looked at the interface (as explained in the new spec files) and I like what I see. I am definitely in favour of bringing this into fortran-lang, and while I’m rather overloaded right now, in principle I’m willing to help maintain it in the future.

1 Like

Yes, that’s exactly what I would like to do.

API-doc website: Fortran-lang/stdlib or maybe Contributing and specs – Fortran-lang/stdlib

1 Like

I think it’s just a miscommunication over exactly what is meant by ‘internal API for stdlib’ (I’m actually not clear either)

You are right, the kind parameters in stdlib are one of the two modules which are currently undocumented in the internal API documentation. It’s a low-hanging fruit for a new contributor to fix or a straight-forward task for the existing developers to take care of, if nobody takes this up, I’ll send a patch by end of the week.

I opened an issue at stdlib to continue this discussion

And there is already support to fix this short-coming in stdlib.

Ah, now I see the issue. Thank you. So if we generate the list of kinds at build time, and then use that list to generate our generic implementations of stdlib functions, how would you propose handling symbolic shorthands like dp, etc? Should we just omit those for the purpose of stdlib, and assume users of the library will define their own based on whatever kinds they need?

And just out of curiosity, the possible five real kinds you’re referring to… are are they the common 32/64/128 bit kinds, plus 16 bit and 80 bit?

Thanks @kargl for the feedback. I think what you are suggesting makes sense. Please keep the messages up for a while until we document your suggestion.

@zoziha, anything else you want to do before we move fftpack under fortran-lang?

No sir. :heart:

1 Like

Ok, I just moved it under fortran-lang! Here is the repository:

GitHub - fortran-lang/fftpack: double precision version of fftpack

2 Likes

The Fast Fourier Transform section of Fortran code on GitHub lists the following repos. Modern_fftpack has 12 stars.

2DECOMP&FFT: scalable library for Fourier transforms, derived from 2decomp&fft, by BenMql. Compatible with in-core Chebyshev transforms.

ffte: computes Discrete Fourier Transforms of 1-, 2- and 3- dimensional sequences of length (2^p) (3^q) (5^r).

fftlog-f90: extended Fortran 90 version of the FFTLog code by Andrew Hamilton to convolve a tabulated function with a Bessel functions, by coccoinomane

fft-overlap: efficient implementations of ffts on multiple GPUs and across multiple nodes, by dappelha. Overlapping data transfer on multiple levels.

fftpack: double precision version of original fftpack, from fortran-lang

FFTPack: easily usable package of functions using wrapping the Fortran 77 FFTPack library, by keurfonluu

modern_fftpack: Fortran 2008 library of fast Fourier transforms – an object-oriented modernization of NCAR’s FFTPACK5.1., by jlokimlin

nufft: code for the 1D, 2D and 3D non-uniform FFTs, from biotrump

2 Likes

@Beliavsky thanks for the list. I also wrote a modern Fortran version of fftpack here: hfsolver/fourier.f90 at master · certik/hfsolver · GitHub. One has to be careful with performance, for example allocating memory in the inner loop: modern_fftpack/complex_forward_1d.f90 at 6909d44988925dcae1ee478c06be31e5605d3974 · jlokimlin/modern_fftpack · GitHub might slow things down. For that reason, I think the best way forward is to “preserve” the original code as much as possible, for benchmarking and reference purposes, but do the following changes:

  • add modern Fortran interfaces
  • fpm, CI, … (mostly done)
  • possibly move from fixed form to free form; and move the global subroutines to modules
  • add single precision version
  • add benchmarks

But changes like rewriting the complex passes to use more modern Fortran as I have done here: hfsolver/fourier.f90 at b4c50c1979fb7e468b1852b144ba756f5a51788d · certik/hfsolver · GitHub should be done in a separate library.

Such ash fft.f. :slight_smile: / See Naming convention for Fortran wrappers packages /

Thank you both for your work on this! I’ve already downloaded it, built it with FPM, and played around with the modern interfaces. Quick question (probably for @zoziha) - the documentation says "The contents of wsave must not be changed between calls of dfftf or dfftb" - can I re-use wsave in between multiple forward-backward transforms, or do I need to call dffti before every forward transform?

1 Like

From my perceptual point of view, keurfonluu/FFTPack did a very great job. The original intention to improve fortran-lang/fftpack also stems from my desire to submit keurfonluu/FFTPack to fpm-registry (see add fftpack package by zoziha · Pull Request #40 · fortran-lang/fpm-registry (github.com))

Some codes based on fftpack5.0/5.1 may have license issues, we can keep a little conservative. (see To update the fftpack package · Issue #1 · fortran-lang/fftpack (github.com))

@Beliavsky Thank you for your sharing.

I may not understand the specific details. I just transferred the instructions in the doc file from netlib/fftpack to fftpack.md. According to the instructions in the doc file, it seems best to call again, of course you can try it boldly. :joy:

In addition, English is not my native language, and some writing in fftpack.md still needs review and improvement. (see Add `(i)rfft/dffti/dfftf/dfftb` interface and ready to move to `fortran-lang` by zoziha · Pull Request #6 · fortran-lang/fftpack · GitHub) :crazy_face:

1 Like

cc: @certik, @awvwgk , @aradi and team working on stdlib

I’ve not studied FYPP but this preprocessor appears to be employed actively in stdlib. With FYPP, is it possible somehow to make use of REAL_KINDS array named constant in ISO_FORTRAN_ENV in the code implementations of algorithms (and containers) in stdlib, rather than the current use of REAL32, REAL64, and REAL128, none of which are guaranteed to be supported by a conforming Fortran processor.

Note with REAL_KINDS constant, a Fortran standard conforming processor can be expected to have its size be at least 2, often 3 (e.g., IFORT), a few others like gfortran and NagFor to be 4, and in the near future, some processors may have it of be of size 5 or more. As you all know, the precise size of this array and the values of its elements would not matter to a truly generic algorithm. And this is something I hope will taken into close consideration during the Generics development with Fortran 202Y.

In the mean time, it will be good if there is a FYPP approach that can be used with this REAL_KINDS constant. For example, consider a (hypothetical and simple-minded) module toward a generic function to compute the statistical mean of rank-1 arrays:

module xx_stats_rank1_mean
! Module for a generic algorithm; should NOT have any dependence on REAL32, REAL64, etc.
   use, intrinsic :: iso_fortran_env, only : rk => real_kinds
   generic :: rank1_mean => rank1_mean_1, rank1_mean_2, .. rank1_mean_n ! n is the size of real_kinds
contains
   function rank1_mean_1( x ) result(mean)
      real(rk(1)), intent(in) :: x(:)
      real(rk(1)) :: mean
      mean = sum(x) / real(size(x), rk(1))
   end function rank1_mean_1
   function rank1_mean_2( x ) result(mean)
      real(rk(2)), intent(in) :: x(:)
      real(rk(2)) :: mean
      mean = sum(x) / real(size(x), rk(2))
   end function rank1_mean_2
   ..
   function rank1_mean_n( x ) result(mean)
      real(rk(n)), intent(in) :: x(:)
      real(rk(n)) :: mean
      mean = sum(x) / real(size(x), rk(n))
   end function rank1_mean_n
end module 

Can the macro language be setup using FYPP to expand out the generic function and setup the generic interface based simply on the size of REAL_KINDS constant array?

Thanks,

1 Like

Oh, yes, the current fortran-lang/fftpack is based on netlib/fftpack4.0 (public domain).
We also have a nacr/fftpack5.1 (FFTPACK license, look like MIT and BSD-3, see comment and FFTPACK5 (the “Software”) License), which natively supports multi-dimensional FFT, but gfortran compilation requires --fallow-argument-mismatch (see Build failure with gfortran 10 · Issue #123 · NCAR/ncl · GitHub) (This issue will be solved by that fpm supports native flag settings, see my discourse help).

My initial idea was to develop certik(netlib)/fftpack4.0 and nacr/fftpack5.1, first do fftpack4.0 and then fftpack5.1. I do not involve john/fftpack5.1 by default, it uses the GPL license, which may be less friendly to users.

But in fact we better maintain only one fortran-lang/fftpack?
My question is: do we need to develop based on nacr/fftpack5.1? Or let’s leave it to time, let’s finish fftpack4.0 first.

My thoughts:
First, I will try my best to improve fftpack4.0, and when it has modern features, I am willing to see whether it is necessary to maintain fftpack5.1 (may consider the needs of multi-dimensional fft), if necessary, I will hope to achieve the modern features of fftpack5.1 in another branch (I believe there is fftpack4.0 experience, this work is not difficult, it is just a matter of time). In general, the two fftpacks try to keep the same modern interface.

P.S.

  1. netlib/dfftpack1.0(fftpack4.0) (License : public domain)
    Fortran77
  2. nacr/fftpack5.1 (License: look like MIT and BSD-3, see comment)
    Fortran77
  3. John Burkardt/fftpack5.1 (License: LGPL/GPL)
    Fortran90 from nacr/fftpack5.1
  4. QcmPlab/SciFortran (License: shall be GPL)
    From John Burkardt fortran90, fftpack5.1 with interface.
  5. keurfonluu/FFTPack (License: shall be GPL)
    From John Burkardt fortran90, fftpack5.1 with interface.
  6. jlokimlin/modern_fftpack (License: FFTPACK license look like MIT and BSD-3 ?)
    Fortran90, look like from nacr/fftpack5.1 ?