| [2759] | 1 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
|---|
| 2 | ! |
|---|
| 3 | ! FFTPACK 5.0 |
|---|
| 4 | ! Copyright (C) 1995-2004, Scientific Computing Division, |
|---|
| 5 | ! University Corporation for Atmospheric Research |
|---|
| 6 | ! Licensed under the GNU General Public License (GPL) |
|---|
| 7 | ! |
|---|
| 8 | ! Authors: Paul N. Swarztrauber and Richard A. Valent |
|---|
| 9 | ! |
|---|
| 10 | ! $Id: mrftb1.f,v 1.2 2004/06/15 21:29:20 rodney Exp $ |
|---|
| 11 | ! |
|---|
| 12 | !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC |
|---|
| 13 | |
|---|
| 14 | SUBROUTINE MRFTB1 (M,IM,N,IN,C,CH,WA,FAC) |
|---|
| 15 | REAL CH(M,*), C(IN,*), WA(N) ,FAC(15) |
|---|
| 16 | ! |
|---|
| 17 | NF = FAC(2) |
|---|
| 18 | NA = 0 |
|---|
| 19 | DO 10 K1=1,NF |
|---|
| 20 | IP = FAC(K1+2) |
|---|
| 21 | NA = 1-NA |
|---|
| 22 | IF(IP .LE. 5) GO TO 10 |
|---|
| 23 | IF(K1 .EQ. NF) GO TO 10 |
|---|
| 24 | NA = 1-NA |
|---|
| 25 | 10 END DO |
|---|
| 26 | HALF = .5 |
|---|
| 27 | HALFM = -.5 |
|---|
| 28 | MODN = MOD(N,2) |
|---|
| 29 | NL = N-2 |
|---|
| 30 | IF(MODN .NE. 0) NL = N-1 |
|---|
| 31 | IF (NA .EQ. 0) GO TO 120 |
|---|
| 32 | M2 = 1-IM |
|---|
| 33 | DO 117 I=1,M |
|---|
| 34 | M2 = M2+IM |
|---|
| 35 | CH(I,1) = C(M2,1) |
|---|
| 36 | CH(I,N) = C(M2,N) |
|---|
| 37 | 117 END DO |
|---|
| 38 | DO 118 J=2,NL,2 |
|---|
| 39 | M2 = 1-IM |
|---|
| 40 | DO 118 I=1,M |
|---|
| 41 | M2 = M2+IM |
|---|
| 42 | CH(I,J) = HALF*C(M2,J) |
|---|
| 43 | CH(I,J+1) = HALFM*C(M2,J+1) |
|---|
| 44 | 118 CONTINUE |
|---|
| 45 | GO TO 124 |
|---|
| 46 | 120 continue |
|---|
| 47 | DO 122 J=2,NL,2 |
|---|
| 48 | M2 = 1-IM |
|---|
| 49 | DO 122 I=1,M |
|---|
| 50 | M2 = M2+IM |
|---|
| 51 | C(M2,J) = HALF*C(M2,J) |
|---|
| 52 | C(M2,J+1) = HALFM*C(M2,J+1) |
|---|
| 53 | 122 CONTINUE |
|---|
| 54 | 124 L1 = 1 |
|---|
| 55 | IW = 1 |
|---|
| 56 | DO 116 K1=1,NF |
|---|
| 57 | IP = FAC(K1+2) |
|---|
| 58 | L2 = IP*L1 |
|---|
| 59 | IDO = N/L2 |
|---|
| 60 | IDL1 = IDO*L1 |
|---|
| 61 | IF (IP .NE. 4) GO TO 103 |
|---|
| 62 | IX2 = IW+IDO |
|---|
| 63 | IX3 = IX2+IDO |
|---|
| 64 | IF (NA .NE. 0) GO TO 101 |
|---|
| 65 | CALL MRADB4 (M,IDO,L1,C,IM,IN,CH,1,M,WA(IW),WA(IX2),WA(IX3)) |
|---|
| 66 | GO TO 102 |
|---|
| 67 | 101 CALL MRADB4 (M,IDO,L1,CH,1,M,C,IM,IN,WA(IW),WA(IX2),WA(IX3)) |
|---|
| 68 | 102 NA = 1-NA |
|---|
| 69 | GO TO 115 |
|---|
| 70 | 103 IF (IP .NE. 2) GO TO 106 |
|---|
| 71 | IF (NA .NE. 0) GO TO 104 |
|---|
| 72 | CALL MRADB2 (M,IDO,L1,C,IM,IN,CH,1,M,WA(IW)) |
|---|
| 73 | GO TO 105 |
|---|
| 74 | 104 CALL MRADB2 (M,IDO,L1,CH,1,M,C,IM,IN,WA(IW)) |
|---|
| 75 | 105 NA = 1-NA |
|---|
| 76 | GO TO 115 |
|---|
| 77 | 106 IF (IP .NE. 3) GO TO 109 |
|---|
| 78 | IX2 = IW+IDO |
|---|
| 79 | IF (NA .NE. 0) GO TO 107 |
|---|
| 80 | CALL MRADB3 (M,IDO,L1,C,IM,IN,CH,1,M,WA(IW),WA(IX2)) |
|---|
| 81 | GO TO 108 |
|---|
| 82 | 107 CALL MRADB3 (M,IDO,L1,CH,1,M,C,IM,IN,WA(IW),WA(IX2)) |
|---|
| 83 | 108 NA = 1-NA |
|---|
| 84 | GO TO 115 |
|---|
| 85 | 109 IF (IP .NE. 5) GO TO 112 |
|---|
| 86 | IX2 = IW+IDO |
|---|
| 87 | IX3 = IX2+IDO |
|---|
| 88 | IX4 = IX3+IDO |
|---|
| 89 | IF (NA .NE. 0) GO TO 110 |
|---|
| 90 | CALL MRADB5 (M,IDO,L1,C,IM,IN,CH,1,M,WA(IW),WA(IX2), & |
|---|
| 91 | & WA(IX3),WA(IX4)) |
|---|
| 92 | GO TO 111 |
|---|
| 93 | 110 CALL MRADB5 (M,IDO,L1,CH,1,M,C,IM,IN,WA(IW),WA(IX2), & |
|---|
| 94 | & WA(IX3),WA(IX4)) |
|---|
| 95 | 111 NA = 1-NA |
|---|
| 96 | GO TO 115 |
|---|
| 97 | 112 IF (NA .NE. 0) GO TO 113 |
|---|
| 98 | CALL MRADBG (M,IDO,IP,L1,IDL1,C,C,C,IM,IN,CH,CH,1, & |
|---|
| 99 | & M,WA(IW)) |
|---|
| 100 | GO TO 114 |
|---|
| 101 | 113 CALL MRADBG (M,IDO,IP,L1,IDL1,CH,CH,CH,1,M,C,C,IM, & |
|---|
| 102 | & IN,WA(IW)) |
|---|
| 103 | 114 IF (IDO .EQ. 1) NA = 1-NA |
|---|
| 104 | 115 L1 = L2 |
|---|
| 105 | IW = IW+(IP-1)*IDO |
|---|
| 106 | 116 END DO |
|---|
| 107 | RETURN |
|---|
| 108 | END |
|---|
