annotate libcruft/fftpack/passf.f @ 7789:82be108cc558

First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents 30c606bec7a8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 subroutine passf (nac,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,wa)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2 dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 1 c1(ido,l1,ip) ,wa(1) ,c2(idl1,ip),
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 2 ch2(idl1,ip)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 idot = ido/2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 nt = ip*idl1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 ipp2 = ip+2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 ipph = (ip+1)/2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 idp = ip*ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 c
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 if (ido .lt. l1) go to 106
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 do 103 j=2,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 jc = ipp2-j
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 do 102 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 do 101 i=1,ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 101 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 102 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 103 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 do 105 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 do 104 i=1,ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 ch(i,k,1) = cc(i,1,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 104 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 105 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 go to 112
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 106 do 109 j=2,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 jc = ipp2-j
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 do 108 i=1,ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
30 do 107 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
31 ch(i,k,j) = cc(i,j,k)+cc(i,jc,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
32 ch(i,k,jc) = cc(i,j,k)-cc(i,jc,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
33 107 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
34 108 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
35 109 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
36 do 111 i=1,ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
37 do 110 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
38 ch(i,k,1) = cc(i,1,k)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
39 110 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
40 111 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
41 112 idl = 2-ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
42 inc = 0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
43 do 116 l=2,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
44 lc = ipp2-l
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
45 idl = idl+ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
46 do 113 ik=1,idl1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
47 c2(ik,l) = ch2(ik,1)+wa(idl-1)*ch2(ik,2)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
48 c2(ik,lc) = -wa(idl)*ch2(ik,ip)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
49 113 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
50 idlj = idl
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
51 inc = inc+ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
52 do 115 j=3,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53 jc = ipp2-j
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 idlj = idlj+inc
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
55 if (idlj .gt. idp) idlj = idlj-idp
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
56 war = wa(idlj-1)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
57 wai = wa(idlj)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
58 do 114 ik=1,idl1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
59 c2(ik,l) = c2(ik,l)+war*ch2(ik,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
60 c2(ik,lc) = c2(ik,lc)-wai*ch2(ik,jc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 114 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
62 115 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
63 116 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
64 do 118 j=2,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
65 do 117 ik=1,idl1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
66 ch2(ik,1) = ch2(ik,1)+ch2(ik,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
67 117 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
68 118 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
69 do 120 j=2,ipph
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
70 jc = ipp2-j
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
71 do 119 ik=2,idl1,2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
72 ch2(ik-1,j) = c2(ik-1,j)-c2(ik,jc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
73 ch2(ik-1,jc) = c2(ik-1,j)+c2(ik,jc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
74 ch2(ik,j) = c2(ik,j)+c2(ik-1,jc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
75 ch2(ik,jc) = c2(ik,j)-c2(ik-1,jc)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
76 119 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
77 120 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
78 nac = 1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
79 if (ido .eq. 2) return
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
80 nac = 0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
81 do 121 ik=1,idl1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
82 c2(ik,1) = ch2(ik,1)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
83 121 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
84 do 123 j=2,ip
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
85 do 122 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
86 c1(1,k,j) = ch(1,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
87 c1(2,k,j) = ch(2,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
88 122 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
89 123 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
90 if (idot .gt. l1) go to 127
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
91 idij = 0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
92 do 126 j=2,ip
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
93 idij = idij+2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
94 do 125 i=4,ido,2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
95 idij = idij+2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
96 do 124 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
97 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
98 c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
99 124 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
100 125 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
101 126 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
102 return
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
103 127 idj = 2-ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
104 do 130 j=2,ip
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
105 idj = idj+ido
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
106 do 129 k=1,l1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
107 idij = idj
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
108 do 128 i=4,ido,2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
109 idij = idij+2
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
110 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)+wa(idij)*ch(i,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
111 c1(i,k,j) = wa(idij-1)*ch(i,k,j)-wa(idij)*ch(i-1,k,j)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
112 128 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
113 129 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
114 130 continue
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
115 return
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
116 end