Variable kind

[hsnyder]

This PR allows the user to change the real/complex kind at compile time by changing the definition of RK in rk.f

My original plan (the plan I discussed with @certik) was to change IMPLICIT DOUBLE PRECISION to IMPLICIT REAL(RK) and to include a file directly about that line, which defined INTEGER, PARAMETER :: RK = KIND(1.0d0), but this didn't work because in that solution I ended up with a variable declaration before an IMPLICIT statement, which is illegal. So instead, I opted to create a module called fftpack_kind, which defines RK, and USE that module in each of the old F77 subroutines. Basically the same solution, but I use a module rather than including a loose file. I hope that's okay.

Most of these changes were done via perl -pi -e, rather than doing them manually.

[hsnyder]

I should mention that so far the only testing I have done is running fpm test --flag "-O0" (the flag is just because the code fails to build with the default flags that fpm passes to gfortran).

I'll do some numerical tests tomorrow on a few examples. Are there any more specific tests that @zoziha or @certik would like?

[zoziha]

I think what you done is the best way (use fftpack_kind), thanks @hsnyder !
Can be submitted later (see #14 (comment)), to extend different real kinds, using fpm test --flag "-g -Dsp":

module fftpack_kind
    use iso_fortran_env, only: sp => real32, dp => real64, qp => real128
#if   defined(sp)
    integer, parameter :: rk = sp
#elif defined(qp)
    integer, parameter :: rk = qp
#else
    integer, parameter :: rk = dp
#endif
end module fftpack_kind

the flag is just because the code fails to build with the default flags that fpm passes to gfortran.

Yes, because of the error reported by -fcheck=bounds.
Later, when fpm supports this feature (fortran-lang/fpm#498), this will not be a problem.
Of course, we can also change (*) and (1) in the original fftpack code to appropriate (:). (see fortran-fans@b40eb69)

I tried to implement multi-dimensional fft early on, such as fft2(see add_fft2), so fftpack changed (*) and (1) to (:) is appropriate and reasonable.

numerical tests

Regarding numerical testing, I don't have any ideas. I don't actually have enough numerical experience.

[jacobwilliams]

Note that lines like this:

DATA TSQRT2 /2.82842712474619009760D0/

are not going to be correct for quad precision. They need to be replaced with something like this:

real(rk),parameter :: tsqrt2 = 2.0_rk * sqrt(2.0_rk) 

But... I would recommend converting all these .f files to .f90 before making changes like this...

[jacobwilliams]

Also the calls to FLOAT(..) need to be changed to real(..,rk)

[hsnyder]

Good points. I fixed the FLOAT calls, but maybe we need to close this PR and think this through more carefully. Your point about the DATA statements and quad precision is well taken, and there are still lots of double precision literals scattered throughout the code (e.g. 2.0D0). There are also calls to "DCMPLX" which I don't even think is standard Fortran.

[jacobwilliams]

how about this: #18

[certik]

Closing in favor of #18.