Changeset 1320 for palm/trunk/SOURCE/singleton.f90
- Timestamp:
- Mar 20, 2014 8:40:49 AM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/singleton.f90 (modified) (34 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/singleton.f90
r484 r1320 4 4 ! Current revisions: 5 5 ! ----------------- 6 ! 6 ! kind-parameters added to all INTEGER and REAL declaration statements, 7 ! kinds are defined in new module kinds, 8 ! revision history before 2012 removed, 7 9 ! 8 10 ! Former revisions: 9 11 ! ----------------- 10 ! $Id$11 ! RCS Log replace by Id keyword, revision history cleaned up12 !13 ! Revision 1.2 2004/04/30 12:52:09 raasch14 ! Shape of arrays is explicitly stored in ishape and handled to the15 ! fft-routines instead of the shape-function (due to a compiler error on16 ! decalpha)17 !18 12 ! Revision 1.1 2002/05/02 18:56:59 raasch 19 13 ! Initial revision … … 158 152 !----------------------------------------------------------------------------- 159 153 154 USE kinds 155 160 156 IMPLICIT NONE 161 157 162 158 PRIVATE 163 PUBLIC:: fft, fftn, fftkind 164 165 INTEGER, PARAMETER:: fftkind = KIND(0.0) ! adjust here for other precisions 166 167 REAL(fftkind), PARAMETER:: sin60 = 0.86602540378443865_fftkind 168 REAL(fftkind), PARAMETER:: cos72 = 0.30901699437494742_fftkind 169 REAL(fftkind), PARAMETER:: sin72 = 0.95105651629515357_fftkind 170 REAL(fftkind), PARAMETER:: pi = 3.14159265358979323_fftkind 159 PUBLIC:: fft, fftn 160 161 REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp 162 REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp 163 REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp 164 REAL(wp), PARAMETER:: pi = 3.14159265358979323_wp 171 165 172 166 INTERFACE fft … … 187 181 ! 188 182 !-- Formal parameters 189 COMPLEX( fftkind),DIMENSION(:), INTENT(IN) :: array190 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim191 LOGICAL, INTENT(IN), OPTIONAL:: inv192 INTEGER, INTENT(OUT), OPTIONAL:: stat183 COMPLEX(wp), DIMENSION(:), INTENT(IN) :: array 184 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 185 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 186 LOGICAL, INTENT(IN), OPTIONAL:: inv 193 187 ! 194 188 !-- Function result 195 COMPLEX( fftkind), DIMENSION(SIZE(array, 1)):: ft196 197 INTEGER 189 COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft 190 191 INTEGER(iwp):: ishape(1) 198 192 199 193 ! … … 211 205 ! 212 206 !-- Formal parameters 213 COMPLEX( fftkind),DIMENSION(:,:), INTENT(IN) :: array214 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim215 LOGICAL, INTENT(IN), OPTIONAL:: inv216 INTEGER, INTENT(OUT), OPTIONAL:: stat207 COMPLEX(wp), DIMENSION(:,:), INTENT(IN) :: array 208 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 209 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 210 LOGICAL, INTENT(IN), OPTIONAL:: inv 217 211 ! 218 212 !-- Function result 219 COMPLEX( fftkind), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft220 221 INTEGER :: ishape(2)213 COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft 214 215 INTEGER(iwp) :: ishape(2) 222 216 ! 223 217 !-- Intrinsics used … … 234 228 ! 235 229 !-- Formal parameters 236 COMPLEX( fftkind),DIMENSION(:,:,:), INTENT(IN) :: array237 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim238 LOGICAL, INTENT(IN), OPTIONAL:: inv239 INTEGER, INTENT(OUT), OPTIONAL:: stat230 COMPLEX(wp), DIMENSION(:,:,:), INTENT(IN) :: array 231 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 232 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 233 LOGICAL, INTENT(IN), OPTIONAL:: inv 240 234 ! 241 235 !-- Function result 242 COMPLEX( fftkind), &236 COMPLEX(wp), & 243 237 DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft 244 238 245 INTEGER :: ishape(3)239 INTEGER(iwp) :: ishape(3) 246 240 247 241 ! … … 259 253 ! 260 254 !-- Formal parameters 261 COMPLEX( fftkind),DIMENSION(:,:,:,:), INTENT(IN) :: array262 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim263 LOGICAL, INTENT(IN), OPTIONAL:: inv264 INTEGER, INTENT(OUT), OPTIONAL:: stat255 COMPLEX(wp), DIMENSION(:,:,:,:), INTENT(IN) :: array 256 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 257 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 258 LOGICAL, INTENT(IN), OPTIONAL:: inv 265 259 ! 266 260 !-- Function result 267 COMPLEX( fftkind), DIMENSION( &261 COMPLEX(wp), DIMENSION( & 268 262 SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft 269 263 270 INTEGER :: ishape(4)264 INTEGER(iwp) :: ishape(4) 271 265 ! 272 266 !-- Intrinsics used … … 283 277 ! 284 278 !-- Formal parameters 285 COMPLEX( fftkind),DIMENSION(:,:,:,:,:), INTENT(IN) :: array286 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim287 LOGICAL, INTENT(IN), OPTIONAL:: inv288 INTEGER, INTENT(OUT), OPTIONAL:: stat279 COMPLEX(wp), DIMENSION(:,:,:,:,:), INTENT(IN) :: array 280 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 281 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 282 LOGICAL, INTENT(IN), OPTIONAL:: inv 289 283 ! 290 284 !-- Function result 291 COMPLEX( fftkind), DIMENSION( &285 COMPLEX(wp), DIMENSION( & 292 286 SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & 293 287 SIZE(array, 5)):: ft 294 288 295 INTEGER :: ishape(5)289 INTEGER(iwp) :: ishape(5) 296 290 297 291 ! … … 309 303 ! 310 304 !-- Formal parameters 311 COMPLEX( fftkind),DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array312 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim313 LOGICAL, INTENT(IN), OPTIONAL:: inv314 INTEGER, INTENT(OUT), OPTIONAL:: stat305 COMPLEX(wp), DIMENSION(:,:,:,:,:,:), INTENT(IN) :: array 306 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 307 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 308 LOGICAL, INTENT(IN), OPTIONAL:: inv 315 309 ! 316 310 !-- Function result 317 COMPLEX( fftkind), DIMENSION( &311 COMPLEX(wp), DIMENSION( & 318 312 SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & 319 313 SIZE(array, 5), SIZE(array, 6)):: ft 320 314 321 INTEGER :: ishape(6)315 INTEGER(iwp) :: ishape(6) 322 316 323 317 ! … … 335 329 ! 336 330 !-- Formal parameters 337 COMPLEX( fftkind), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array338 INTEGER , DIMENSION(:),INTENT(IN), OPTIONAL:: dim339 LOGICAL, INTENT(IN), OPTIONAL:: inv340 INTEGER, INTENT(OUT), OPTIONAL:: stat331 COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: array 332 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 333 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 334 LOGICAL, INTENT(IN), OPTIONAL:: inv 341 335 ! 342 336 !-- Function result 343 COMPLEX( fftkind), DIMENSION( &337 COMPLEX(wp), DIMENSION( & 344 338 SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), & 345 339 SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft 346 340 347 INTEGER :: ishape(7)341 INTEGER(iwp) :: ishape(7) 348 342 349 343 ! … … 361 355 ! 362 356 !-- Formal parameters 363 COMPLEX( fftkind),DIMENSION(*), INTENT(INOUT) :: array364 INTEGER ,DIMENSION(:), INTENT(IN) :: shape365 INTEGER ,DIMENSION(:), INTENT(IN), OPTIONAL:: dim366 LOGICAL, INTENT(IN), OPTIONAL:: inv367 INTEGER, INTENT(OUT), OPTIONAL:: stat357 COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array 358 INTEGER(iwp), DIMENSION(:), INTENT(IN) :: shape 359 INTEGER(iwp), DIMENSION(:), INTENT(IN), OPTIONAL:: dim 360 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 361 LOGICAL, INTENT(IN), OPTIONAL:: inv 368 362 ! 369 363 !-- Local arrays 370 INTEGER , DIMENSION(SIZE(shape)):: d364 INTEGER(iwp), DIMENSION(SIZE(shape)):: d 371 365 ! 372 366 !-- Local scalars 373 367 LOGICAL :: inverse 374 INTEGER 375 REAL( fftkind):: scale368 INTEGER(iwp) :: i, ndim, ntotal 369 REAL(wp):: scale 376 370 ! 377 371 !-- Intrinsics used … … 394 388 395 389 ntotal = PRODUCT(shape) 396 scale = SQRT(1.0_ fftkind/ PRODUCT(shape(d(1:ndim))))390 scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim)))) 397 391 DO i = 1, ntotal 398 392 array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, & 399 KIND= fftkind)393 KIND=wp) 400 394 END DO 401 395 … … 414 408 ! 415 409 !-- Formal parameters 416 COMPLEX( fftkind),DIMENSION(*), INTENT(INOUT) :: array417 INTEGER ,INTENT(IN) :: ntotal, npass, nspan418 LOGICAL, INTENT(IN) :: inv419 INTEGER, INTENT(OUT), OPTIONAL:: stat410 COMPLEX(wp), DIMENSION(*), INTENT(INOUT) :: array 411 INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan 412 INTEGER(iwp), INTENT(OUT), OPTIONAL:: stat 413 LOGICAL, INTENT(IN) :: inv 420 414 ! 421 415 !-- Local arrays 422 INTEGER, DIMENSION(BIT_SIZE(0)) :: factor423 COMPLEX(fftkind), DIMENSION(:), ALLOCATABLE :: ctmp424 REAL(fftkind), DIMENSION(:), ALLOCATABLE :: sine, cosine425 INTEGER, DIMENSION(:), ALLOCATABLE :: perm416 COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: ctmp 417 INTEGER(iwp), DIMENSION(BIT_SIZE(0)) :: factor 418 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: perm 419 REAL(wp), DIMENSION(:), ALLOCATABLE :: sine, cosine 426 420 ! 427 421 !-- Local scalars 428 INTEGER :: maxfactor, nfactor, nsquare, nperm422 INTEGER(iwp) :: maxfactor, nfactor, nsquare, nperm 429 423 ! 430 424 !-- Intrinsics used … … 476 470 ! 477 471 !-- Formal parameters 478 INTEGER , INTENT(IN) :: npass479 INTEGER , DIMENSION(*), INTENT(OUT):: factor480 INTEGER , INTENT(OUT):: nfactor, nsquare472 INTEGER(iwp), INTENT(IN) :: npass 473 INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor 474 INTEGER(iwp), INTENT(OUT):: nfactor, nsquare 481 475 ! 482 476 !-- Local scalars 483 INTEGER :: j, jj, k477 INTEGER(iwp):: j, jj, k 484 478 485 479 nfactor = 0 … … 541 535 ! 542 536 !-- Formal parameters 543 COMPLEX( fftkind),DIMENSION(*), INTENT(IN OUT):: array544 INTEGER, INTENT(IN) :: ntotal, npass, nspan545 INTEGER , DIMENSION(*), INTENT(IN) :: factor546 INTEGER , INTENT(IN) :: nfactor547 COMPLEX(fftkind), DIMENSION(*), INTENT(OUT) :: ctmp548 REAL(fftkind), DIMENSION(*), INTENT(OUT) :: sine, cosine549 LOGICAL, INTENT(IN) :: inv537 COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array 538 COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp 539 INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan 540 INTEGER(iwp), DIMENSION(*), INTENT(IN) :: factor 541 INTEGER(iwp), INTENT(IN) :: nfactor 542 LOGICAL, INTENT(IN) :: inv 543 REAL(wp), DIMENSION(*), INTENT(OUT) :: sine, cosine 550 544 ! 551 545 !-- Local scalars 552 INTEGER 553 INTEGER 554 INTEGER 555 INTEGER 556 REAL( fftkind):: s60, c72, s72, pi2, radf557 REAL( fftkind):: c1, s1, c2, s2, c3, s3, cd, sd, ak558 COMPLEX( fftkind):: cc, cj, ck, cjp, cjm, ckp, ckm546 INTEGER(iwp):: ii, ispan 547 INTEGER(iwp):: j, jc, jf, jj 548 INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4 549 INTEGER(iwp):: nn, nt 550 REAL(wp) :: s60, c72, s72, pi2, radf 551 REAL(wp) :: c1, s1, c2, s2, c3, s3, cd, sd, ak 552 COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm 559 553 560 554 c72 = cos72 … … 574 568 jc = nspan / npass 575 569 radf = pi2 * jc 576 pi2 = pi2 * 2.0_ fftkind!-- use 2 PI from here on570 pi2 = pi2 * 2.0_wp !-- use 2 PI from here on 577 571 578 572 ii = 0 … … 581 575 sd = radf / kspan 582 576 cd = SIN(sd) 583 cd = 2.0_ fftkind* cd * cd577 cd = 2.0_wp * cd * cd 584 578 sd = SIN(sd + sd) 585 579 kk = 1 … … 606 600 IF (kk > kspan) RETURN 607 601 DO 608 c1 = 1.0_ fftkind- cd602 c1 = 1.0_wp - cd 609 603 s1 = sd 610 604 DO … … 614 608 ck = array(kk) - array(k2) 615 609 array(kk) = array(kk) + array(k2) 616 array(k2) = ck * CMPLX(c1, s1, KIND= fftkind)610 array(k2) = ck * CMPLX(c1, s1, KIND=wp) 617 611 kk = k2 + kspan 618 612 IF (kk >= nt) EXIT … … 625 619 ak = c1 - (cd * c1 + sd * s1) 626 620 s1 = sd * c1 - cd * s1 + s1 627 c1 = 2.0_ fftkind- (ak * ak + s1 * s1)621 c1 = 2.0_wp - (ak * ak + s1 * s1) 628 622 s1 = s1 * c1 629 623 c1 = c1 * ak … … 641 635 642 636 DO 643 c1 = 1.0_ fftkind644 s1 = 0.0_ fftkind637 c1 = 1.0_wp 638 s1 = 0.0_wp 645 639 DO 646 640 DO … … 655 649 cjp = ckp - cjp 656 650 IF (inv) THEN 657 ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND= fftkind)658 ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND= fftkind)651 ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) 652 ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) 659 653 ELSE 660 ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND= fftkind)661 ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND= fftkind)654 ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp) 655 ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp) 662 656 END IF 663 657 ! 664 658 !-- Avoid useless multiplies 665 IF (s1 == 0.0_ fftkind) THEN659 IF (s1 == 0.0_wp) THEN 666 660 array(k1) = ckp 667 661 array(k2) = cjp 668 662 array(k3) = ckm 669 663 ELSE 670 array(k1) = ckp * CMPLX(c1, s1, KIND= fftkind)671 array(k2) = cjp * CMPLX(c2, s2, KIND= fftkind)672 array(k3) = ckm * CMPLX(c3, s3, KIND= fftkind)664 array(k1) = ckp * CMPLX(c1, s1, KIND=wp) 665 array(k2) = cjp * CMPLX(c2, s2, KIND=wp) 666 array(k3) = ckm * CMPLX(c3, s3, KIND=wp) 673 667 END IF 674 668 kk = k3 + kspan … … 678 672 c2 = c1 - (cd * c1 + sd * s1) 679 673 s1 = sd * c1 - cd * s1 + s1 680 c1 = 2.0_ fftkind- (c2 * c2 + s1 * s1)674 c1 = 2.0_wp - (c2 * c2 + s1 * s1) 681 675 s1 = s1 * c1 682 676 c1 = c1 * c2 … … 684 678 !-- Values of c2, c3, s2, s3 that will get used next time 685 679 c2 = c1 * c1 - s1 * s1 686 s2 = 2.0_ fftkind* c1 * s1680 s2 = 2.0_wp * c1 * s1 687 681 c3 = c2 * c1 - s2 * s1 688 682 s3 = c2 * s1 + s2 * c1 … … 712 706 array(kk) = ck + cj 713 707 ck = ck - CMPLX( & 714 0.5_ fftkind* REAL (cj), &715 0.5_ fftkind* AIMAG(cj), &716 KIND= fftkind)708 0.5_wp * REAL (cj), & 709 0.5_wp * AIMAG(cj), & 710 KIND=wp) 717 711 cj = CMPLX( & 718 712 (REAL (array(k1)) - REAL (array(k2))) * s60, & 719 713 (AIMAG(array(k1)) - AIMAG(array(k2))) * s60, & 720 KIND= fftkind)721 array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND= fftkind)722 array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND= fftkind)714 KIND=wp) 715 array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) 716 array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) 723 717 kk = k2 + kspan 724 718 IF (kk >= nn) EXIT … … 730 724 CASE (5) !-- transform for factor of 5 (optional code) 731 725 c2 = c72 * c72 - s72 * s72 732 s2 = 2.0_ fftkind* c72 * s72726 s2 = 2.0_wp * c72 * s72 733 727 DO 734 728 DO … … 744 738 array(kk) = cc + ckp + cjp 745 739 ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, & 746 KIND= fftkind) + &740 KIND=wp) + & 747 741 CMPLX(REAL(cjp) * c2, AIMAG(cjp) * c2, & 748 KIND= fftkind) + cc742 KIND=wp) + cc 749 743 cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, & 750 KIND= fftkind) + &744 KIND=wp) + & 751 745 CMPLX(REAL(cjm) * s2, AIMAG(cjm) * s2, & 752 KIND= fftkind)753 array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND= fftkind)754 array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND= fftkind)746 KIND=wp) 747 array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) 748 array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) 755 749 ck = CMPLX(REAL(ckp) * c2, AIMAG(ckp) * c2, & 756 KIND= fftkind) + &750 KIND=wp) + & 757 751 CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, & 758 KIND= fftkind) + cc752 KIND=wp) + cc 759 753 cj = CMPLX(REAL(ckm) * s2, AIMAG(ckm) * s2, & 760 KIND= fftkind) - &754 KIND=wp) - & 761 755 CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, & 762 KIND= fftkind)763 array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND= fftkind)764 array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND= fftkind)756 KIND=wp) 757 array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp) 758 array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp) 765 759 kk = k4 + kspan 766 760 IF (kk >= nn) EXIT … … 776 770 c1 = COS(s1) 777 771 s1 = SIN(s1) 778 cosine (jf) = 1.0_ fftkind779 sine (jf) = 0.0_ fftkind772 cosine (jf) = 1.0_wp 773 sine (jf) = 0.0_wp 780 774 j = 1 781 775 DO … … 816 810 jj = j 817 811 ck = cc 818 cj = (0.0_ fftkind, 0.0_fftkind)812 cj = (0.0_wp, 0.0_wp) 819 813 k = 1 820 814 DO … … 822 816 ck = ck + CMPLX( & 823 817 REAL (ctmp(k)) * cosine(jj), & 824 AIMAG(ctmp(k)) * cosine(jj), KIND= fftkind)818 AIMAG(ctmp(k)) * cosine(jj), KIND=wp) 825 819 k = k + 1 826 820 cj = cj + CMPLX( & 827 821 REAL (ctmp(k)) * sine(jj), & 828 AIMAG(ctmp(k)) * sine(jj), KIND= fftkind)822 AIMAG(ctmp(k)) * sine(jj), KIND=wp) 829 823 jj = jj + j 830 824 IF (jj > jf) jj = jj - jf … … 833 827 k = jf - j 834 828 array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), & 835 KIND= fftkind)829 KIND=wp) 836 830 array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), & 837 KIND= fftkind)831 KIND=wp) 838 832 j = j + 1 839 833 IF (j >= k) EXIT … … 852 846 kk = jc + 1 853 847 DO 854 c2 = 1.0_ fftkind- cd848 c2 = 1.0_wp - cd 855 849 s1 = sd 856 850 DO … … 860 854 DO 861 855 DO 862 array(kk) = CMPLX(c2, s2, KIND= fftkind) * array(kk)856 array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk) 863 857 kk = kk + ispan 864 858 IF (kk > nt) EXIT … … 872 866 c2 = c1 - (cd * c1 + sd * s1) 873 867 s1 = s1 + sd * c1 - cd * s1 874 c1 = 2.0_ fftkind- (c2 * c2 + s1 * s1)868 c1 = 2.0_wp - (c2 * c2 + s1 * s1) 875 869 s1 = s1 * c1 876 870 c2 = c2 * c1 … … 892 886 ! 893 887 !-- Formal parameters 894 COMPLEX( fftkind),DIMENSION(*), INTENT(IN OUT):: array895 INTEGER, INTENT(IN) :: ntotal, npass, nspan896 INTEGER , DIMENSION(*), INTENT(IN OUT):: factor897 INTEGER , INTENT(IN) :: nfactor, nsquare898 INTEGER , INTENT(IN) :: maxfactor899 COMPLEX(fftkind), DIMENSION(*), INTENT(OUT) :: ctmp900 INTEGER ,DIMENSION(*), INTENT(OUT) :: perm888 COMPLEX(wp), DIMENSION(*), INTENT(IN OUT):: array 889 COMPLEX(wp), DIMENSION(*), INTENT(OUT) :: ctmp 890 INTEGER(iwp), INTENT(IN) :: ntotal, npass, nspan 891 INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor 892 INTEGER(iwp), INTENT(IN) :: nfactor, nsquare 893 INTEGER(iwp), INTENT(IN) :: maxfactor 894 INTEGER(iwp), DIMENSION(*), INTENT(OUT) :: perm 901 895 ! 902 896 !-- Local scalars 903 INTEGER :: ii, ispan 904 INTEGER :: j, jc, jj 905 INTEGER :: k, kk, kspan, kt, k1, k2, k3 906 INTEGER :: nn, nt 907 COMPLEX(fftkind):: ck 908 897 COMPLEX(wp) :: ck 898 INTEGER(iwp):: ii, ispan 899 INTEGER(iwp):: j, jc, jj 900 INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3 901 INTEGER(iwp):: nn, nt 909 902 ! 910 903 !-- Permute the results to normal order---done in two stages
Note: See TracChangeset
for help on using the changeset viewer.
