Initial commit

1 files changed, 179 insertions, 0 deletions

diff --git a/math/ieee/chap1/inishl.f b/math/ieee/chap1/inishl.f
new file mode 100644
index 00000000..925cd657
--- /dev/null
+++ b/math/ieee/chap1/inishl.f
@@ -0,0 +1,179 @@
+c
+c-----------------------------------------------------------------------
+c subroutine: inishl
+c   this subroutine initializes the wfta routine for a given
+c   value of the transform length n.  the factors of n are
+c   determined, the multiplication coefficients are calculated
+c   and stored in the array coef(.), the input and output
+c   permutation vectors are computed and stored in the arrays
+c   indx1(.) and indx2(.)
+c
+c-----------------------------------------------------------------------
+c
+      subroutine inishl(n,coef,xr,xi,indx1,indx2,ierr)
+      dimension coef(1),xr(1),xi(1)
+      integer s1,s2,s3,s4,indx1(1),indx2(1),p1
+      dimension co3(3),co4(4),co8(8),co9(11),co16(18),cda(18),cdb(11),
+     1cdc(9),cdd(6)
+      common na,nb,nc,nd,nd1,nd2,nd3,nd4
+c
+c   data statements assign short dft coefficients.
+c
+      data co4(1),co4(2),co4(3),co4(4)/4*1.0/
+c
+      data cda(1),cda(2),cda(3),cda(4),cda(5),cda(6),cda(7),
+     1 cda(8),cda(9),cda(10),cda(11),cda(12),cda(13),cda(14),
+     2 cda(15),cda(16),cda(17),cda(18)/18*1.0/
+c
+      data cdb(1),cdb(2),cdb(3),cdb(4),cdb(5),cdb(6),cdb(7),cdb(8),
+     1 cdb(9),cdb(10),cdb(11)/11*1.0/
+c
+      data ionce/1/
+c
+c   get multiplier constants
+c
+      if(ionce.ne.1) go to 20
+      call const(co3,co8,co16,co9,cdc,cdd)
+20    ionce=-1
+c
+c   following segment determines factors of n and chooses
+c   the appropriate short dft coefficients.
+c
+      iout=i1mach(2)
+      ierr=0
+      nd1=1
+      na=1
+      nb=1
+      nd2=1
+      nc=1
+      nd3=1
+      nd=1
+      nd4=1
+      if(n.le.0 .or. n.gt.5040) go to 190
+      if(16*(n/16).eq.n) go to 30
+      if(8*(n/8).eq.n) go to 40
+      if(4*(n/4).eq.n) go to 50
+      if(2*(n/2).ne.n) go to 70
+      nd1=2
+      na=2
+      cda(2)=1.0
+      go to 70
+30    nd1=18
+      na=16
+      do 31 j=1,18
+31    cda(j)=co16(j)
+      go to 70
+40    nd1=8
+      na=8
+      do 41 j=1,8
+41    cda(j)=co8(j)
+      go to 70
+50    nd1=4
+      na=4
+      do 51 j=1,4
+51    cda(j)=co4(j)
+70    if(3*(n/3).ne.n) go to 120
+      if(9*(n/9).eq.n) go to 100
+      nd2=3
+      nb=3
+      do 71 j=1,3
+71    cdb(j)=co3(j)
+      go to 120
+100   nd2=11
+      nb=9
+      do 110 j=1,11
+110   cdb(j)=co9(j)
+120   if(7*(n/7).ne.n) go to 160
+      nd3=9
+      nc=7
+160   if(5*(n/5).ne.n) go to 190
+      nd4=6
+      nd=5
+190   m=na*nb*nc*nd
+      if(m.eq.n) go to 250
+      write(iout,210)
+210   format(21h this n does not work)
+      ierr=-1
+      return
+c
+c   next segment generates the dft coefficients by
+c   multiplying together the short dft coefficients
+c
+250   j=1
+      do 300 n4=1,nd4
+      do 300 n3=1,nd3
+      do 300 n2=1,nd2
+      do 300 n1=1,nd1
+      coef(j)=cda(n1)*cdb(n2)*cdc(n3)*cdd(n4)
+      j=j+1
+300   continue
+c
+c   following segment forms the input indexing vector
+c
+      j=1
+      nu=nb*nc*nd
+      nv=na*nc*nd
+      nw=na*nb*nd
+      ny=na*nb*nc
+      k=1
+      do 440 n4=1,nd
+      do 430 n3=1,nc
+      do 420 n2=1,nb
+      do 410 n1=1,na
+405   if(k.le.n) go to 408
+      k=k-n
+      go to 405
+408   indx1(j)=k
+      j=j+1
+410   k=k+nu
+420   k=k+nv
+430   k=k+nw
+440   k=k+ny
+c
+c   following segment forms the output indexing vector
+c
+      m=1
+      s1=0
+      s2=0
+      s3=0
+      s4=0
+      if(na.eq.1) go to 530
+520   p1=m*nu-1
+      if((p1/na)*na.eq.p1) go to 510
+      m=m+1
+      go to 520
+510   s1=p1+1
+530   if(nb.eq.1) go to 540
+      m=1
+550   p1=m*nv-1
+      if((p1/nb)*nb.eq.p1) go to 560
+      m=m+1
+      go to 550
+560   s2=p1+1
+540   if(nc.eq.1) go to 630
+      m=1
+620   p1=m*nw-1
+      if((p1/nc)*nc.eq.p1) go to 610
+      m=m+1
+      go to 620
+610   s3=p1+1
+630   if(nd.eq.1) go to 660
+      m=1
+640   p1=m*ny-1
+      if((p1/nd)*nd.eq.p1) go to 650
+      m=m+1
+      go to 640
+650   s4=p1+1
+660   j=1
+      do 810 n4=1,nd
+      do 810 n3=1,nc
+      do 810 n2=1,nb
+      do 810 n1=1,na
+      indx2(j)=s1*(n1-1)+s2*(n2-1)+s3*(n3-1)+s4*(n4-1)+1
+900   if(indx2(j).le.n) go to 910
+      indx2(j)=indx2(j)-n
+      go to 900
+910   j=j+1
+810   continue
+      return
+      end