[1850] | 1 | !> @file singleton_mod.f90 |
---|
[2000] | 2 | !------------------------------------------------------------------------------! |
---|
[484] | 3 | ! Current revisions: |
---|
[1] | 4 | ! ----------------- |
---|
[1683] | 5 | ! |
---|
[2001] | 6 | ! |
---|
[1321] | 7 | ! Former revisions: |
---|
| 8 | ! ----------------- |
---|
| 9 | ! $Id: singleton_mod.f90 4180 2019-08-21 14:37:54Z scharf $ |
---|
[3761] | 10 | ! statement added to prevent compiler warning about unused variables |
---|
| 11 | ! |
---|
[1321] | 12 | ! |
---|
[1] | 13 | ! Description: |
---|
| 14 | ! ------------ |
---|
[1682] | 15 | !> Multivariate Fast Fourier Transform |
---|
| 16 | !> |
---|
| 17 | !> Fortran 90 Implementation of Singleton's mixed-radix algorithm, |
---|
| 18 | !> RC Singleton, Stanford Research Institute, Sept. 1968. |
---|
| 19 | !> |
---|
| 20 | !> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and |
---|
| 21 | !> John Beale. |
---|
| 22 | !> |
---|
| 23 | !> Fourier transforms can be computed either in place, using assumed size |
---|
| 24 | !> arguments, or by generic function, using assumed shape arguments. |
---|
| 25 | !> |
---|
| 26 | !> |
---|
| 27 | !> Public: |
---|
| 28 | !> |
---|
| 29 | !> fftkind kind parameter of complex arguments |
---|
| 30 | !> and function results. |
---|
| 31 | !> |
---|
| 32 | !> fft(array, dim, inv, stat) generic transform function |
---|
| 33 | !> COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array |
---|
| 34 | !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 35 | !> LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
| 36 | !> INTEGER, INTENT(OUT), OPTIONAL:: stat |
---|
| 37 | !> |
---|
| 38 | !> fftn(array, shape, dim, inv, stat) in place transform subroutine |
---|
| 39 | !> COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array |
---|
| 40 | !> INTEGER, DIMENSION(:), INTENT(IN) :: shape |
---|
| 41 | !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 42 | !> LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
| 43 | !> INTEGER, INTENT(OUT), OPTIONAL:: stat |
---|
| 44 | !> |
---|
| 45 | !> |
---|
| 46 | !> Formal Parameters: |
---|
| 47 | !> |
---|
| 48 | !> array The complex array to be transformed. array can be of arbitrary |
---|
| 49 | !> rank (i.e. up to seven). |
---|
| 50 | !> |
---|
| 51 | !> shape With subroutine fftn, the shape of the array to be transformed |
---|
| 52 | !> has to be passed separately, since fftradix - the internal trans- |
---|
| 53 | !> formation routine - will treat array always as one dimensional. |
---|
| 54 | !> The product of elements in shape must be the number of |
---|
| 55 | !> elements in array. |
---|
| 56 | !> Although passing array with assumed shape would have been nicer, |
---|
| 57 | !> I prefered assumed size in order to prevent the compiler from |
---|
| 58 | !> using a copy-in-copy-out mechanism. That would generally be |
---|
| 59 | !> necessary with fftn passing array to fftradix and with fftn |
---|
| 60 | !> being prepared for accepting non consecutive array sections. |
---|
| 61 | !> Using assumed size, it's up to the user to pass an array argu- |
---|
| 62 | !> ment, that can be addressed as continous one dimensional array |
---|
| 63 | !> without copying. Otherwise, transformation will not really be |
---|
| 64 | !> performed in place. |
---|
| 65 | !> On the other hand, since the rank of array and the size of |
---|
| 66 | !> shape needn't match, fftn is appropriate for handling more than |
---|
| 67 | !> seven dimensions. |
---|
| 68 | !> As far as function fft is concerned all this doesn't matter, |
---|
| 69 | !> because the argument will be copied anyway. Thus no extra |
---|
| 70 | !> shape argument is needed for fft. |
---|
| 71 | !> |
---|
| 72 | !> Optional Parameters: |
---|
| 73 | !> |
---|
| 74 | !> dim One dimensional integer array, containing the dimensions to be |
---|
| 75 | !> transformed. Default is (/1,...,N/) with N being the rank of |
---|
| 76 | !> array, i.e. complete transform. dim can restrict transformation |
---|
| 77 | !> to a subset of available dimensions. Its size must not exceed the |
---|
| 78 | !> rank of array or the size of shape respectivly. |
---|
| 79 | !> |
---|
| 80 | !> inv If .true., inverse transformation will be performed. Default is |
---|
| 81 | !> .false., i.e. forward transformation. |
---|
| 82 | !> |
---|
| 83 | !> stat If present, a system dependent nonzero status value will be |
---|
| 84 | !> returned in stat, if allocation of temporary storage failed. |
---|
| 85 | !> |
---|
| 86 | !> |
---|
| 87 | !> Scaling: |
---|
| 88 | !> |
---|
| 89 | !> Transformation results will always be scaled by the square root of the |
---|
| 90 | !> product of sizes of each dimension in dim. (See examples below) |
---|
| 91 | !> |
---|
| 92 | !> |
---|
| 93 | !> Examples: |
---|
| 94 | !> |
---|
| 95 | !> Let A be a L*M*N three dimensional complex array. Then |
---|
| 96 | !> |
---|
| 97 | !> result = fft(A) |
---|
| 98 | !> |
---|
| 99 | !> will produce a three dimensional transform, scaled by sqrt(L*M*N), while |
---|
| 100 | !> |
---|
| 101 | !> call fftn(A, SHAPE(A)) |
---|
| 102 | !> |
---|
| 103 | !> will do the same in place. |
---|
| 104 | !> |
---|
| 105 | !> result = fft(A, dim=(/1,3/)) |
---|
| 106 | !> |
---|
| 107 | !> will transform with respect to the first and the third dimension, scaled |
---|
| 108 | !> by sqrt(L*N). |
---|
| 109 | !> |
---|
| 110 | !> result = fft(fft(A), inv=.true.) |
---|
| 111 | !> |
---|
| 112 | !> should (approximately) reproduce A. |
---|
| 113 | !> With B having the same shape as A |
---|
| 114 | !> |
---|
| 115 | !> result = fft(fft(A) * CONJG(fft(B)), inv=.true.) |
---|
| 116 | !> |
---|
| 117 | !> will correlate A and B. |
---|
| 118 | !> |
---|
| 119 | !> |
---|
| 120 | !> Remarks: |
---|
| 121 | !> |
---|
| 122 | !> Following changes have been introduced with respect to fftn.c: |
---|
| 123 | !> - complex arguments and results are of type complex, rather than |
---|
| 124 | !> real an imaginary part separately. |
---|
| 125 | !> - increment parameter (magnitude of isign) has been dropped, |
---|
| 126 | !> inc is always one, direction of transform is given by inv. |
---|
| 127 | !> - maxf and maxp have been dropped. The amount of temporary storage |
---|
| 128 | !> needed is determined by the fftradix routine. Both fftn and fft |
---|
| 129 | !> can handle any size of array. (Maybe they take a lot of time and |
---|
| 130 | !> memory, but they will do it) |
---|
| 131 | !> |
---|
| 132 | !> Redesigning fftradix in a way, that it handles assumed shape arrays |
---|
| 133 | !> would have been desirable. However, I found it rather hard to do this |
---|
| 134 | !> in an efficient way. Problems were: |
---|
| 135 | !> - to prevent stride multiplications when indexing arrays. At least our |
---|
| 136 | !> compiler was not clever enough to discover that in fact additions |
---|
| 137 | !> would do the job as well. On the other hand, I haven't been clever |
---|
| 138 | !> enough to find an implementation using array operations. |
---|
| 139 | !> - fftradix is rather large and different versions would be necessaray |
---|
| 140 | !> for each possible rank of array. |
---|
| 141 | !> Consequently, in place transformation still needs the argument stored |
---|
| 142 | !> in a consecutive bunch of memory and can't be performed on array |
---|
| 143 | !> sections like A(100:199:-3, 50:1020). Calling fftn with such sections |
---|
| 144 | !> will most probably imply copy-in-copy-out. However, the function fft |
---|
| 145 | !> works with everything it gets and should be convenient to use. |
---|
| 146 | !> |
---|
| 147 | !> Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de> |
---|
| 148 | !> Restructured fftradix for better optimization. M. Steffens, 4 June 1997 |
---|
[2000] | 149 | !------------------------------------------------------------------------------! |
---|
[1682] | 150 | MODULE singleton |
---|
| 151 | |
---|
[1] | 152 | |
---|
[1320] | 153 | USE kinds |
---|
| 154 | |
---|
[1] | 155 | IMPLICIT NONE |
---|
| 156 | |
---|
| 157 | PRIVATE |
---|
[1320] | 158 | PUBLIC:: fft, fftn |
---|
[1] | 159 | |
---|
[1320] | 160 | REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp |
---|
| 161 | REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp |
---|
| 162 | REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp |
---|
| 163 | REAL(wp), PARAMETER:: pi = 3.14159265358979323_wp |
---|
[1] | 164 | |
---|
| 165 | INTERFACE fft |
---|
| 166 | MODULE PROCEDURE fft1d |
---|
| 167 | MODULE PROCEDURE fft2d |
---|
| 168 | MODULE PROCEDURE fft3d |
---|
| 169 | MODULE PROCEDURE fft4d |
---|
| 170 | MODULE PROCEDURE fft5d |
---|
| 171 | MODULE PROCEDURE fft6d |
---|
| 172 | MODULE PROCEDURE fft7d |
---|
| 173 | END INTERFACE |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | CONTAINS |
---|
| 177 | |
---|
| 178 | |
---|
[1682] | 179 | !------------------------------------------------------------------------------! |
---|
| 180 | ! Description: |
---|
| 181 | ! ------------ |
---|
| 182 | !> @todo Missing function description. |
---|
| 183 | !------------------------------------------------------------------------------! |
---|
[1] | 184 | FUNCTION fft1d(array, dim, inv, stat) RESULT(ft) |
---|
| 185 | ! |
---|
| 186 | !-- Formal parameters |
---|
[1320] | 187 | COMPLEX(wp), DIMENSION(:), INTENT(IN) :: array |
---|
| 188 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 189 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 190 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 191 | ! |
---|
| 192 | !-- Function result |
---|
[1320] | 193 | COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft |
---|
[1] | 194 | |
---|
[3761] | 195 | INTEGER(iwp):: idum |
---|
[1320] | 196 | INTEGER(iwp):: ishape(1) |
---|
[1] | 197 | |
---|
| 198 | ! |
---|
| 199 | !-- Intrinsics used |
---|
| 200 | INTRINSIC SIZE, SHAPE |
---|
| 201 | |
---|
| 202 | ft = array |
---|
| 203 | ishape = SHAPE( array ) |
---|
| 204 | CALL fftn(ft, ishape, inv = inv, stat = stat) |
---|
[3761] | 205 | ! |
---|
| 206 | !-- Next statement to prevent compiler warning about unused variable |
---|
| 207 | IF ( PRESENT( dim ) ) idum = 1 |
---|
[1] | 208 | |
---|
| 209 | END FUNCTION fft1d |
---|
| 210 | |
---|
| 211 | |
---|
[1682] | 212 | !------------------------------------------------------------------------------! |
---|
| 213 | ! Description: |
---|
| 214 | ! ------------ |
---|
| 215 | !> @todo Missing function description. |
---|
| 216 | !------------------------------------------------------------------------------! |
---|
[1] | 217 | FUNCTION fft2d(array, dim, inv, stat) RESULT(ft) |
---|
| 218 | ! |
---|
| 219 | !-- Formal parameters |
---|
[1320] | 220 | COMPLEX(wp), DIMENSION(:,:), INTENT(IN) :: array |
---|
| 221 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 222 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 223 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 224 | ! |
---|
| 225 | !-- Function result |
---|
[1320] | 226 | COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft |
---|
[1] | 227 | |
---|
[1320] | 228 | INTEGER(iwp) :: ishape(2) |
---|
[1] | 229 | ! |
---|
| 230 | !-- Intrinsics used |
---|
| 231 | INTRINSIC SIZE, SHAPE |
---|
| 232 | |
---|
| 233 | ft = array |
---|
| 234 | ishape = SHAPE( array ) |
---|
| 235 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 236 | |
---|
| 237 | END FUNCTION fft2d |
---|
| 238 | |
---|
| 239 | |
---|
[1682] | 240 | !------------------------------------------------------------------------------! |
---|
| 241 | ! Description: |
---|
| 242 | ! ------------ |
---|
| 243 | !> @todo Missing function description. |
---|
| 244 | !------------------------------------------------------------------------------! |
---|
[1] | 245 | FUNCTION fft3d(array, dim, inv, stat) RESULT(ft) |
---|
| 246 | ! |
---|
| 247 | !-- Formal parameters |
---|
[1320] | 248 | COMPLEX(wp), DIMENSION(:,:,:), INTENT(IN) :: array |
---|
| 249 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 250 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 251 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 252 | ! |
---|
| 253 | !-- Function result |
---|
[1320] | 254 | COMPLEX(wp), & |
---|
[1] | 255 | DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft |
---|
| 256 | |
---|
[1320] | 257 | INTEGER(iwp) :: ishape(3) |
---|
[1] | 258 | |
---|
| 259 | ! |
---|
| 260 | !-- Intrinsics used |
---|
| 261 | INTRINSIC SIZE, SHAPE |
---|
| 262 | |
---|
| 263 | ft = array |
---|
| 264 | ishape = SHAPE( array) |
---|
| 265 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 266 | |
---|
| 267 | END FUNCTION fft3d |
---|
| 268 | |
---|
| 269 | |
---|
[1682] | 270 | !------------------------------------------------------------------------------! |
---|
| 271 | ! Description: |
---|
| 272 | ! ------------ |
---|
| 273 | !> @todo Missing function description. |
---|
| 274 | !------------------------------------------------------------------------------! |
---|
[1] | 275 | FUNCTION fft4d(array, dim, inv, stat) RESULT(ft) |
---|
| 276 | ! |
---|
| 277 | !-- Formal parameters |
---|
[1320] | 278 | COMPLEX(wp), DIMENSION(:,:,:,:), INTENT(IN) :: array |
---|
| 279 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 280 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 281 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 282 | ! |
---|
| 283 | !-- Function result |
---|
[1320] | 284 | COMPLEX(wp), DIMENSION( & |
---|
[1] | 285 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft |
---|
| 286 | |
---|
[1320] | 287 | INTEGER(iwp) :: ishape(4) |
---|
[1] | 288 | ! |
---|
| 289 | !-- Intrinsics used |
---|
| 290 | INTRINSIC SIZE, SHAPE |
---|
| 291 | |
---|
| 292 | ft = array |
---|
| 293 | ishape = SHAPE( array ) |
---|
| 294 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 295 | |
---|
| 296 | END FUNCTION fft4d |
---|
| 297 | |
---|
| 298 | |
---|
[1682] | 299 | !------------------------------------------------------------------------------! |
---|
| 300 | ! Description: |
---|
| 301 | ! ------------ |
---|
| 302 | !> @todo Missing function description. |
---|
| 303 | !------------------------------------------------------------------------------! |
---|
[1] | 304 | FUNCTION fft5d(array, dim, inv, stat) RESULT(ft) |
---|
| 305 | ! |
---|
| 306 | !-- Formal parameters |
---|
[1320] | 307 | COMPLEX(wp), DIMENSION(:,:,:,:,:), INTENT(IN) :: array |
---|
| 308 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 309 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 310 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 311 | ! |
---|
| 312 | !-- Function result |
---|
[1320] | 313 | COMPLEX(wp), DIMENSION( & |
---|
[1] | 314 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
| 315 | SIZE(array, 5)):: ft |
---|
| 316 | |
---|
[1320] | 317 | INTEGER(iwp) :: ishape(5) |
---|
[1] | 318 | |
---|
| 319 | ! |
---|
| 320 | !-- Intrinsics used |
---|
| 321 | INTRINSIC SIZE, SHAPE |
---|
| 322 | |
---|
| 323 | ft = array |
---|
| 324 | ishape = SHAPE( array ) |
---|
| 325 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 326 | |
---|
| 327 | END FUNCTION fft5d |
---|
| 328 | |
---|
| 329 | |
---|
[1682] | 330 | !------------------------------------------------------------------------------! |
---|
| 331 | ! Description: |
---|
| 332 | ! ------------ |
---|
| 333 | !> @todo Missing function description. |
---|
| 334 | !------------------------------------------------------------------------------! |
---|
[1] | 335 | FUNCTION fft6d(array, dim, inv, stat) RESULT(ft) |
---|
| 336 | ! |
---|
| 337 | !-- Formal parameters |
---|
[1320] | 338 | COMPLEX(wp), DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array |
---|
| 339 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 340 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 341 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 342 | ! |
---|
| 343 | !-- Function result |
---|
[1320] | 344 | COMPLEX(wp), DIMENSION( & |
---|
[1] | 345 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
| 346 | SIZE(array, 5), SIZE(array, 6)):: ft |
---|
| 347 | |
---|
[1320] | 348 | INTEGER(iwp) :: ishape(6) |
---|
[1] | 349 | |
---|
| 350 | ! |
---|
| 351 | !-- Intrinsics used |
---|
| 352 | INTRINSIC SIZE, SHAPE |
---|
| 353 | |
---|
| 354 | ft = array |
---|
| 355 | ishape = SHAPE( array ) |
---|
| 356 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 357 | |
---|
| 358 | END FUNCTION fft6d |
---|
| 359 | |
---|
| 360 | |
---|
[1682] | 361 | !------------------------------------------------------------------------------! |
---|
| 362 | ! Description: |
---|
| 363 | ! ------------ |
---|
| 364 | !> @todo Missing function description. |
---|
| 365 | !------------------------------------------------------------------------------! |
---|
[1] | 366 | FUNCTION fft7d(array, dim, inv, stat) RESULT(ft) |
---|
| 367 | ! |
---|
| 368 | !-- Formal parameters |
---|
[1320] | 369 | COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array |
---|
| 370 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 371 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 372 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 373 | ! |
---|
| 374 | !-- Function result |
---|
[1320] | 375 | COMPLEX(wp), DIMENSION( & |
---|
[1] | 376 | SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & |
---|
| 377 | SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft |
---|
| 378 | |
---|
[1320] | 379 | INTEGER(iwp) :: ishape(7) |
---|
[1] | 380 | |
---|
| 381 | ! |
---|
| 382 | !-- Intrinsics used |
---|
| 383 | INTRINSIC SIZE, SHAPE |
---|
| 384 | |
---|
| 385 | ft = array |
---|
| 386 | ishape = SHAPE( array ) |
---|
| 387 | CALL fftn(ft, ishape, dim, inv, stat) |
---|
| 388 | |
---|
| 389 | END FUNCTION fft7d |
---|
| 390 | |
---|
| 391 | |
---|
[1682] | 392 | !------------------------------------------------------------------------------! |
---|
| 393 | ! Description: |
---|
| 394 | ! ------------ |
---|
| 395 | !> @todo Missing subroutine description. |
---|
| 396 | !------------------------------------------------------------------------------! |
---|
[1] | 397 | SUBROUTINE fftn(array, shape, dim, inv, stat) |
---|
| 398 | ! |
---|
| 399 | !-- Formal parameters |
---|
[1320] | 400 | COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array |
---|
| 401 | INTEGER(iwp), DIMENSION(:), INTENT(IN) :: shape |
---|
| 402 | INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim |
---|
| 403 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 404 | LOGICAL, INTENT(IN), OPTIONAL:: inv |
---|
[1] | 405 | ! |
---|
| 406 | !-- Local arrays |
---|
[1320] | 407 | INTEGER(iwp), DIMENSION(SIZE(shape)):: d |
---|
[1] | 408 | ! |
---|
| 409 | !-- Local scalars |
---|
| 410 | LOGICAL :: inverse |
---|
[1320] | 411 | INTEGER(iwp) :: i, ndim, ntotal |
---|
| 412 | REAL(wp):: scale |
---|
[1] | 413 | ! |
---|
| 414 | !-- Intrinsics used |
---|
| 415 | INTRINSIC PRESENT, MIN, PRODUCT, SIZE, SQRT |
---|
| 416 | |
---|
| 417 | ! |
---|
| 418 | !-- Optional parameter settings |
---|
| 419 | IF (PRESENT(inv)) THEN |
---|
| 420 | inverse = inv |
---|
| 421 | ELSE |
---|
| 422 | inverse = .FALSE. |
---|
| 423 | END IF |
---|
| 424 | IF (PRESENT(dim)) THEN |
---|
| 425 | ndim = MIN(SIZE(dim), SIZE(d)) |
---|
| 426 | d(1:ndim) = DIM(1:ndim) |
---|
| 427 | ELSE |
---|
| 428 | ndim = SIZE(d) |
---|
| 429 | d = (/(i, i = 1, SIZE(d))/) |
---|
| 430 | END IF |
---|
| 431 | |
---|
| 432 | ntotal = PRODUCT(shape) |
---|
[1320] | 433 | scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim)))) |
---|
[1] | 434 | DO i = 1, ntotal |
---|
| 435 | array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, & |
---|
[1320] | 436 | KIND=wp) |
---|
[1] | 437 | END DO |
---|
| 438 | |
---|
| 439 | DO i = 1, ndim |
---|
| 440 | CALL fftradix(array, ntotal, shape(d(i)), PRODUCT(shape(1:d(i))), & |
---|
| 441 | inverse, stat) |
---|
| 442 | IF (PRESENT(stat)) THEN |
---|
| 443 | IF (stat /=0) RETURN |
---|
| 444 | END IF |
---|
| 445 | END DO |
---|
| 446 | |
---|
| 447 | END SUBROUTINE fftn |
---|
| 448 | |
---|
| 449 | |
---|
[1682] | 450 | !------------------------------------------------------------------------------! |
---|
| 451 | ! Description: |
---|
| 452 | ! ------------ |
---|
| 453 | !> @todo Missing subroutine description. |
---|
| 454 | !------------------------------------------------------------------------------! |
---|
[1] | 455 | SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat) |
---|
| 456 | ! |
---|
| 457 | !-- Formal parameters |
---|
[1320] | 458 | COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array |
---|
| 459 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
| 460 | INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat |
---|
| 461 | LOGICAL, INTENT(IN) :: inv |
---|
[1] | 462 | ! |
---|
| 463 | !-- Local arrays |
---|
[1320] | 464 | COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: ctmp |
---|
| 465 | INTEGER(iwp), DIMENSION(BIT_SIZE(0)) :: factor |
---|
| 466 | INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: perm |
---|
| 467 | REAL(wp), DIMENSION(:), ALLOCATABLE :: sine, cosine |
---|
[1] | 468 | ! |
---|
| 469 | !-- Local scalars |
---|
[1320] | 470 | INTEGER(iwp) :: maxfactor, nfactor, nsquare, nperm |
---|
[1] | 471 | ! |
---|
| 472 | !-- Intrinsics used |
---|
| 473 | INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, & |
---|
| 474 | CMPLX, REAL, AIMAG |
---|
| 475 | |
---|
| 476 | IF (npass <= 1) RETURN |
---|
| 477 | |
---|
| 478 | CALL factorize(npass, factor, nfactor, nsquare) |
---|
| 479 | |
---|
| 480 | maxfactor = MAXVAL(factor(:nfactor)) |
---|
| 481 | IF (nfactor - ISHFT(nsquare, 1) > 0) THEN |
---|
| 482 | nperm = MAX(nfactor + 1, PRODUCT(factor(nsquare+1: nfactor-nsquare)) - 1) |
---|
| 483 | ELSE |
---|
| 484 | nperm = nfactor + 1 |
---|
| 485 | END IF |
---|
| 486 | |
---|
| 487 | IF (PRESENT(stat)) THEN |
---|
| 488 | ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT=stat) |
---|
| 489 | IF (stat /= 0) RETURN |
---|
| 490 | CALL transform(array, ntotal, npass, nspan, & |
---|
| 491 | factor, nfactor, ctmp, sine, cosine, inv) |
---|
| 492 | DEALLOCATE(sine, cosine, STAT=stat) |
---|
| 493 | IF (stat /= 0) RETURN |
---|
| 494 | ALLOCATE(perm(nperm), STAT=stat) |
---|
| 495 | IF (stat /= 0) RETURN |
---|
| 496 | CALL permute(array, ntotal, npass, nspan, & |
---|
| 497 | factor, nfactor, nsquare, maxfactor, & |
---|
| 498 | ctmp, perm) |
---|
| 499 | DEALLOCATE(perm, ctmp, STAT=stat) |
---|
| 500 | IF (stat /= 0) RETURN |
---|
| 501 | ELSE |
---|
| 502 | ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor)) |
---|
| 503 | CALL transform(array, ntotal, npass, nspan, & |
---|
| 504 | factor, nfactor, ctmp, sine, cosine, inv) |
---|
| 505 | DEALLOCATE(sine, cosine) |
---|
| 506 | ALLOCATE(perm(nperm)) |
---|
| 507 | CALL permute(array, ntotal, npass, nspan, & |
---|
| 508 | factor, nfactor, nsquare, maxfactor, & |
---|
| 509 | ctmp, perm) |
---|
| 510 | DEALLOCATE(perm, ctmp) |
---|
| 511 | END IF |
---|
| 512 | |
---|
| 513 | |
---|
| 514 | CONTAINS |
---|
| 515 | |
---|
| 516 | |
---|
[1682] | 517 | !------------------------------------------------------------------------------! |
---|
| 518 | ! Description: |
---|
| 519 | ! ------------ |
---|
| 520 | !> @todo Missing subroutine description. |
---|
| 521 | !------------------------------------------------------------------------------! |
---|
[1] | 522 | SUBROUTINE factorize(npass, factor, nfactor, nsquare) |
---|
| 523 | ! |
---|
| 524 | !-- Formal parameters |
---|
[1320] | 525 | INTEGER(iwp), INTENT(IN) :: npass |
---|
| 526 | INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor |
---|
| 527 | INTEGER(iwp), INTENT(OUT):: nfactor, nsquare |
---|
[1] | 528 | ! |
---|
| 529 | !-- Local scalars |
---|
[1320] | 530 | INTEGER(iwp):: j, jj, k |
---|
[1] | 531 | |
---|
| 532 | nfactor = 0 |
---|
| 533 | k = npass |
---|
| 534 | DO WHILE (MOD(k, 16) == 0) |
---|
| 535 | nfactor = nfactor + 1 |
---|
| 536 | factor(nfactor) = 4 |
---|
| 537 | k = k / 16 |
---|
| 538 | END DO |
---|
| 539 | j = 3 |
---|
| 540 | jj = 9 |
---|
| 541 | DO |
---|
| 542 | DO WHILE (MOD(k, jj) == 0) |
---|
| 543 | nfactor = nfactor + 1 |
---|
| 544 | factor(nfactor) = j |
---|
| 545 | k = k / jj |
---|
| 546 | END DO |
---|
| 547 | j = j + 2 |
---|
| 548 | jj = j * j |
---|
| 549 | IF (jj > k) EXIT |
---|
| 550 | END DO |
---|
| 551 | IF (k <= 4) THEN |
---|
| 552 | nsquare = nfactor |
---|
| 553 | factor(nfactor + 1) = k |
---|
| 554 | IF (k /= 1) nfactor = nfactor + 1 |
---|
| 555 | ELSE |
---|
| 556 | IF (k - ISHFT(k / 4, 2) == 0) THEN |
---|
| 557 | nfactor = nfactor + 1 |
---|
| 558 | factor(nfactor) = 2 |
---|
| 559 | k = k / 4 |
---|
| 560 | END IF |
---|
| 561 | nsquare = nfactor |
---|
| 562 | j = 2 |
---|
| 563 | DO |
---|
| 564 | IF (MOD(k, j) == 0) THEN |
---|
| 565 | nfactor = nfactor + 1 |
---|
| 566 | factor(nfactor) = j |
---|
| 567 | k = k / j |
---|
| 568 | END IF |
---|
| 569 | j = ISHFT((j + 1) / 2, 1) + 1 |
---|
| 570 | IF (j > k) EXIT |
---|
| 571 | END DO |
---|
| 572 | END IF |
---|
| 573 | IF (nsquare > 0) THEN |
---|
| 574 | j = nsquare |
---|
| 575 | DO |
---|
| 576 | nfactor = nfactor + 1 |
---|
| 577 | factor(nfactor) = factor(j) |
---|
| 578 | j = j - 1 |
---|
| 579 | IF (j==0) EXIT |
---|
| 580 | END DO |
---|
| 581 | END IF |
---|
| 582 | |
---|
| 583 | END SUBROUTINE factorize |
---|
| 584 | |
---|
| 585 | |
---|
[1682] | 586 | !------------------------------------------------------------------------------! |
---|
| 587 | ! Description: |
---|
| 588 | ! ------------ |
---|
| 589 | !> @todo Missing subroutine description. |
---|
| 590 | !------------------------------------------------------------------------------! |
---|
[1] | 591 | SUBROUTINE transform(array, ntotal, npass, nspan, & |
---|
| 592 | factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform |
---|
| 593 | ! |
---|
| 594 | !-- Formal parameters |
---|
[1320] | 595 | COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array |
---|
| 596 | COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp |
---|
| 597 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
| 598 | INTEGER(iwp), DIMENSION(*), INTENT(IN) :: factor |
---|
| 599 | INTEGER(iwp), INTENT(IN) :: nfactor |
---|
| 600 | LOGICAL, INTENT(IN) :: inv |
---|
| 601 | REAL(wp), DIMENSION(*), INTENT(OUT) :: sine, cosine |
---|
[1] | 602 | ! |
---|
| 603 | !-- Local scalars |
---|
[1320] | 604 | INTEGER(iwp):: ii, ispan |
---|
| 605 | INTEGER(iwp):: j, jc, jf, jj |
---|
| 606 | INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4 |
---|
| 607 | INTEGER(iwp):: nn, nt |
---|
| 608 | REAL(wp) :: s60, c72, s72, pi2, radf |
---|
| 609 | REAL(wp) :: c1, s1, c2, s2, c3, s3, cd, sd, ak |
---|
| 610 | COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm |
---|
[1] | 611 | |
---|
| 612 | c72 = cos72 |
---|
| 613 | IF (inv) THEN |
---|
| 614 | s72 = sin72 |
---|
| 615 | s60 = sin60 |
---|
| 616 | pi2 = pi |
---|
| 617 | ELSE |
---|
| 618 | s72 = -sin72 |
---|
| 619 | s60 = -sin60 |
---|
| 620 | pi2 = -pi |
---|
| 621 | END IF |
---|
| 622 | |
---|
| 623 | nt = ntotal |
---|
| 624 | nn = nt - 1 |
---|
| 625 | kspan = nspan |
---|
| 626 | jc = nspan / npass |
---|
| 627 | radf = pi2 * jc |
---|
[1320] | 628 | pi2 = pi2 * 2.0_wp !-- use 2 PI from here on |
---|
[1] | 629 | |
---|
| 630 | ii = 0 |
---|
| 631 | jf = 0 |
---|
| 632 | DO |
---|
| 633 | sd = radf / kspan |
---|
| 634 | cd = SIN(sd) |
---|
[1320] | 635 | cd = 2.0_wp * cd * cd |
---|
[1] | 636 | sd = SIN(sd + sd) |
---|
| 637 | kk = 1 |
---|
| 638 | ii = ii + 1 |
---|
| 639 | |
---|
| 640 | SELECT CASE (factor(ii)) |
---|
| 641 | CASE (2) |
---|
| 642 | ! |
---|
| 643 | !-- Transform for factor of 2 (including rotation factor) |
---|
| 644 | kspan = kspan / 2 |
---|
| 645 | k1 = kspan + 2 |
---|
| 646 | DO |
---|
| 647 | DO |
---|
| 648 | k2 = kk + kspan |
---|
| 649 | ck = array(k2) |
---|
| 650 | array(k2) = array(kk)-ck |
---|
| 651 | array(kk) = array(kk) + ck |
---|
| 652 | kk = k2 + kspan |
---|
| 653 | IF (kk > nn) EXIT |
---|
| 654 | END DO |
---|
| 655 | kk = kk - nn |
---|
| 656 | IF (kk > jc) EXIT |
---|
| 657 | END DO |
---|
| 658 | IF (kk > kspan) RETURN |
---|
| 659 | DO |
---|
[1320] | 660 | c1 = 1.0_wp - cd |
---|
[1] | 661 | s1 = sd |
---|
| 662 | DO |
---|
| 663 | DO |
---|
| 664 | DO |
---|
| 665 | k2 = kk + kspan |
---|
| 666 | ck = array(kk) - array(k2) |
---|
| 667 | array(kk) = array(kk) + array(k2) |
---|
[1320] | 668 | array(k2) = ck * CMPLX(c1, s1, KIND=wp) |
---|
[1] | 669 | kk = k2 + kspan |
---|
| 670 | IF (kk >= nt) EXIT |
---|
| 671 | END DO |
---|
| 672 | k2 = kk - nt |
---|
| 673 | c1 = -c1 |
---|
| 674 | kk = k1 - k2 |
---|
| 675 | IF (kk <= k2) EXIT |
---|
| 676 | END DO |
---|
| 677 | ak = c1 - (cd * c1 + sd * s1) |
---|
| 678 | s1 = sd * c1 - cd * s1 + s1 |
---|
[1320] | 679 | c1 = 2.0_wp - (ak * ak + s1 * s1) |
---|
[1] | 680 | s1 = s1 * c1 |
---|
| 681 | c1 = c1 * ak |
---|
| 682 | kk = kk + jc |
---|
| 683 | IF (kk >= k2) EXIT |
---|
| 684 | END DO |
---|
| 685 | k1 = k1 + 1 + 1 |
---|
| 686 | kk = (k1 - kspan) / 2 + jc |
---|
| 687 | IF (kk > jc + jc) EXIT |
---|
| 688 | END DO |
---|
| 689 | |
---|
| 690 | CASE (4) !-- transform for factor of 4 |
---|
| 691 | ispan = kspan |
---|
| 692 | kspan = kspan / 4 |
---|
| 693 | |
---|
| 694 | DO |
---|
[1320] | 695 | c1 = 1.0_wp |
---|
| 696 | s1 = 0.0_wp |
---|
[1] | 697 | DO |
---|
| 698 | DO |
---|
| 699 | k1 = kk + kspan |
---|
| 700 | k2 = k1 + kspan |
---|
| 701 | k3 = k2 + kspan |
---|
| 702 | ckp = array(kk) + array(k2) |
---|
| 703 | ckm = array(kk) - array(k2) |
---|
| 704 | cjp = array(k1) + array(k3) |
---|
| 705 | cjm = array(k1) - array(k3) |
---|
| 706 | array(kk) = ckp + cjp |
---|
| 707 | cjp = ckp - cjp |
---|
| 708 | IF (inv) THEN |
---|
[1320] | 709 | ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) |
---|
| 710 | ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) |
---|
[1] | 711 | ELSE |
---|
[1320] | 712 | ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) |
---|
| 713 | ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) |
---|
[1] | 714 | END IF |
---|
| 715 | ! |
---|
| 716 | !-- Avoid useless multiplies |
---|
[1320] | 717 | IF (s1 == 0.0_wp) THEN |
---|
[1] | 718 | array(k1) = ckp |
---|
| 719 | array(k2) = cjp |
---|
| 720 | array(k3) = ckm |
---|
| 721 | ELSE |
---|
[1320] | 722 | array(k1) = ckp * CMPLX(c1, s1, KIND=wp) |
---|
| 723 | array(k2) = cjp * CMPLX(c2, s2, KIND=wp) |
---|
| 724 | array(k3) = ckm * CMPLX(c3, s3, KIND=wp) |
---|
[1] | 725 | END IF |
---|
| 726 | kk = k3 + kspan |
---|
| 727 | IF (kk > nt) EXIT |
---|
| 728 | END DO |
---|
| 729 | |
---|
| 730 | c2 = c1 - (cd * c1 + sd * s1) |
---|
| 731 | s1 = sd * c1 - cd * s1 + s1 |
---|
[1320] | 732 | c1 = 2.0_wp - (c2 * c2 + s1 * s1) |
---|
[1] | 733 | s1 = s1 * c1 |
---|
| 734 | c1 = c1 * c2 |
---|
| 735 | ! |
---|
| 736 | !-- Values of c2, c3, s2, s3 that will get used next time |
---|
| 737 | c2 = c1 * c1 - s1 * s1 |
---|
[1320] | 738 | s2 = 2.0_wp * c1 * s1 |
---|
[1] | 739 | c3 = c2 * c1 - s2 * s1 |
---|
| 740 | s3 = c2 * s1 + s2 * c1 |
---|
| 741 | kk = kk - nt + jc |
---|
| 742 | IF (kk > kspan) EXIT |
---|
| 743 | END DO |
---|
| 744 | kk = kk - kspan + 1 |
---|
| 745 | IF (kk > jc) EXIT |
---|
| 746 | END DO |
---|
| 747 | IF (kspan == jc) RETURN |
---|
| 748 | |
---|
| 749 | CASE default |
---|
| 750 | ! |
---|
| 751 | !-- Transform for odd factors |
---|
| 752 | k = factor(ii) |
---|
| 753 | ispan = kspan |
---|
| 754 | kspan = kspan / k |
---|
| 755 | |
---|
| 756 | SELECT CASE (k) |
---|
| 757 | CASE (3) !-- transform for factor of 3 (optional code) |
---|
| 758 | DO |
---|
| 759 | DO |
---|
| 760 | k1 = kk + kspan |
---|
| 761 | k2 = k1 + kspan |
---|
| 762 | ck = array(kk) |
---|
| 763 | cj = array(k1) + array(k2) |
---|
| 764 | array(kk) = ck + cj |
---|
| 765 | ck = ck - CMPLX( & |
---|
[1320] | 766 | 0.5_wp * REAL (cj), & |
---|
| 767 | 0.5_wp * AIMAG(cj), & |
---|
| 768 | KIND=wp) |
---|
[1] | 769 | cj = CMPLX( & |
---|
| 770 | (REAL (array(k1)) - REAL (array(k2))) * s60, & |
---|
| 771 | (AIMAG(array(k1)) - AIMAG(array(k2))) * s60, & |
---|
[1320] | 772 | KIND=wp) |
---|
| 773 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
| 774 | array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
[1] | 775 | kk = k2 + kspan |
---|
| 776 | IF (kk >= nn) EXIT |
---|
| 777 | END DO |
---|
| 778 | kk = kk - nn |
---|
| 779 | IF (kk > kspan) EXIT |
---|
| 780 | END DO |
---|
| 781 | |
---|
| 782 | CASE (5) !-- transform for factor of 5 (optional code) |
---|
| 783 | c2 = c72 * c72 - s72 * s72 |
---|
[1320] | 784 | s2 = 2.0_wp * c72 * s72 |
---|
[1] | 785 | DO |
---|
| 786 | DO |
---|
| 787 | k1 = kk + kspan |
---|
| 788 | k2 = k1 + kspan |
---|
| 789 | k3 = k2 + kspan |
---|
| 790 | k4 = k3 + kspan |
---|
| 791 | ckp = array(k1) + array(k4) |
---|
| 792 | ckm = array(k1) - array(k4) |
---|
| 793 | cjp = array(k2) + array(k3) |
---|
| 794 | cjm = array(k2) - array(k3) |
---|
| 795 | cc = array(kk) |
---|
| 796 | array(kk) = cc + ckp + cjp |
---|
| 797 | ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, & |
---|
[1320] | 798 | KIND=wp) + & |
---|
[1] | 799 | CMPLX(REAL(cjp) * c2, AIMAG(cjp) * c2, & |
---|
[1320] | 800 | KIND=wp) + cc |
---|
[1] | 801 | cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, & |
---|
[1320] | 802 | KIND=wp) + & |
---|
[1] | 803 | CMPLX(REAL(cjm) * s2, AIMAG(cjm) * s2, & |
---|
[1320] | 804 | KIND=wp) |
---|
| 805 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
| 806 | array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
[1] | 807 | ck = CMPLX(REAL(ckp) * c2, AIMAG(ckp) * c2, & |
---|
[1320] | 808 | KIND=wp) + & |
---|
[1] | 809 | CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, & |
---|
[1320] | 810 | KIND=wp) + cc |
---|
[1] | 811 | cj = CMPLX(REAL(ckm) * s2, AIMAG(ckm) * s2, & |
---|
[1320] | 812 | KIND=wp) - & |
---|
[1] | 813 | CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, & |
---|
[1320] | 814 | KIND=wp) |
---|
| 815 | array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) |
---|
| 816 | array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) |
---|
[1] | 817 | kk = k4 + kspan |
---|
| 818 | IF (kk >= nn) EXIT |
---|
| 819 | END DO |
---|
| 820 | kk = kk - nn |
---|
| 821 | IF (kk > kspan) EXIT |
---|
| 822 | END DO |
---|
| 823 | |
---|
| 824 | CASE default |
---|
| 825 | IF (k /= jf) THEN |
---|
| 826 | jf = k |
---|
| 827 | s1 = pi2 / k |
---|
| 828 | c1 = COS(s1) |
---|
| 829 | s1 = SIN(s1) |
---|
[1320] | 830 | cosine (jf) = 1.0_wp |
---|
| 831 | sine (jf) = 0.0_wp |
---|
[1] | 832 | j = 1 |
---|
| 833 | DO |
---|
| 834 | cosine (j) = cosine (k) * c1 + sine (k) * s1 |
---|
| 835 | sine (j) = cosine (k) * s1 - sine (k) * c1 |
---|
| 836 | k = k-1 |
---|
| 837 | cosine (k) = cosine (j) |
---|
| 838 | sine (k) = -sine (j) |
---|
| 839 | j = j + 1 |
---|
| 840 | IF (j >= k) EXIT |
---|
| 841 | END DO |
---|
| 842 | END IF |
---|
| 843 | DO |
---|
| 844 | DO |
---|
| 845 | k1 = kk |
---|
| 846 | k2 = kk + ispan |
---|
| 847 | cc = array(kk) |
---|
| 848 | ck = cc |
---|
| 849 | j = 1 |
---|
| 850 | k1 = k1 + kspan |
---|
| 851 | DO |
---|
| 852 | k2 = k2 - kspan |
---|
| 853 | j = j + 1 |
---|
| 854 | ctmp(j) = array(k1) + array(k2) |
---|
| 855 | ck = ck + ctmp(j) |
---|
| 856 | j = j + 1 |
---|
| 857 | ctmp(j) = array(k1) - array(k2) |
---|
| 858 | k1 = k1 + kspan |
---|
| 859 | IF (k1 >= k2) EXIT |
---|
| 860 | END DO |
---|
| 861 | array(kk) = ck |
---|
| 862 | k1 = kk |
---|
| 863 | k2 = kk + ispan |
---|
| 864 | j = 1 |
---|
| 865 | DO |
---|
| 866 | k1 = k1 + kspan |
---|
| 867 | k2 = k2 - kspan |
---|
| 868 | jj = j |
---|
| 869 | ck = cc |
---|
[1320] | 870 | cj = (0.0_wp, 0.0_wp) |
---|
[1] | 871 | k = 1 |
---|
| 872 | DO |
---|
| 873 | k = k + 1 |
---|
| 874 | ck = ck + CMPLX( & |
---|
| 875 | REAL (ctmp(k)) * cosine(jj), & |
---|
[1320] | 876 | AIMAG(ctmp(k)) * cosine(jj), KIND=wp) |
---|
[1] | 877 | k = k + 1 |
---|
| 878 | cj = cj + CMPLX( & |
---|
| 879 | REAL (ctmp(k)) * sine(jj), & |
---|
[1320] | 880 | AIMAG(ctmp(k)) * sine(jj), KIND=wp) |
---|
[1] | 881 | jj = jj + j |
---|
| 882 | IF (jj > jf) jj = jj - jf |
---|
| 883 | IF (k >= jf) EXIT |
---|
| 884 | END DO |
---|
| 885 | k = jf - j |
---|
| 886 | array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), & |
---|
[1320] | 887 | KIND=wp) |
---|
[1] | 888 | array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), & |
---|
[1320] | 889 | KIND=wp) |
---|
[1] | 890 | j = j + 1 |
---|
| 891 | IF (j >= k) EXIT |
---|
| 892 | END DO |
---|
| 893 | kk = kk + ispan |
---|
| 894 | IF (kk > nn) EXIT |
---|
| 895 | END DO |
---|
| 896 | kk = kk - nn |
---|
| 897 | IF (kk > kspan) EXIT |
---|
| 898 | END DO |
---|
| 899 | |
---|
| 900 | END SELECT |
---|
| 901 | ! |
---|
| 902 | !-- Multiply by rotation factor (except for factors of 2 and 4) |
---|
| 903 | IF (ii == nfactor) RETURN |
---|
| 904 | kk = jc + 1 |
---|
| 905 | DO |
---|
[1320] | 906 | c2 = 1.0_wp - cd |
---|
[1] | 907 | s1 = sd |
---|
| 908 | DO |
---|
| 909 | c1 = c2 |
---|
| 910 | s2 = s1 |
---|
| 911 | kk = kk + kspan |
---|
| 912 | DO |
---|
| 913 | DO |
---|
[1320] | 914 | array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk) |
---|
[1] | 915 | kk = kk + ispan |
---|
| 916 | IF (kk > nt) EXIT |
---|
| 917 | END DO |
---|
| 918 | ak = s1 * s2 |
---|
| 919 | s2 = s1 * c2 + c1 * s2 |
---|
| 920 | c2 = c1 * c2 - ak |
---|
| 921 | kk = kk - nt + kspan |
---|
| 922 | IF (kk > ispan) EXIT |
---|
| 923 | END DO |
---|
| 924 | c2 = c1 - (cd * c1 + sd * s1) |
---|
| 925 | s1 = s1 + sd * c1 - cd * s1 |
---|
[1320] | 926 | c1 = 2.0_wp - (c2 * c2 + s1 * s1) |
---|
[1] | 927 | s1 = s1 * c1 |
---|
| 928 | c2 = c2 * c1 |
---|
| 929 | kk = kk - ispan + jc |
---|
| 930 | IF (kk > kspan) EXIT |
---|
| 931 | END DO |
---|
| 932 | kk = kk - kspan + jc + 1 |
---|
| 933 | IF (kk > jc + jc) EXIT |
---|
| 934 | END DO |
---|
| 935 | |
---|
| 936 | END SELECT |
---|
| 937 | END DO |
---|
| 938 | END SUBROUTINE transform |
---|
| 939 | |
---|
| 940 | |
---|
[1682] | 941 | !------------------------------------------------------------------------------! |
---|
| 942 | ! Description: |
---|
| 943 | ! ------------ |
---|
| 944 | !> @todo Missing subroutine description. |
---|
| 945 | !------------------------------------------------------------------------------! |
---|
[1] | 946 | SUBROUTINE permute(array, ntotal, npass, nspan, & |
---|
| 947 | factor, nfactor, nsquare, maxfactor, & |
---|
| 948 | ctmp, perm) |
---|
| 949 | ! |
---|
| 950 | !-- Formal parameters |
---|
[1320] | 951 | COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array |
---|
| 952 | COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp |
---|
| 953 | INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan |
---|
| 954 | INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor |
---|
| 955 | INTEGER(iwp), INTENT(IN) :: nfactor, nsquare |
---|
| 956 | INTEGER(iwp), INTENT(IN) :: maxfactor |
---|
| 957 | INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: perm |
---|
[1] | 958 | ! |
---|
| 959 | !-- Local scalars |
---|
[1320] | 960 | COMPLEX(wp) :: ck |
---|
| 961 | INTEGER(iwp):: ii, ispan |
---|
| 962 | INTEGER(iwp):: j, jc, jj |
---|
| 963 | INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3 |
---|
| 964 | INTEGER(iwp):: nn, nt |
---|
[1] | 965 | ! |
---|
| 966 | !-- Permute the results to normal order---done in two stages |
---|
| 967 | !-- Permutation for square factors of n |
---|
| 968 | |
---|
| 969 | nt = ntotal |
---|
| 970 | nn = nt - 1 |
---|
| 971 | kt = nsquare |
---|
| 972 | kspan = nspan |
---|
| 973 | jc = nspan / npass |
---|
| 974 | |
---|
| 975 | perm (1) = nspan |
---|
| 976 | IF (kt > 0) THEN |
---|
| 977 | k = kt + kt + 1 |
---|
| 978 | IF (nfactor < k) k = k - 1 |
---|
| 979 | j = 1 |
---|
| 980 | perm (k + 1) = jc |
---|
| 981 | DO |
---|
| 982 | perm (j + 1) = perm (j) / factor(j) |
---|
| 983 | perm (k) = perm (k + 1) * factor(j) |
---|
| 984 | j = j + 1 |
---|
| 985 | k = k - 1 |
---|
| 986 | IF (j >= k) EXIT |
---|
| 987 | END DO |
---|
| 988 | k3 = perm (k + 1) |
---|
| 989 | kspan = perm (2) |
---|
| 990 | kk = jc + 1 |
---|
| 991 | k2 = kspan + 1 |
---|
| 992 | j = 1 |
---|
| 993 | |
---|
| 994 | IF (npass /= ntotal) THEN |
---|
| 995 | permute_multi: DO |
---|
| 996 | DO |
---|
| 997 | DO |
---|
| 998 | k = kk + jc |
---|
| 999 | DO |
---|
| 1000 | ! |
---|
| 1001 | !-- Swap array(kk) <> array(k2) |
---|
| 1002 | ck = array(kk) |
---|
| 1003 | array(kk) = array(k2) |
---|
| 1004 | array(k2) = ck |
---|
| 1005 | kk = kk + 1 |
---|
| 1006 | k2 = k2 + 1 |
---|
| 1007 | IF (kk >= k) EXIT |
---|
| 1008 | END DO |
---|
| 1009 | kk = kk + nspan - jc |
---|
| 1010 | k2 = k2 + nspan - jc |
---|
| 1011 | IF (kk >= nt) EXIT |
---|
| 1012 | END DO |
---|
| 1013 | kk = kk - nt + jc |
---|
| 1014 | k2 = k2 - nt + kspan |
---|
| 1015 | IF (k2 >= nspan) EXIT |
---|
| 1016 | END DO |
---|
| 1017 | DO |
---|
| 1018 | DO |
---|
| 1019 | k2 = k2 - perm (j) |
---|
| 1020 | j = j + 1 |
---|
| 1021 | k2 = perm (j + 1) + k2 |
---|
| 1022 | IF (k2 <= perm (j)) EXIT |
---|
| 1023 | END DO |
---|
| 1024 | j = 1 |
---|
| 1025 | DO |
---|
| 1026 | IF (kk < k2) CYCLE permute_multi |
---|
| 1027 | kk = kk + jc |
---|
| 1028 | k2 = k2 + kspan |
---|
| 1029 | IF (k2 >= nspan) EXIT |
---|
| 1030 | END DO |
---|
| 1031 | IF (kk >= nspan) EXIT |
---|
| 1032 | END DO |
---|
| 1033 | EXIT |
---|
| 1034 | END DO permute_multi |
---|
| 1035 | ELSE |
---|
| 1036 | permute_single: DO |
---|
| 1037 | DO |
---|
| 1038 | ! |
---|
| 1039 | !-- Swap array(kk) <> array(k2) |
---|
| 1040 | ck = array(kk) |
---|
| 1041 | array(kk) = array(k2) |
---|
| 1042 | array(k2) = ck |
---|
| 1043 | kk = kk + 1 |
---|
| 1044 | k2 = k2 + kspan |
---|
| 1045 | IF (k2 >= nspan) EXIT |
---|
| 1046 | END DO |
---|
| 1047 | DO |
---|
| 1048 | DO |
---|
| 1049 | k2 = k2 - perm (j) |
---|
| 1050 | j = j + 1 |
---|
| 1051 | k2 = perm (j + 1) + k2 |
---|
| 1052 | IF (k2 <= perm (j)) EXIT |
---|
| 1053 | END DO |
---|
| 1054 | j = 1 |
---|
| 1055 | DO |
---|
| 1056 | IF (kk < k2) CYCLE permute_single |
---|
| 1057 | kk = kk + 1 |
---|
| 1058 | k2 = k2 + kspan |
---|
| 1059 | IF (k2 >= nspan) EXIT |
---|
| 1060 | END DO |
---|
| 1061 | IF (kk >= nspan) EXIT |
---|
| 1062 | END DO |
---|
| 1063 | EXIT |
---|
| 1064 | END DO permute_single |
---|
| 1065 | END IF |
---|
| 1066 | jc = k3 |
---|
| 1067 | END IF |
---|
| 1068 | |
---|
| 1069 | IF (ISHFT(kt, 1) + 1 >= nfactor) RETURN |
---|
| 1070 | |
---|
| 1071 | ispan = perm (kt + 1) |
---|
| 1072 | ! |
---|
| 1073 | !-- Permutation for square-free factors of n |
---|
| 1074 | j = nfactor - kt |
---|
| 1075 | factor(j + 1) = 1 |
---|
| 1076 | DO |
---|
| 1077 | factor(j) = factor(j) * factor(j+1) |
---|
| 1078 | j = j - 1 |
---|
| 1079 | IF (j == kt) EXIT |
---|
| 1080 | END DO |
---|
| 1081 | kt = kt + 1 |
---|
| 1082 | nn = factor(kt) - 1 |
---|
| 1083 | j = 0 |
---|
| 1084 | jj = 0 |
---|
| 1085 | DO |
---|
| 1086 | k = kt + 1 |
---|
| 1087 | k2 = factor(kt) |
---|
| 1088 | kk = factor(k) |
---|
| 1089 | j = j + 1 |
---|
| 1090 | IF (j > nn) EXIT !-- exit infinite loop |
---|
| 1091 | jj = jj + kk |
---|
| 1092 | DO WHILE (jj >= k2) |
---|
| 1093 | jj = jj - k2 |
---|
| 1094 | k2 = kk |
---|
| 1095 | k = k + 1 |
---|
| 1096 | kk = factor(k) |
---|
| 1097 | jj = jj + kk |
---|
| 1098 | END DO |
---|
| 1099 | perm (j) = jj |
---|
| 1100 | END DO |
---|
| 1101 | ! |
---|
| 1102 | !-- Determine the permutation cycles of length greater than 1 |
---|
| 1103 | j = 0 |
---|
| 1104 | DO |
---|
| 1105 | DO |
---|
| 1106 | j = j + 1 |
---|
| 1107 | kk = perm(j) |
---|
| 1108 | IF (kk >= 0) EXIT |
---|
| 1109 | END DO |
---|
| 1110 | IF (kk /= j) THEN |
---|
| 1111 | DO |
---|
| 1112 | k = kk |
---|
| 1113 | kk = perm (k) |
---|
| 1114 | perm (k) = -kk |
---|
| 1115 | IF (kk == j) EXIT |
---|
| 1116 | END DO |
---|
| 1117 | k3 = kk |
---|
| 1118 | ELSE |
---|
| 1119 | perm (j) = -j |
---|
| 1120 | IF (j == nn) EXIT !-- exit infinite loop |
---|
| 1121 | END IF |
---|
| 1122 | END DO |
---|
| 1123 | ! |
---|
| 1124 | !-- Reorder a and b, following the permutation cycles |
---|
| 1125 | DO |
---|
| 1126 | j = k3 + 1 |
---|
| 1127 | nt = nt - ispan |
---|
| 1128 | ii = nt - 1 + 1 |
---|
| 1129 | IF (nt < 0) EXIT !-- exit infinite loop |
---|
| 1130 | DO |
---|
| 1131 | DO |
---|
| 1132 | j = j-1 |
---|
| 1133 | IF (perm(j) >= 0) EXIT |
---|
| 1134 | END DO |
---|
| 1135 | jj = jc |
---|
| 1136 | DO |
---|
| 1137 | kspan = jj |
---|
| 1138 | IF (jj > maxfactor) kspan = maxfactor |
---|
| 1139 | jj = jj - kspan |
---|
| 1140 | k = perm(j) |
---|
| 1141 | kk = jc * k + ii + jj |
---|
| 1142 | k1 = kk + kspan |
---|
| 1143 | k2 = 0 |
---|
| 1144 | DO |
---|
| 1145 | k2 = k2 + 1 |
---|
| 1146 | ctmp(k2) = array(k1) |
---|
| 1147 | k1 = k1 - 1 |
---|
| 1148 | IF (k1 == kk) EXIT |
---|
| 1149 | END DO |
---|
| 1150 | DO |
---|
| 1151 | k1 = kk + kspan |
---|
| 1152 | k2 = k1 - jc * (k + perm(k)) |
---|
| 1153 | k = -perm(k) |
---|
| 1154 | DO |
---|
| 1155 | array(k1) = array(k2) |
---|
| 1156 | k1 = k1 - 1 |
---|
| 1157 | k2 = k2 - 1 |
---|
| 1158 | IF (k1 == kk) EXIT |
---|
| 1159 | END DO |
---|
| 1160 | kk = k2 |
---|
| 1161 | IF (k == j) EXIT |
---|
| 1162 | END DO |
---|
| 1163 | k1 = kk + kspan |
---|
| 1164 | k2 = 0 |
---|
| 1165 | DO |
---|
| 1166 | k2 = k2 + 1 |
---|
| 1167 | array(k1) = ctmp(k2) |
---|
| 1168 | k1 = k1 - 1 |
---|
| 1169 | IF (k1 == kk) EXIT |
---|
| 1170 | END DO |
---|
| 1171 | IF (jj == 0) EXIT |
---|
| 1172 | END DO |
---|
| 1173 | IF (j == 1) EXIT |
---|
| 1174 | END DO |
---|
| 1175 | END DO |
---|
| 1176 | |
---|
| 1177 | END SUBROUTINE permute |
---|
| 1178 | |
---|
| 1179 | END SUBROUTINE fftradix |
---|
| 1180 | |
---|
| 1181 | END MODULE singleton |
---|