Changeset 1682 for palm/trunk/SOURCE/singleton.f90
- Timestamp:
- Oct 7, 2015 11:56:08 PM (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/singleton.f90
r1321 r1682 1 MODULE singleton 2 1 !> @file singleton.f90 3 2 !----------------------------------------------------------------------------- 4 3 ! Current revisions: 5 4 ! ----------------- 6 ! 5 ! Code annotations made doxygen readable 7 6 ! 8 7 ! Former revisions: … … 21 20 ! Description: 22 21 ! ------------ 23 ! Multivariate Fast Fourier Transform24 ! 25 ! Fortran 90 Implementation of Singleton's mixed-radix algorithm,26 ! RC Singleton, Stanford Research Institute, Sept. 1968.27 ! 28 ! Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and29 ! John Beale.30 ! 31 ! Fourier transforms can be computed either in place, using assumed size32 ! arguments, or by generic function, using assumed shape arguments.33 ! 34 ! 35 ! Public:36 ! 37 ! fftkind kind parameter of complex arguments38 ! and function results.39 ! 40 ! fft(array, dim, inv, stat) generic transform function41 ! COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array42 ! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim43 ! LOGICAL, INTENT(IN), OPTIONAL:: inv44 ! INTEGER, INTENT(OUT), OPTIONAL:: stat45 ! 46 ! fftn(array, shape, dim, inv, stat) in place transform subroutine47 ! COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array48 ! INTEGER, DIMENSION(:), INTENT(IN) :: shape49 ! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim50 ! LOGICAL, INTENT(IN), OPTIONAL:: inv51 ! INTEGER, INTENT(OUT), OPTIONAL:: stat52 ! 53 ! 54 ! Formal Parameters:55 ! 56 ! array The complex array to be transformed. array can be of arbitrary57 ! rank (i.e. up to seven).58 ! 59 ! shape With subroutine fftn, the shape of the array to be transformed60 ! has to be passed separately, since fftradix - the internal trans-61 ! formation routine - will treat array always as one dimensional.62 ! The product of elements in shape must be the number of63 ! elements in array.64 ! Although passing array with assumed shape would have been nicer,65 ! I prefered assumed size in order to prevent the compiler from66 ! using a copy-in-copy-out mechanism. That would generally be67 ! necessary with fftn passing array to fftradix and with fftn68 ! being prepared for accepting non consecutive array sections.69 ! Using assumed size, it's up to the user to pass an array argu-70 ! ment, that can be addressed as continous one dimensional array71 ! without copying. Otherwise, transformation will not really be72 ! performed in place.73 ! On the other hand, since the rank of array and the size of74 ! shape needn't match, fftn is appropriate for handling more than75 ! seven dimensions.76 ! As far as function fft is concerned all this doesn't matter,77 ! because the argument will be copied anyway. Thus no extra78 ! shape argument is needed for fft.79 ! 80 ! Optional Parameters:81 ! 82 ! dim One dimensional integer array, containing the dimensions to be83 ! transformed. Default is (/1,...,N/) with N being the rank of84 ! array, i.e. complete transform. dim can restrict transformation85 ! to a subset of available dimensions. Its size must not exceed the86 ! rank of array or the size of shape respectivly.87 ! 88 ! inv If .true., inverse transformation will be performed. Default is89 ! .false., i.e. forward transformation.90 ! 91 ! stat If present, a system dependent nonzero status value will be92 ! returned in stat, if allocation of temporary storage failed.93 ! 94 ! 95 ! Scaling:96 ! 97 ! Transformation results will always be scaled by the square root of the98 ! product of sizes of each dimension in dim. (See examples below)99 ! 100 ! 101 ! Examples:102 ! 103 ! Let A be a L*M*N three dimensional complex array. Then104 ! 105 ! result = fft(A)106 ! 107 ! will produce a three dimensional transform, scaled by sqrt(L*M*N), while108 ! 109 ! call fftn(A, SHAPE(A))110 ! 111 ! will do the same in place.112 ! 113 ! result = fft(A, dim=(/1,3/))114 ! 115 ! will transform with respect to the first and the third dimension, scaled116 ! by sqrt(L*N).117 ! 118 ! result = fft(fft(A), inv=.true.)119 ! 120 ! should (approximately) reproduce A.121 ! With B having the same shape as A122 ! 123 ! result = fft(fft(A) * CONJG(fft(B)), inv=.true.)124 ! 125 ! will correlate A and B.126 ! 127 ! 128 ! Remarks:129 ! 130 ! Following changes have been introduced with respect to fftn.c:131 ! - complex arguments and results are of type complex, rather than132 ! real an imaginary part separately.133 ! - increment parameter (magnitude of isign) has been dropped,134 ! inc is always one, direction of transform is given by inv.135 ! - maxf and maxp have been dropped. The amount of temporary storage136 ! needed is determined by the fftradix routine. Both fftn and fft137 ! can handle any size of array. (Maybe they take a lot of time and138 ! memory, but they will do it)139 ! 140 ! Redesigning fftradix in a way, that it handles assumed shape arrays141 ! would have been desirable. However, I found it rather hard to do this142 ! in an efficient way. Problems were:143 ! - to prevent stride multiplications when indexing arrays. At least our144 ! compiler was not clever enough to discover that in fact additions145 ! would do the job as well. On the other hand, I haven't been clever146 ! enough to find an implementation using array operations.147 ! - fftradix is rather large and different versions would be necessaray148 ! for each possible rank of array.149 ! Consequently, in place transformation still needs the argument stored150 ! in a consecutive bunch of memory and can't be performed on array151 ! sections like A(100:199:-3, 50:1020). Calling fftn with such sections152 ! will most probably imply copy-in-copy-out. However, the function fft153 ! works with everything it gets and should be convenient to use.154 ! 155 ! Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de>156 ! Restructured fftradix for better optimization. M. Steffens, 4 June 199722 !> Multivariate Fast Fourier Transform 23 !> 24 !> Fortran 90 Implementation of Singleton's mixed-radix algorithm, 25 !> RC Singleton, Stanford Research Institute, Sept. 1968. 26 !> 27 !> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and 28 !> John Beale. 29 !> 30 !> Fourier transforms can be computed either in place, using assumed size 31 !> arguments, or by generic function, using assumed shape arguments. 32 !> 33 !> 34 !> Public: 35 !> 36 !> fftkind kind parameter of complex arguments 37 !> and function results. 38 !> 39 !> fft(array, dim, inv, stat) generic transform function 40 !> COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array 41 !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim 42 !> LOGICAL, INTENT(IN), OPTIONAL:: inv 43 !> INTEGER, INTENT(OUT), OPTIONAL:: stat 44 !> 45 !> fftn(array, shape, dim, inv, stat) in place transform subroutine 46 !> COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array 47 !> INTEGER, DIMENSION(:), INTENT(IN) :: shape 48 !> INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim 49 !> LOGICAL, INTENT(IN), OPTIONAL:: inv 50 !> INTEGER, INTENT(OUT), OPTIONAL:: stat 51 !> 52 !> 53 !> Formal Parameters: 54 !> 55 !> array The complex array to be transformed. array can be of arbitrary 56 !> rank (i.e. up to seven). 57 !> 58 !> shape With subroutine fftn, the shape of the array to be transformed 59 !> has to be passed separately, since fftradix - the internal trans- 60 !> formation routine - will treat array always as one dimensional. 61 !> The product of elements in shape must be the number of 62 !> elements in array. 63 !> Although passing array with assumed shape would have been nicer, 64 !> I prefered assumed size in order to prevent the compiler from 65 !> using a copy-in-copy-out mechanism. That would generally be 66 !> necessary with fftn passing array to fftradix and with fftn 67 !> being prepared for accepting non consecutive array sections. 68 !> Using assumed size, it's up to the user to pass an array argu- 69 !> ment, that can be addressed as continous one dimensional array 70 !> without copying. Otherwise, transformation will not really be 71 !> performed in place. 72 !> On the other hand, since the rank of array and the size of 73 !> shape needn't match, fftn is appropriate for handling more than 74 !> seven dimensions. 75 !> As far as function fft is concerned all this doesn't matter, 76 !> because the argument will be copied anyway. Thus no extra 77 !> shape argument is needed for fft. 78 !> 79 !> Optional Parameters: 80 !> 81 !> dim One dimensional integer array, containing the dimensions to be 82 !> transformed. Default is (/1,...,N/) with N being the rank of 83 !> array, i.e. complete transform. dim can restrict transformation 84 !> to a subset of available dimensions. Its size must not exceed the 85 !> rank of array or the size of shape respectivly. 86 !> 87 !> inv If .true., inverse transformation will be performed. Default is 88 !> .false., i.e. forward transformation. 89 !> 90 !> stat If present, a system dependent nonzero status value will be 91 !> returned in stat, if allocation of temporary storage failed. 92 !> 93 !> 94 !> Scaling: 95 !> 96 !> Transformation results will always be scaled by the square root of the 97 !> product of sizes of each dimension in dim. (See examples below) 98 !> 99 !> 100 !> Examples: 101 !> 102 !> Let A be a L*M*N three dimensional complex array. Then 103 !> 104 !> result = fft(A) 105 !> 106 !> will produce a three dimensional transform, scaled by sqrt(L*M*N), while 107 !> 108 !> call fftn(A, SHAPE(A)) 109 !> 110 !> will do the same in place. 111 !> 112 !> result = fft(A, dim=(/1,3/)) 113 !> 114 !> will transform with respect to the first and the third dimension, scaled 115 !> by sqrt(L*N). 116 !> 117 !> result = fft(fft(A), inv=.true.) 118 !> 119 !> should (approximately) reproduce A. 120 !> With B having the same shape as A 121 !> 122 !> result = fft(fft(A) * CONJG(fft(B)), inv=.true.) 123 !> 124 !> will correlate A and B. 125 !> 126 !> 127 !> Remarks: 128 !> 129 !> Following changes have been introduced with respect to fftn.c: 130 !> - complex arguments and results are of type complex, rather than 131 !> real an imaginary part separately. 132 !> - increment parameter (magnitude of isign) has been dropped, 133 !> inc is always one, direction of transform is given by inv. 134 !> - maxf and maxp have been dropped. The amount of temporary storage 135 !> needed is determined by the fftradix routine. Both fftn and fft 136 !> can handle any size of array. (Maybe they take a lot of time and 137 !> memory, but they will do it) 138 !> 139 !> Redesigning fftradix in a way, that it handles assumed shape arrays 140 !> would have been desirable. However, I found it rather hard to do this 141 !> in an efficient way. Problems were: 142 !> - to prevent stride multiplications when indexing arrays. At least our 143 !> compiler was not clever enough to discover that in fact additions 144 !> would do the job as well. On the other hand, I haven't been clever 145 !> enough to find an implementation using array operations. 146 !> - fftradix is rather large and different versions would be necessaray 147 !> for each possible rank of array. 148 !> Consequently, in place transformation still needs the argument stored 149 !> in a consecutive bunch of memory and can't be performed on array 150 !> sections like A(100:199:-3, 50:1020). Calling fftn with such sections 151 !> will most probably imply copy-in-copy-out. However, the function fft 152 !> works with everything it gets and should be convenient to use. 153 !> 154 !> Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de> 155 !> Restructured fftradix for better optimization. M. Steffens, 4 June 1997 157 156 !----------------------------------------------------------------------------- 157 MODULE singleton 158 158 159 159 160 USE kinds … … 183 184 184 185 186 !------------------------------------------------------------------------------! 187 ! Description: 188 ! ------------ 189 !> @todo Missing function description. 190 !------------------------------------------------------------------------------! 185 191 FUNCTION fft1d(array, dim, inv, stat) RESULT(ft) 186 192 ! … … 207 213 208 214 215 !------------------------------------------------------------------------------! 216 ! Description: 217 ! ------------ 218 !> @todo Missing function description. 219 !------------------------------------------------------------------------------! 209 220 FUNCTION fft2d(array, dim, inv, stat) RESULT(ft) 210 221 ! … … 230 241 231 242 243 !------------------------------------------------------------------------------! 244 ! Description: 245 ! ------------ 246 !> @todo Missing function description. 247 !------------------------------------------------------------------------------! 232 248 FUNCTION fft3d(array, dim, inv, stat) RESULT(ft) 233 249 ! … … 255 271 256 272 273 !------------------------------------------------------------------------------! 274 ! Description: 275 ! ------------ 276 !> @todo Missing function description. 277 !------------------------------------------------------------------------------! 257 278 FUNCTION fft4d(array, dim, inv, stat) RESULT(ft) 258 279 ! … … 279 300 280 301 302 !------------------------------------------------------------------------------! 303 ! Description: 304 ! ------------ 305 !> @todo Missing function description. 306 !------------------------------------------------------------------------------! 281 307 FUNCTION fft5d(array, dim, inv, stat) RESULT(ft) 282 308 ! … … 305 331 306 332 333 !------------------------------------------------------------------------------! 334 ! Description: 335 ! ------------ 336 !> @todo Missing function description. 337 !------------------------------------------------------------------------------! 307 338 FUNCTION fft6d(array, dim, inv, stat) RESULT(ft) 308 339 ! … … 331 362 332 363 364 !------------------------------------------------------------------------------! 365 ! Description: 366 ! ------------ 367 !> @todo Missing function description. 368 !------------------------------------------------------------------------------! 333 369 FUNCTION fft7d(array, dim, inv, stat) RESULT(ft) 334 370 ! … … 357 393 358 394 395 !------------------------------------------------------------------------------! 396 ! Description: 397 ! ------------ 398 !> @todo Missing subroutine description. 399 !------------------------------------------------------------------------------! 359 400 SUBROUTINE fftn(array, shape, dim, inv, stat) 360 401 ! … … 410 451 411 452 453 !------------------------------------------------------------------------------! 454 ! Description: 455 ! ------------ 456 !> @todo Missing subroutine description. 457 !------------------------------------------------------------------------------! 412 458 SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat) 413 459 ! … … 472 518 473 519 520 !------------------------------------------------------------------------------! 521 ! Description: 522 ! ------------ 523 !> @todo Missing subroutine description. 524 !------------------------------------------------------------------------------! 474 525 SUBROUTINE factorize(npass, factor, nfactor, nsquare) 475 526 ! … … 536 587 537 588 589 !------------------------------------------------------------------------------! 590 ! Description: 591 ! ------------ 592 !> @todo Missing subroutine description. 593 !------------------------------------------------------------------------------! 538 594 SUBROUTINE transform(array, ntotal, npass, nspan, & 539 595 factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform … … 886 942 887 943 944 !------------------------------------------------------------------------------! 945 ! Description: 946 ! ------------ 947 !> @todo Missing subroutine description. 948 !------------------------------------------------------------------------------! 888 949 SUBROUTINE permute(array, ntotal, npass, nspan, & 889 950 factor, nfactor, nsquare, maxfactor, &
Note: See TracChangeset
for help on using the changeset viewer.