|
XBeach
|
00001 module math_tools 00002 ! MODULE CONTAINS: 00003 ! 1) Singleton fft transform 00004 ! 2) Hilbert transform 00005 ! 3) Matrix flip up to down and/or left to right 00006 ! 4) xerf 00007 ! 5) Random function 00008 00009 ! 00010 ! 1) Singleton fft transform: 00011 ! 00012 !----------------------------------------------------------------------------- 00013 ! Multivariate Fast Fourier Transform 00014 ! 00015 ! Fortran 90 Implementation of Singleton's mixed-radix algorithm, 00016 ! RC Singleton, Stanford Research Institute, Sept. 1968. 00017 ! 00018 ! Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and 00019 ! John Beale. 00020 ! 00021 ! Fourier transforms can be computed either in place, using assumed size 00022 ! arguments, or by generic function, using assumed shape arguments. 00023 ! 00024 ! 00025 ! Public: 00026 ! 00027 ! fftkind kind parameter of complex arguments 00028 ! and function results. 00029 ! 00030 ! fft(array, dim, inv, stat) generic transform function 00031 ! COMPLEX(fftkind), DIMENSION(:,...,:), INTENT(IN) :: array 00032 ! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim 00033 ! LOGICAL, INTENT(IN), OPTIONAL:: inv 00034 ! INTEGER, INTENT(OUT), OPTIONAL:: stat 00035 ! 00036 ! fftn(array, shape, dim, inv, stat) in place transform subroutine 00037 ! COMPLEX(fftkind), DIMENSION(*), INTENT(INOUT) :: array 00038 ! INTEGER, DIMENSION(:), INTENT(IN) :: shape 00039 ! INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL:: dim 00040 ! LOGICAL, INTENT(IN), OPTIONAL:: inv 00041 ! INTEGER, INTENT(OUT), OPTIONAL:: stat 00042 ! 00043 ! 00044 ! Formal Parameters: 00045 ! 00046 ! array The complex array to be transformed. array can be of arbitrary 00047 ! rank (i.e. up to seven). 00048 ! 00049 ! shape With subroutine fftn, the shape of the array to be transformed 00050 ! has to be passed separately, since fftradix - the internal trans- 00051 ! formation routine - will treat array always as one dimensional. 00052 ! The product of elements in shape must be the number of 00053 ! elements in array. 00054 ! Although passing array with assumed shape would have been nicer, 00055 ! I prefered assumed size in order to prevent the compiler from 00056 ! using a copy-in-copy-out mechanism. That would generally be 00057 ! necessary with fftn passing array to fftradix and with fftn 00058 ! being prepared for accepting non consecutive array sections. 00059 ! Using assumed size, it's up to the user to pass an array argu- 00060 ! ment, that can be addressed as continous one dimensional array 00061 ! without copying. Otherwise, transformation will not really be 00062 ! performed in place. 00063 ! On the other hand, since the rank of array and the size of 00064 ! shape needn't match, fftn is appropriate for handling more than 00065 ! seven dimensions. 00066 ! As far as function fft is concerned all this doesn't matter, 00067 ! because the argument will be copied anyway. Thus no extra 00068 ! shape argument is needed for fft. 00069 ! 00070 ! Optional Parameters: 00071 ! 00072 ! dim One dimensional integer array, containing the dimensions to be 00073 ! transformed. Default is (/1,...,N/) with N being the rank of 00074 ! array, i.e. complete transform. dim can restrict transformation 00075 ! to a subset of available dimensions. Its size must not exceed the 00076 ! rank of array or the size of shape respectivly. 00077 ! 00078 ! inv If .true., inverse transformation will be performed. Default is 00079 ! .false., i.e. forward transformation. 00080 ! 00081 ! stat If present, a system dependent nonzero status value will be 00082 ! returned in stat, if allocation of temporary storage failed. 00083 ! 00084 ! 00085 ! Scaling: 00086 ! 00087 ! Transformation results will always be scaled by the square root of the 00088 ! product of sizes of each dimension in dim. (See examples below) 00089 ! 00090 ! 00091 ! Examples: 00092 ! 00093 ! Let A be a L*M*N three dimensional complex array. Then 00094 ! 00095 ! result = fft(A) 00096 ! 00097 ! will produce a three dimensional transform, scaled by sqrt(L*M*N), while 00098 ! 00099 ! call fftn(A, SHAPE(A)) 00100 ! 00101 ! will do the same in place. 00102 ! 00103 ! result = fft(A, dim=(/1,3/)) 00104 ! 00105 ! will transform with respect to the first and the third dimension, scaled 00106 ! by sqrt(L*N). 00107 ! 00108 ! result = fft(fft(A), inv=.true.) 00109 ! 00110 ! should (approximately) reproduce A. 00111 ! With B having the same shape as A 00112 ! 00113 ! result = fft(fft(A) * CONJG(fft(B)), inv=.true.) 00114 ! 00115 ! will correlate A and B. 00116 ! 00117 ! 00118 ! Remarks: 00119 ! 00120 ! Following changes have been introduced with respect to fftn.c: 00121 ! - complex arguments and results are of type complex, rather than 00122 ! real an imaginary part separately. 00123 ! - increment parameter (magnitude of isign) has been dropped, 00124 ! inc is always one, direction of transform is given by inv. 00125 ! - maxf and maxp have been dropped. The amount of temporary storage 00126 ! needed is determined by the fftradix routine. Both fftn and fft 00127 ! can handle any size of array. (Maybe they take a lot of time and 00128 ! memory, but they will do it) 00129 ! 00130 ! Redesigning fftradix in a way, that it handles assumed shape arrays 00131 ! would have been desirable. However, I found it rather hard to do this 00132 ! in an efficient way. Problems were: 00133 ! - to prevent stride multiplications when indexing arrays. At least our 00134 ! compiler was not clever enough to discover that in fact additions 00135 ! would do the job as well. On the other hand, I haven't been clever 00136 ! enough to find an implementation using array operations. 00137 ! - fftradix is rather large and different versions would be necessaray 00138 ! for each possible rank of array. 00139 ! Consequently, in place transformation still needs the argument stored 00140 ! in a consecutive bunch of memory and can't be performed on array 00141 ! sections like A(100:199:-3, 50:1020). Calling fftn with such sections 00142 ! will most probably imply copy-in-copy-out. However, the function fft 00143 ! works with everything it gets and should be convenient to use. 00144 ! 00145 ! Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de> 00146 !----------------------------------------------------------------------------- 00147 implicit none 00148 save 00149 00150 private 00151 public:: fft, fftn, fftkind, realkind, flipa, flipv, hilbert, flipiv, xerf, random 00152 00153 integer, parameter:: fftkind = kind(0.0d0) !selected_real_kind(14,307) ! !--- adjust here for other precisions 00154 integer, parameter:: realkind = kind(0.0d0) 00155 real(fftkind), parameter:: sin60 = 0.86602540378443865_fftkind 00156 real(fftkind), parameter:: cos72 = 0.30901699437494742_fftkind 00157 real(fftkind), parameter:: sin72 = 0.95105651629515357_fftkind 00158 real(fftkind), parameter:: pi = 3.14159265358979323_fftkind 00159 00160 interface fft 00161 module procedure fft1d 00162 module procedure fft2d 00163 module procedure fft3d 00164 module procedure fft4d 00165 module procedure fft5d 00166 module procedure fft6d 00167 module procedure fft7d 00168 end interface fft 00169 00170 contains 00171 00172 function fft1d(array, dim, inv, stat) result(ft) 00173 ! wwvv not used: dim 00174 ! this whole function is not used, but i put dim in the fftn call 00175 !--- formal parameters 00176 complex(fftkind), dimension(:), intent(in) :: array 00177 integer, dimension(:), intent(in), optional:: dim 00178 logical, intent(in), optional:: inv 00179 integer, intent(out), optional:: stat 00180 !--- function result 00181 complex(fftkind), dimension(size(array, 1)):: ft 00182 !--- intrinsics used 00183 intrinsic size, shape 00184 ft = array 00185 call fftn(ft, shape(array), dim = dim, inv = inv, stat = stat) 00186 00187 end function fft1d 00188 00189 00190 function fft2d(array, dim, inv, stat) result(ft) 00191 !--- formal parameters 00192 complex(fftkind), dimension(:,:), intent(in) :: array 00193 integer, dimension(:), intent(in), optional:: dim 00194 logical, intent(in), optional:: inv 00195 integer, intent(out), optional:: stat 00196 !--- function result 00197 complex(fftkind), dimension(size(array, 1), size(array, 2)):: ft 00198 !--- intrinsics used 00199 intrinsic size, shape 00200 00201 ft = array 00202 call fftn(ft, shape(array), dim, inv, stat) 00203 00204 end function fft2d 00205 00206 00207 function fft3d(array, dim, inv, stat) result(ft) 00208 !--- formal parameters 00209 complex(fftkind), dimension(:,:,:), intent(in) :: array 00210 integer, dimension(:), intent(in), optional:: dim 00211 logical, intent(in), optional:: inv 00212 integer, intent(out), optional:: stat 00213 !--- function result 00214 complex(fftkind), 00215 dimension(size(array, 1), size(array, 2), size(array, 3)):: ft 00216 !--- intrinsics used 00217 intrinsic size, shape 00218 00219 ft = array 00220 call fftn(ft, shape(array), dim, inv, stat) 00221 00222 end function fft3d 00223 00224 00225 function fft4d(array, dim, inv, stat) result(ft) 00226 !--- formal parameters 00227 complex(fftkind), dimension(:,:,:,:), intent(in) :: array 00228 integer, dimension(:), intent(in), optional:: dim 00229 logical, intent(in), optional:: inv 00230 integer, intent(out), optional:: stat 00231 !--- function result 00232 complex(fftkind), dimension( & size(array, 1), size(array, 2), size(array, 3), size(array, 4)):: ft 00233 !--- intrinsics used 00234 intrinsic size, shape 00235 00236 ft = array 00237 call fftn(ft, shape(array), dim, inv, stat) 00238 00239 end function fft4d 00240 00241 00242 function fft5d(array, dim, inv, stat) result(ft) 00243 !--- formal parameters 00244 complex(fftkind), dimension(:,:,:,:,:), intent(in) :: array 00245 integer, dimension(:), intent(in), optional:: dim 00246 logical, intent(in), optional:: inv 00247 integer, intent(out), optional:: stat 00248 !--- function result 00249 complex(fftkind), dimension( & size(array, 1), size(array, 2), size(array, 3), size(array, 4), & size(array, 5)):: ft 00250 !--- intrinsics used 00251 intrinsic size, shape 00252 00253 ft = array 00254 call fftn(ft, shape(array), dim, inv, stat) 00255 00256 end function fft5d 00257 00258 00259 function fft6d(array, dim, inv, stat) result(ft) 00260 !--- formal parameters 00261 complex(fftkind), dimension(:,:,:,:,:,:), intent(in) :: array 00262 integer, dimension(:), intent(in), optional:: dim 00263 logical, intent(in), optional:: inv 00264 integer, intent(out), optional:: stat 00265 !--- function result 00266 complex(fftkind), dimension( & size(array, 1), size(array, 2), size(array, 3), size(array, 4), & size(array, 5), size(array, 6)):: ft 00267 !--- intrinsics used 00268 intrinsic size, shape 00269 00270 ft = array 00271 call fftn(ft, shape(array), dim, inv, stat) 00272 00273 end function fft6d 00274 00275 00276 function fft7d(array, dim, inv, stat) result(ft) 00277 !--- formal parameters 00278 complex(fftkind), dimension(:,:,:,:,:,:,:), intent(in) :: array 00279 integer, dimension(:), intent(in), optional:: dim 00280 logical, intent(in), optional:: inv 00281 integer, intent(out), optional:: stat 00282 !--- function result 00283 complex(fftkind), dimension( & size(array, 1), size(array, 2), size(array, 3), size(array, 4), & size(array, 5), size(array, 6), size(array, 7)):: ft 00284 !--- intrinsics used 00285 intrinsic size, shape 00286 00287 ft = array 00288 call fftn(ft, shape(array), dim, inv, stat) 00289 00290 end function fft7d 00291 00292 00293 subroutine fftn(array, shape, dim, inv, stat) 00294 !--- formal parameters 00295 complex(fftkind), dimension(*), intent(inout) :: array 00296 integer, dimension(:), intent(in) :: shape 00297 integer, dimension(:), intent(in), optional:: dim 00298 logical, intent(in), optional:: inv 00299 integer, intent(out), optional:: stat 00300 !--- local arrays 00301 integer, dimension(size(shape)):: d 00302 !--- local scalars 00303 logical :: inverse 00304 integer :: i, ndim, ntotal 00305 real(fftkind):: scale 00306 !--- intrinsics used 00307 intrinsic present, min, product, size, sqrt 00308 00309 i=0 ! wwvv to quiet the compiler 00310 !--- optional parameter settings 00311 if (present(inv)) then 00312 inverse = inv 00313 else 00314 inverse = .false. 00315 end if 00316 if (present(dim)) then 00317 ndim = min(size(dim), size(d)) 00318 d(1:ndim) = dim(1:ndim) 00319 else 00320 ndim = size(d) 00321 d = (/(i, i = 1, size(d))/) 00322 end if 00323 00324 ntotal = product(shape) 00325 scale = sqrt(1.0_fftkind / product(shape(d(1:ndim)))) 00326 ! forall (i = 1: ntotal) array(i) = array(i) * scale 00327 do i = 1, ntotal 00328 array(i) = array(i) * scale 00329 end do 00330 do i = 1, ndim 00331 call fftradix(array, ntotal, shape(d(i)), product(shape(1:d(i))), & 00332 inverse, stat) 00333 if (present(stat)) then 00334 if (stat /=0) return 00335 end if 00336 end do 00337 00338 end subroutine fftn 00339 00340 00341 subroutine fftradix(array, ntotal, npass, nspan, inv, stat) 00342 !--- formal parameters 00343 integer, intent(in) :: ntotal, npass, nspan 00344 complex(fftkind), dimension(*), intent(inout) :: array 00345 logical, intent(in) :: inv 00346 integer, intent(out), optional:: stat 00347 !--- local arrays 00348 integer, dimension(bit_size(0)) :: factor 00349 complex(fftkind), dimension(:), allocatable :: ctmp 00350 real(fftkind), dimension(:), allocatable :: sine, cosine 00351 integer, dimension(:), allocatable :: perm 00352 !--- local scalars 00353 integer :: ii, kspan, ispan 00354 integer :: j, jc, jf, jj, k, k1, k2, k3, k4, kk, kt, nn, ns, nt 00355 integer :: maxfactor, nfactor, nperm 00356 real(fftkind) :: s60, c72, s72, pi2 00357 real(fftkind) :: radf 00358 real(fftkind) :: c1, c2, c3, cd, ak 00359 real(fftkind) :: s1, s2, s3, sd 00360 complex(fftkind):: cc, cj, ck, cjp, cjm, ckp, ckm 00361 !--- intrinsics used 00362 intrinsic maxval, mod, present, ishft, bit_size, sin, cos, & 00363 cmplx, real, aimag 00364 00365 if (npass <= 1) return 00366 00367 c72 = cos72 00368 if (inv) then 00369 s72 = sin72 00370 s60 = sin60 00371 pi2 = pi 00372 else 00373 s72 = -sin72 00374 s60 = -sin60 00375 pi2 = -pi 00376 end if 00377 00378 nt = ntotal 00379 ns = nspan 00380 kspan = ns 00381 00382 nn = nt - 1 00383 jc = ns / npass 00384 radf = pi2 * jc 00385 pi2 = pi2 * 2.0_fftkind !-- use 2 PI from here on 00386 00387 call factorize 00388 00389 maxfactor = maxval(factor (:nfactor)) 00390 if (nfactor - ishft(kt, 1) > 0) then 00391 nperm = max(nfactor + 1, product(factor(kt+1: nfactor-kt)) - 1) 00392 else 00393 nperm = nfactor + 1 00394 end if 00395 00396 if (present(stat)) then 00397 allocate(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), stat=stat) 00398 if (stat /= 0) return 00399 call transform 00400 deallocate(sine, cosine, stat=stat) 00401 if (stat /= 0) return 00402 allocate(perm(nperm), stat=stat) 00403 if (stat /= 0) return 00404 call permute 00405 deallocate(perm, ctmp, stat=stat) 00406 if (stat /= 0) return 00407 else 00408 allocate(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor)) 00409 call transform 00410 deallocate(sine, cosine) 00411 allocate(perm(nperm)) 00412 call permute 00413 deallocate(perm, ctmp) 00414 end if 00415 00416 contains 00417 00418 subroutine factorize 00419 nfactor = 0 00420 k = npass 00421 do while (mod(k, 16) == 0) 00422 nfactor = nfactor + 1 00423 factor (nfactor) = 4 00424 k = k / 16 00425 end do 00426 j = 3 00427 jj = 9 00428 do 00429 do while (mod(k, jj) == 0) 00430 nfactor=nfactor + 1 00431 factor (nfactor) = j 00432 k = k / jj 00433 end do 00434 j = j + 2 00435 jj = j * j 00436 if (jj > k) exit 00437 end do 00438 if (k <= 4) then 00439 kt = nfactor 00440 factor (nfactor + 1) = k 00441 if (k /= 1) nfactor = nfactor + 1 00442 else 00443 if (k - ishft(k / 4, 2) == 0) then 00444 nfactor = nfactor + 1 00445 factor (nfactor) = 2 00446 k = k / 4 00447 end if 00448 kt = nfactor 00449 j = 2 00450 do 00451 if (mod(k, j) == 0) then 00452 nfactor = nfactor + 1 00453 factor (nfactor) = j 00454 k = k / j 00455 end if 00456 j = ishft((j + 1)/2, 1) + 1 00457 if (j > k) exit 00458 end do 00459 end if 00460 if (kt > 0) then 00461 j = kt 00462 do 00463 nfactor = nfactor + 1 00464 factor (nfactor) = factor (j) 00465 j = j - 1 00466 if (j==0) exit 00467 end do 00468 end if 00469 end subroutine factorize 00470 00471 00472 subroutine transform !-- compute fourier transform 00473 ii = 0 00474 jf = 0 00475 do 00476 sd = radf / kspan 00477 cd = sin(sd) 00478 cd = 2.0_fftkind * cd * cd 00479 sd = sin(sd + sd) 00480 kk = 1 00481 ii = ii + 1 00482 00483 select case (factor (ii)) 00484 case (2) 00485 00486 !-- transform for factor of 2 (including rotation factor) 00487 kspan = kspan / 2 00488 k1 = kspan + 2 00489 do 00490 do 00491 k2 = kk + kspan 00492 ck = array(k2) 00493 array(k2) = array(kk)-ck 00494 array(kk) = array(kk) + ck 00495 kk = k2 + kspan 00496 if (kk > nn) exit 00497 end do 00498 kk = kk - nn 00499 if (kk > jc) exit 00500 end do 00501 if (kk > kspan) return 00502 do 00503 c1 = 1.0_fftkind - cd 00504 s1 = sd 00505 do 00506 do 00507 do 00508 k2 = kk + kspan 00509 ck = array(kk) - array(k2) 00510 array(kk) = array(kk) + array(k2) 00511 array(k2) = ck * cmplx(c1, s1, kind=fftkind) 00512 kk = k2 + kspan 00513 if (kk >= nt) exit 00514 end do 00515 k2 = kk - nt 00516 c1 = -c1 00517 kk = k1 - k2 00518 if (kk <= k2) exit 00519 end do 00520 ak = c1 - (cd * c1 + sd * s1) 00521 s1 = sd * c1 - cd * s1 + s1 00522 c1 = 2.0_fftkind - (ak * ak + s1 * s1) 00523 s1 = s1 * c1 00524 c1 = c1 * ak 00525 kk = kk + jc 00526 if (kk >= k2) exit 00527 end do 00528 k1 = k1 + 1 + 1 00529 kk = (k1 - kspan) / 2 + jc 00530 if (kk > jc + jc) exit 00531 end do 00532 00533 00534 case (4) !-- transform for factor of 4 00535 ispan = kspan 00536 kspan = kspan / 4 00537 00538 do 00539 c1 = 1.0_fftkind 00540 s1 = 0.0_fftkind 00541 do 00542 do 00543 k1 = kk + kspan 00544 k2 = k1 + kspan 00545 k3 = k2 + kspan 00546 ckp = array(kk) + array(k2) 00547 ckm = array(kk) - array(k2) 00548 cjp = array(k1) + array(k3) 00549 cjm = array(k1) - array(k3) 00550 array(kk) = ckp + cjp 00551 cjp = ckp - cjp 00552 if (inv) then 00553 ckp = ckm + cmplx(-aimag(cjm), real(cjm), kind=fftkind) 00554 ckm = ckm + cmplx(aimag(cjm), -real(cjm), kind=fftkind) 00555 else 00556 ckp = ckm + cmplx(aimag(cjm), -real(cjm), kind=fftkind) 00557 ckm = ckm + cmplx(-aimag(cjm), real(cjm), kind=fftkind) 00558 end if 00559 !-- avoid useless multiplies 00560 if (s1 == 0.0_fftkind) then 00561 array(k1) = ckp 00562 array(k2) = cjp 00563 array(k3) = ckm 00564 else 00565 array(k1) = ckp * cmplx(c1, s1, kind=fftkind) 00566 array(k2) = cjp * cmplx(c2, s2, kind=fftkind) 00567 array(k3) = ckm * cmplx(c3, s3, kind=fftkind) 00568 end if 00569 kk = k3 + kspan 00570 if (kk > nt) exit 00571 end do 00572 00573 c2 = c1 - (cd * c1 + sd * s1) 00574 s1 = sd * c1 - cd * s1 + s1 00575 c1 = 2.0_fftkind - (c2 * c2 + s1 * s1) 00576 s1 = s1 * c1 00577 c1 = c1 * c2 00578 !-- values of c2, c3, s2, s3 that will get used next time 00579 c2 = c1 * c1 - s1 * s1 00580 s2 = 2.0_fftkind * c1 * s1 00581 c3 = c2 * c1 - s2 * s1 00582 s3 = c2 * s1 + s2 * c1 00583 kk = kk - nt + jc 00584 if (kk > kspan) exit 00585 end do 00586 kk = kk - kspan + 1 00587 if (kk > jc) exit 00588 end do 00589 if (kspan == jc) return 00590 00591 case default 00592 !-- transform for odd factors 00593 k = factor (ii) 00594 ispan = kspan 00595 kspan = kspan / k 00596 00597 00598 select case (k) 00599 case (3) !-- transform for factor of 3 (optional code) 00600 do 00601 do 00602 k1 = kk + kspan 00603 k2 = k1 + kspan 00604 ck = array(kk) 00605 cj = array(k1) + array(k2) 00606 array(kk) = ck + cj 00607 ck = ck - 0.5_fftkind * cj 00608 cj = (array(k1) - array(k2)) * s60 00609 array(k1) = ck + cmplx(-aimag(cj), real(cj), kind=fftkind) 00610 array(k2) = ck + cmplx(aimag(cj), -real(cj), kind=fftkind) 00611 kk = k2 + kspan 00612 if (kk >= nn) exit 00613 end do 00614 kk = kk - nn 00615 if (kk > kspan) exit 00616 end do 00617 00618 case (5) !-- transform for factor of 5 (optional code) 00619 c2 = c72 * c72 - s72 * s72 00620 s2 = 2.0_fftkind * c72 * s72 00621 do 00622 do 00623 k1 = kk + kspan 00624 k2 = k1 + kspan 00625 k3 = k2 + kspan 00626 k4 = k3 + kspan 00627 ckp = array(k1) + array(k4) 00628 ckm = array(k1) - array(k4) 00629 cjp = array(k2) + array(k3) 00630 cjm = array(k2) - array(k3) 00631 cc = array(kk) 00632 array(kk) = cc + ckp + cjp 00633 ck = ckp * c72 + cjp * c2 + cc 00634 cj = ckm * s72 + cjm * s2 00635 array(k1) = ck + cmplx(-aimag(cj), real(cj), kind=fftkind) 00636 array(k4) = ck + cmplx(aimag(cj), -real(cj), kind=fftkind) 00637 ck = ckp * c2 + cjp * c72 + cc 00638 cj = ckm * s2 - cjm * s72 00639 array(k2) = ck + cmplx(-aimag(cj), real(cj), kind=fftkind) 00640 array(k3) = ck + cmplx(aimag(cj), -real(cj), kind=fftkind) 00641 kk = k4 + kspan 00642 if (kk >= nn) exit 00643 end do 00644 kk = kk - nn 00645 if (kk > kspan) exit 00646 end do 00647 00648 case default 00649 00650 if (k /= jf) then 00651 jf = k 00652 s1 = pi2 / k 00653 c1 = cos(s1) 00654 s1 = sin(s1) 00655 cosine (jf) = 1.0_fftkind 00656 sine (jf) = 0.0_fftkind 00657 j = 1 00658 do 00659 cosine (j) = cosine (k) * c1 + sine (k) * s1 00660 sine (j) = cosine (k) * s1 - sine (k) * c1 00661 k = k-1 00662 cosine (k) = cosine (j) 00663 sine (k) = -sine (j) 00664 j = j + 1 00665 if (j >= k) exit 00666 end do 00667 end if 00668 do 00669 do 00670 k1 = kk 00671 k2 = kk + ispan 00672 cc = array(kk) 00673 ck = cc 00674 j = 1 00675 k1 = k1 + kspan 00676 do 00677 k2 = k2 - kspan 00678 j = j + 1 00679 ctmp(j) = array(k1) + array(k2) 00680 ck = ck + ctmp(j) 00681 j = j + 1 00682 ctmp(j) = array(k1) - array(k2) 00683 k1 = k1 + kspan 00684 if (k1 >= k2) exit 00685 end do 00686 array(kk) = ck 00687 k1 = kk 00688 k2 = kk + ispan 00689 j = 1 00690 do 00691 k1 = k1 + kspan 00692 k2 = k2 - kspan 00693 jj = j 00694 ck = cc 00695 cj = (0.0_fftkind, 0.0_fftkind) 00696 k = 1 00697 do 00698 k = k + 1 00699 ck = ck + ctmp(k) * cosine (jj) 00700 k = k + 1 00701 cj = cj + ctmp(k) * sine (jj) 00702 jj = jj + j 00703 if (jj > jf) jj = jj - jf 00704 if (k >= jf) exit 00705 end do 00706 k = jf - j 00707 array(k1) = ck + cmplx(-aimag(cj), real(cj), kind=fftkind) 00708 array(k2) = ck + cmplx(aimag(cj), -real(cj), kind=fftkind) 00709 j = j + 1 00710 if (j >= k) exit 00711 end do 00712 kk = kk + ispan 00713 if (kk > nn) exit 00714 end do 00715 kk = kk - nn 00716 if (kk > kspan) exit 00717 end do 00718 00719 end select 00720 !-- multiply by rotation factor (except for factors of 2 and 4) 00721 if (ii == nfactor) return 00722 kk = jc + 1 00723 do 00724 c2 = 1.0_fftkind - cd 00725 s1 = sd 00726 do 00727 c1 = c2 00728 s2 = s1 00729 kk = kk + kspan 00730 do 00731 do 00732 array(kk) = cmplx(c2, s2, kind=fftkind) * array(kk) 00733 kk = kk + ispan 00734 if (kk > nt) exit 00735 end do 00736 ak = s1 * s2 00737 s2 = s1 * c2 + c1 * s2 00738 c2 = c1 * c2 - ak 00739 kk = kk - nt + kspan 00740 if (kk > ispan) exit 00741 end do 00742 c2 = c1 - (cd * c1 + sd * s1) 00743 s1 = s1 + sd * c1 - cd * s1 00744 c1 = 2.0_fftkind - (c2 * c2 + s1 * s1) 00745 s1 = s1 * c1 00746 c2 = c2 * c1 00747 kk = kk - ispan + jc 00748 if (kk > kspan) exit 00749 end do 00750 kk = kk - kspan + jc + 1 00751 if (kk > jc + jc) exit 00752 end do 00753 00754 end select 00755 end do 00756 end subroutine transform 00757 00758 00759 subroutine permute 00760 !-- permute the results to normal order---done in two stages 00761 !-- permutation for square factors of n 00762 perm (1) = ns 00763 if (kt > 0) then 00764 k = kt + kt + 1 00765 if (nfactor < k) k = k - 1 00766 j = 1 00767 perm (k + 1) = jc 00768 do 00769 perm (j + 1) = perm (j) / factor (j) 00770 perm (k) = perm (k + 1) * factor (j) 00771 j = j + 1 00772 k = k - 1 00773 if (j >= k) exit 00774 end do 00775 k3 = perm (k + 1) 00776 kspan = perm (2) 00777 kk = jc + 1 00778 k2 = kspan + 1 00779 j = 1 00780 00781 if (npass /= ntotal) then 00782 permute_multi: do 00783 do 00784 do 00785 k = kk + jc 00786 do 00787 !-- swap array(kk) <> array(k2) 00788 ck = array(kk) 00789 array(kk) = array(k2) 00790 array(k2) = ck 00791 kk = kk + 1 00792 k2 = k2 + 1 00793 if (kk >= k) exit 00794 end do 00795 kk = kk + ns - jc 00796 k2 = k2 + ns - jc 00797 if (kk >= nt) exit 00798 end do 00799 kk = kk - nt + jc 00800 k2 = k2 - nt + kspan 00801 if (k2 >= ns) exit 00802 end do 00803 do 00804 do 00805 k2 = k2 - perm (j) 00806 j = j + 1 00807 k2 = perm (j + 1) + k2 00808 if (k2 <= perm (j)) exit 00809 end do 00810 j = 1 00811 do 00812 if (kk < k2) cycle permute_multi 00813 kk = kk + jc 00814 k2 = k2 + kspan 00815 if (k2 >= ns) exit 00816 end do 00817 if (kk >= ns) exit 00818 end do 00819 exit 00820 end do permute_multi 00821 else 00822 permute_single: do 00823 do 00824 !-- swap array(kk) <> array(k2) 00825 ck = array(kk) 00826 array(kk) = array(k2) 00827 array(k2) = ck 00828 kk = kk + 1 00829 k2 = k2 + kspan 00830 if (k2 >= ns) exit 00831 end do 00832 do 00833 do 00834 k2 = k2 - perm (j) 00835 j = j + 1 00836 k2 = perm (j + 1) + k2 00837 if (k2 <= perm (j)) exit 00838 end do 00839 j = 1 00840 do 00841 if (kk < k2) cycle permute_single 00842 kk = kk + 1 00843 k2 = k2 + kspan 00844 if (k2 >= ns) exit 00845 end do 00846 if (kk >= ns) exit 00847 end do 00848 exit 00849 end do permute_single 00850 end if 00851 jc = k3 00852 end if 00853 00854 if (ishft(kt, 1) + 1 >= nfactor) return 00855 00856 ispan = perm (kt + 1) 00857 !-- permutation for square-free factors of n 00858 j = nfactor - kt 00859 factor (j + 1) = 1 00860 do 00861 factor(j) = factor(j) * factor(j+1) 00862 j = j - 1 00863 if (j == kt) exit 00864 end do 00865 kt = kt + 1 00866 nn = factor(kt) - 1 00867 j = 0 00868 jj = 0 00869 do 00870 k = kt + 1 00871 k2 = factor(kt) 00872 kk = factor(k) 00873 j = j + 1 00874 if (j > nn) exit !-- exit infinite loop 00875 jj = jj + kk 00876 do while (jj >= k2) 00877 jj = jj - k2 00878 k2 = kk 00879 k = k + 1 00880 kk = factor(k) 00881 jj = jj + kk 00882 end do 00883 perm (j) = jj 00884 end do 00885 !-- determine the permutation cycles of length greater than 1 00886 j = 0 00887 do 00888 do 00889 j = j + 1 00890 kk = perm(j) 00891 if (kk >= 0) exit 00892 end do 00893 if (kk /= j) then 00894 do 00895 k = kk 00896 kk = perm (k) 00897 perm (k) = -kk 00898 if (kk == j) exit 00899 end do 00900 k3 = kk 00901 else 00902 perm (j) = -j 00903 if (j == nn) exit !-- exit infinite loop 00904 end if 00905 end do 00906 !-- reorder a and b, following the permutation cycles 00907 do 00908 j = k3 + 1 00909 nt = nt - ispan 00910 ii = nt - 1 + 1 00911 if (nt < 0) exit !-- exit infinite loop 00912 do 00913 do 00914 j = j-1 00915 if (perm(j) >= 0) exit 00916 end do 00917 jj = jc 00918 do 00919 kspan = jj 00920 if (jj > maxfactor) kspan = maxfactor 00921 jj = jj - kspan 00922 k = perm(j) 00923 kk = jc * k + ii + jj 00924 k1 = kk + kspan 00925 k2 = 0 00926 do 00927 k2 = k2 + 1 00928 ctmp(k2) = array(k1) 00929 k1 = k1 - 1 00930 if (k1 == kk) exit 00931 end do 00932 do 00933 k1 = kk + kspan 00934 k2 = k1 - jc * (k + perm(k)) 00935 k = -perm(k) 00936 do 00937 array(k1) = array(k2) 00938 k1 = k1 - 1 00939 k2 = k2 - 1 00940 if (k1 == kk) exit 00941 end do 00942 kk = k2 00943 if (k == j) exit 00944 end do 00945 k1 = kk + kspan 00946 k2 = 0 00947 do 00948 k2 = k2 + 1 00949 array(k1) = ctmp(k2) 00950 k1 = k1 - 1 00951 if (k1 == kk) exit 00952 end do 00953 if (jj == 0) exit 00954 end do 00955 if (j == 1) exit 00956 end do 00957 end do 00958 00959 end subroutine permute 00960 00961 end subroutine fftradix 00962 00963 ! 00964 ! 2) Hilbert transform 00965 ! 00966 00967 !%HILBERT Hilbert transform. 00968 00969 subroutine hilbert(xi,m) 00970 00971 ! use math_tools 00972 00973 implicit none 00974 00975 integer :: m,n 00976 complex(fftkind),dimension(m) :: xi 00977 complex(fftkind),dimension(:),allocatable :: x 00978 integer,dimension(:),allocatable :: h 00979 real :: p 00980 real :: two, logtwo 00981 00982 ! Avoid a bug in gfortran...http://stackoverflow.com/questions/10673701/can-i-call-the-fortran-log-function-with-a-number 00983 two = 2.0 00984 logtwo = log(two) 00985 00986 p=log(real(m))/logtwo 00987 n=ceiling(p) 00988 n=2**n 00989 00990 allocate(x(n)) 00991 x=0 00992 x(1:m)=real(xi) 00993 00994 x=fft(x, inv=.false.) 00995 x=x*sqrt(real(n)) ! Scale factor 00996 00997 allocate(h(n)) 00998 00999 h(1)=1 01000 h(2:(n/2))=2 01001 h((n/2)+1)=1 01002 h((n/2)+2:n)=0 01003 01004 x=x*h 01005 01006 deallocate(h) 01007 01008 x=fft(x,inv=.true.) 01009 x=x/sqrt(real(n)) ! Scale factor 01010 01011 xi=x(1:m) 01012 deallocate(x) 01013 01014 return 01015 01016 end subroutine hilbert 01017 01018 ! 01019 ! 3) Matrix flip up to down and/or left to right 01020 ! 01021 subroutine flipv(x,M) 01022 01023 implicit none 01024 01025 ! x is input vector 01026 ! M is size vector 01027 01028 integer :: M, i 01029 real*8, dimension(M) :: x 01030 real*8, dimension(:),allocatable :: temp1, temp3 01031 integer,dimension(:),allocatable :: temp2 01032 01033 allocate(temp1(M)) 01034 allocate(temp2(M)) 01035 allocate(temp3(M)) 01036 01037 i=0 ! wwvv to quiet the compiler 01038 01039 temp1=0 01040 temp2=0 01041 temp3=0 01042 01043 temp1=x 01044 temp2=(/(i,i=0,M-1)/) 01045 temp3=temp1(M-temp2) 01046 x=temp3 01047 01048 deallocate(temp1) 01049 deallocate(temp2) 01050 deallocate(temp3) 01051 01052 return 01053 01054 end subroutine flipv 01055 01056 subroutine flipiv(x,M) 01057 01058 implicit none 01059 01060 ! x is input vector 01061 ! M is size vector 01062 01063 integer :: M, i 01064 complex(fftkind),dimension(M) :: x 01065 real(realkind),allocatable,dimension(:) :: temp1r, temp3r, temp1i, temp3i 01066 complex(realkind), parameter :: compi=(0.,1.) 01067 allocate(temp1r(M)) 01068 allocate(temp3r(M)) 01069 allocate(temp1i(M)) 01070 allocate(temp3i(M)) 01071 temp1r=0.d0 01072 temp1i=0.d0 01073 temp3r=0.d0 01074 temp3i=0.d0 01075 temp1r=real(x) 01076 temp1i=aimag(x) 01077 do i=0,M-1 01078 temp3r(i+1)=temp1r(M-i) 01079 temp3i(i+1)=temp1i(M-i) 01080 enddo 01081 01082 x=temp3r+compi*temp3i 01083 deallocate(temp1r) 01084 deallocate(temp3r) 01085 deallocate(temp1i) 01086 deallocate(temp3i) 01087 01088 return 01089 01090 end subroutine flipiv 01091 01092 subroutine flipa(x,M1,M2,flip) 01093 use xmpi_module 01094 01095 implicit none 01096 01097 ! x is 2D array 01098 ! flip (1 or 2) is flip rows(1) or columns(2) 01099 01100 integer :: M1, M2, flip, ii 01101 real*8, dimension(M1,M2) :: x 01102 real*8, dimension(:,:),allocatable :: temp1a, temp3a 01103 01104 allocate(temp1a(M1,M2)) 01105 allocate(temp3a(M1,M2)) 01106 01107 temp1a=0 01108 temp3a=0 01109 01110 if (flip==1) then ! flip rows 01111 temp1a=x 01112 do ii=1,M1 01113 temp3a(:,ii)=temp1a(:,M1+1-ii) 01114 end do 01115 x=temp3a 01116 else if (flip==2) then ! flip columns 01117 temp1a=x 01118 do ii=1,M2 01119 temp3a(:,ii)=temp1a(:,M2+1-ii) 01120 end do 01121 x=temp3a 01122 else 01123 write(*,*) 'Catastrophe in flipa: flip must be set to 1 or 2' 01124 call halt_program 01125 end if 01126 01127 deallocate(temp1a) 01128 deallocate(temp3a) 01129 01130 return 01131 01132 end subroutine flipa 01133 01134 function xerf(x) result (y) 01135 01136 implicit none 01137 01138 integer :: i 01139 real*8, dimension(:) :: x 01140 real*8, dimension(size(x)) :: w,y 01141 integer :: k 01142 real*8 :: t 01143 real*8, dimension(0:64) :: a,b 01144 01145 ! based on derf.f from http://www.kurims.kyoto-u.ac.jp/~ooura/ 01146 01147 data (a(i), i = 0, 12) / & 01148 0.00000000005958930743d0, -0.00000000113739022964d0, & 01149 0.00000001466005199839d0, -0.00000016350354461960d0, & 01150 0.00000164610044809620d0, -0.00001492559551950604d0, & 01151 0.00012055331122299265d0, -0.00085483269811296660d0, & 01152 0.00522397762482322257d0, -0.02686617064507733420d0, & 01153 0.11283791670954881569d0, -0.37612638903183748117d0, & 01154 1.12837916709551257377d0 / 01155 data (a(i), i = 13, 25) / & 01156 0.00000000002372510631d0, -0.00000000045493253732d0, & 01157 0.00000000590362766598d0, -0.00000006642090827576d0, & 01158 0.00000067595634268133d0, -0.00000621188515924000d0, & 01159 0.00005103883009709690d0, -0.00037015410692956173d0, & 01160 0.00233307631218880978d0, -0.01254988477182192210d0, & 01161 0.05657061146827041994d0, -0.21379664776456006580d0, & 01162 0.84270079294971486929d0 / 01163 data (a(i), i = 26, 38) / & 01164 0.00000000000949905026d0, -0.00000000018310229805d0, & 01165 0.00000000239463074000d0, -0.00000002721444369609d0, & 01166 0.00000028045522331686d0, -0.00000261830022482897d0, & 01167 0.00002195455056768781d0, -0.00016358986921372656d0, & 01168 0.00107052153564110318d0, -0.00608284718113590151d0, & 01169 0.02986978465246258244d0, -0.13055593046562267625d0, & 01170 0.67493323603965504676d0 / 01171 data (a(i), i = 39, 51) / & 01172 0.00000000000382722073d0, -0.00000000007421598602d0, & 01173 0.00000000097930574080d0, -0.00000001126008898854d0, & 01174 0.00000011775134830784d0, -0.00000111992758382650d0, & 01175 0.00000962023443095201d0, -0.00007404402135070773d0, & 01176 0.00050689993654144881d0, -0.00307553051439272889d0, & 01177 0.01668977892553165586d0, -0.08548534594781312114d0, & 01178 0.56909076642393639985d0 / 01179 data (a(i), i = 52, 64) / & 01180 0.00000000000155296588d0, -0.00000000003032205868d0, & 01181 0.00000000040424830707d0, -0.00000000471135111493d0, & 01182 0.00000005011915876293d0, -0.00000048722516178974d0, & 01183 0.00000430683284629395d0, -0.00003445026145385764d0, & 01184 0.00024879276133931664d0, -0.00162940941748079288d0, & 01185 0.00988786373932350462d0, -0.05962426839442303805d0, & 01186 0.49766113250947636708d0 / 01187 data (b(i), i = 0, 12) / & 01188 -0.00000000029734388465d0, 0.00000000269776334046d0, & 01189 -0.00000000640788827665d0, -0.00000001667820132100d0, & 01190 -0.00000021854388148686d0, 0.00000266246030457984d0, & 01191 0.00001612722157047886d0, -0.00025616361025506629d0, & 01192 0.00015380842432375365d0, 0.00815533022524927908d0, & 01193 -0.01402283663896319337d0, -0.19746892495383021487d0, & 01194 0.71511720328842845913d0 / 01195 data (b(i), i = 13, 25) / & 01196 -0.00000000001951073787d0, -0.00000000032302692214d0, & 01197 0.00000000522461866919d0, 0.00000000342940918551d0, & 01198 -0.00000035772874310272d0, 0.00000019999935792654d0, & 01199 0.00002687044575042908d0, -0.00011843240273775776d0, & 01200 -0.00080991728956032271d0, 0.00661062970502241174d0, & 01201 0.00909530922354827295d0, -0.20160072778491013140d0, & 01202 0.51169696718727644908d0 / 01203 data (b(i), i = 26, 38) / & 01204 0.00000000003147682272d0, -0.00000000048465972408d0, & 01205 0.00000000063675740242d0, 0.00000003377623323271d0, & 01206 -0.00000015451139637086d0, -0.00000203340624738438d0, & 01207 0.00001947204525295057d0, 0.00002854147231653228d0, & 01208 -0.00101565063152200272d0, 0.00271187003520095655d0, & 01209 0.02328095035422810727d0, -0.16725021123116877197d0, & 01210 0.32490054966649436974d0 / 01211 data (b(i), i = 39, 51) / & 01212 0.00000000002319363370d0, -0.00000000006303206648d0, & 01213 -0.00000000264888267434d0, 0.00000002050708040581d0, & 01214 0.00000011371857327578d0, -0.00000211211337219663d0, & 01215 0.00000368797328322935d0, 0.00009823686253424796d0, & 01216 -0.00065860243990455368d0, -0.00075285814895230877d0, & 01217 0.02585434424202960464d0, -0.11637092784486193258d0, & 01218 0.18267336775296612024d0 / 01219 data (b(i), i = 52, 64) / & 01220 -0.00000000000367789363d0, 0.00000000020876046746d0, & 01221 -0.00000000193319027226d0, -0.00000000435953392472d0, & 01222 0.00000018006992266137d0, -0.00000078441223763969d0, & 01223 -0.00000675407647949153d0, 0.00008428418334440096d0, & 01224 -0.00017604388937031815d0, -0.00239729611435071610d0, & 01225 0.02064129023876022970d0, -0.06905562880005864105d0, & 01226 0.09084526782065478489d0 / 01227 01228 w = abs(x) 01229 01230 do i = 1,size(x) 01231 if (w(i) .lt. 2.2d0) then 01232 t = w(i) * w(i) 01233 k = int(t) 01234 t = t - k 01235 k = k * 13 01236 y(i) = ((((((((((((a(k) * t + a(k + 1)) * t + & 01237 a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + & 01238 a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + & 01239 a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + & 01240 a(k + 11)) * t + a(k + 12)) * w(i) 01241 else if (w(i) .lt. 6.9d0) then 01242 k = int(w(i)) 01243 t = w(i) - k 01244 k = 13 * (k - 2) 01245 y(i) = (((((((((((b(k) * t + b(k + 1)) * t + & 01246 b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & 01247 b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & 01248 b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & 01249 b(k + 11)) * t + b(k + 12) 01250 y(i) = 1 - y(i)**16 01251 else 01252 y(i) = 1 01253 end if 01254 01255 if (x(i) .lt. 0) y(i) = -y(i) 01256 enddo 01257 01258 end function xerf 01259 01260 01261 real*8 function random(j) 01262 implicit none 01263 ! 01264 ! From Communications of the ACM, Vol 31 Oct 1988 number 10 01265 ! pp 1192..1201 01266 ! Stephen K. Park and Keith W. Miller 01267 ! 01268 ! Input parameter j: 01269 ! j .ne. 0: random sequence is initialised, using iabs(j) as seed 01270 ! first random number of new series is returned 01271 ! j .eq. 0: function returns next random number 01272 ! 01273 ! Not thread-safe because of seed 01274 ! 01275 01276 integer a,m,q,r,lo,hi,test,seed,j 01277 real*8 rm 01278 parameter(a=16807,m=2147483647,q=127773,r=2836,rm=m) 01279 save seed 01280 data seed/1/ 01281 if (j.ne.0) then 01282 seed=mod(iabs(j),m) 01283 ! Seed may not be zero for random number generation 01284 if(seed==0) then 01285 seed=1 01286 endif 01287 endif 01288 hi=seed/q 01289 lo=mod(seed,q) 01290 test=a*lo-r*hi 01291 if (test .gt. 0) then 01292 seed=test 01293 else 01294 seed=test+m 01295 endif 01296 random=seed/rm 01297 return 01298 end function random 01299 01300 end module math_tools 01301
1.7.6.1
