1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
c
c-----------------------------------------------------------------------
c subroutine: wfta
c winograd fourier transform algorithm
c-----------------------------------------------------------------------
c
subroutine wfta(xr,xi,n,invrs,init,ierr)
dimension xr(1),xi(1)
c
c inputs:
c n-- transform length. must be formed as the product of
c relatively prime integers from the set:
c 2,3,4,5,7,8,9,16
c thus the largest possible value of n is 5040.
c xr(.)-- array that holds the real part of the data
c to be transformed.
c xi(.)-- array that holds the imaginary part of the
c data to be transformed.
c invrs-- parameter that flags whether or not the inverse
c transform is to be calculated. a division by n
c is included in the inverse.
c invrs = 1 yields inverse transform
c invrs .ne. 1 gives forward transform
c init-- parameter that flags whether or not the program
c is to be initialized for this value of n. the
c initialization is performed only once in order to
c to speed up the computation on succeeding calls
c to the wfta routine, when n is held fixed.
c init = 0 results in initialization.
c ierr-- error code that is negative when the wfta
c terminates incorrectly.
c 0 = successful completion
c -1 = this value of n does not factor properly
c -2 = an initialization has not been done for
c this value of n.
c
c
c the following two cards may be changed if the maximum
c desired transform length is less than 5040
c
c *********************************************************************
dimension sr(1782),si(1782),coef(1782)
integer indx1(1008),indx2(1008)
c *********************************************************************
c
common na,nb,nc,nd,nd1,nd2,nd3,nd4
c
c test for initial run
c
if(init.eq.0) call inishl(n,coef,xr,xi,indx1,indx2,ierr)
c
if(ierr.lt.0) return
m=na*nb*nc*nd
if(m.eq.n) go to 100
ierr=-2
return
c
c error(-2)-- program not initialized for this value of n
c
100 nmult=nd1*nd2*nd3*nd4
c
c the following code maps the data arrays xr and xi to
c the working arrays sr and si via the mapping vector
c indx1(.). the permutation of the data follows the
c sino correspondence of the chinese remainder theorem.
c
j=1
k=1
inc1=nd1-na
inc2=nd1*(nd2-nb)
inc3=nd1*nd2*(nd3-nc)
do 140 n4=1,nd
do 130 n3=1,nc
do 120 n2=1,nb
do 110 n1=1,na
ind=indx1(k)
sr(j)=xr(ind)
si(j)=xi(ind)
k=k+1
110 j=j+1
120 j=j+inc1
130 j=j+inc2
140 j=j+inc3
c
c do the pre-weave modules
c
call weave1(sr,si)
c
c the following loop performs all the multiplications of the
c winograd fourier transform algorithm. the multiplication
c coefficients are stored on the initialization pass in the
c array coef(.).
c
do 200 j=1,nmult
sr(j)=sr(j)*coef(j)
si(j)=si(j)*coef(j)
200 continue
c
c do the post-weave modules
c
call weave2(sr,si)
c
c
c the following code maps the working arrays sr and si
c to the data arrays xr and xi via the mapping vector
c indx2(.). the permutation of the data follows the
c chinese remainder theorem.
c
j=1
k=1
inc1=nd1-na
inc2=nd1*(nd2-nb)
inc3=nd1*nd2*(nd3-nc)
c
c check for inverse
c
if(invrs.eq.1) go to 400
do 340 n4=1,nd
do 330 n3=1,nc
do 320 n2=1,nb
do 310 n1=1,na
kndx=indx2(k)
xr(kndx)=sr(j)
xi(kndx)=si(j)
k=k+1
310 j=j+1
320 j=j+inc1
330 j=j+inc2
340 j=j+inc3
return
c
c different permutation for the inverse
c
400 fn=float(n)
np2=n+2
indx2(1)=n+1
do 440 n4=1,nd
do 430 n3=1,nc
do 420 n2=1,nb
do 410 n1=1,na
kndx=np2-indx2(k)
xr(kndx)=sr(j)/fn
xi(kndx)=si(j)/fn
k=k+1
410 j=j+1
420 j=j+inc1
430 j=j+inc2
440 j=j+inc3
return
end
|
