c $Id: convsetsf,v 1.1.2.1 2003/03/13 18:56:54 muchomas Exp $ c $Log: convsetsf,v $ c Revision 1.1.2.1 2003/03/13 18:56:54 muchomas c Move convsets.f to fortran-source text file convsetsf for use with f2c. c Put f2c result in convsets.c, so as to eliminate need for fortran compilers. c c Revision 1.1.2.2 2002/10/28 22:07:29 matthew c Non-orthogonal Laplaciancvs add coeff_d.c! Works for orthorhombic cells, still in need of debugging for other unit cells. c c Revision 1.1.1.1 2001/04/18 17:15:14 daykov c starting project c c Revision 1.1 1999/09/14 19:40:30 torkel c Initial revision c c Revision 1.12 1999/07/04 20:48:04 torkel c Full O implemented. c c Revision 1.11 1999/07/02 07:04:05 torkel c the convolution in Osame seems to work for this version. c c Revision 1.10 1999/07/02 03:11:48 torkel c In the middle of making Osame work c c Revision 1.9 1999/06/23 22:20:05 torkel c revised version where Obelow is dagger of Oabove for grid without c >> siblings. c c Revision 1.8 1999/06/21 18:56:14 torkel c updated version for multilevel operators c c Revision 1.7 1999/06/15 23:24:45 torkel c debugging, trying to make Obelow work c c Revision 1.6 1999/06/11 12:53:50 torkel c including convolutions for obelow c c Revision 1.4 1999/04/25 22:32:56 muchomas c *** empty log message *** c c Revision 1.3 1999/04/25 18:35:39 muchomas c Contains new routines for same-level operator convolutions c Working toward more arbitrary length convolution routines c c Revision 1.2 1999/04/20 15:59:37 muchomas c Fully documented code c c Revision 1.1 1999/04/19 18:12:11 muchomas c Initial revision c c Revision 1.6 1999/04/17 22:00:09 muchomas c Working dagger operators, still slow... c c Revision 1.5 1999/04/14 19:48:37 muchomas c Working transpose transform, but doesn't take input/output arrays -- yet! c c Revision 1.4 1999/04/13 15:11:36 muchomas c Does ghosts, making compact version, but still in larger array c c Revision 1.3 1999/04/13 14:32:53 muchomas c Idag that computes only only ghosts, leaving sparse data distribution c c Revision 1.2 1999/04/11 19:53:46 muchomas c Working forward/inverse xform package c c Revision 1.1 1999/04/10 23:04:29 muchomas c Initial revision c c Revision 1.3 1999/04/10 22:40:41 muchomas c Now using log functionality c c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 3 (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsets3(out,outm3,outm2,outm1,outp1,outp2,outp3, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3 integer i,j,out2 c First, make sure pad is zero do i=n,npad-1 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) do i=npad-1,3,-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) enddo enddo out2=2*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j) enddo out2=0 do j=1+out2,area+out2 out(j)= & t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=0,n-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsets5(out,outm5,outm4,outm3,outm2,outm1, & outp1,outp2,outp3,outp4,outp5, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm5(*),outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*),outp5(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,out2 c First, make sure pad is zero do i=n,npad-1 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) do i=npad-1,5,-1 out2=i*stride do j=1+out2,area+out2 out(j)= $ t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) $ +t4*outm4(j)+t5*outm5(j) enddo enddo out2=4*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j)+t4*outm4(j) enddo out2=3*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) enddo out2=2*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=1*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j) enddo out2=0*stride do j=1+out2,area+out2 out(j)= & t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) $ +t4*outp4(j)+t5*outp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsper5(out,outm5,outm4,outm3,outm2,outm1, & outp1,outp2,outp3,outp4,outp5, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm5(*),outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*),outp5(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,out2 c First, make sure pad is zero do i=n,npad-1 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) i=npad-1 out2=i*stride do j=1+out2,area+out2 out(j)=t0*out(j-(i-4)*stride)+t1*out(j-(i-3)*stride) $ +t2*out(j-(i-2)*stride)+t3*out(j-(i-1)*stride) $ +t4*out(j-i*stride)+t5*outm5(j) enddo i=npad-2 out2=i*stride do j=1+out2,area+out2 out(j)=t0*out(j-(i-3)*stride)+t1*out(j-(i-2)*stride) $ +t2*out(j-(i-1)*stride)+t3*out(j-i*stride) $ +t4*outm4(j)+t5*outm5(j) enddo i=npad-3 out2=i*stride do j=1+out2,area+out2 out(j)=t0*out(j-(i-2)*stride)+t1*out(j-(i-1)*stride) $ +t2*out(j-i*stride)+t3*outm3(j) $ +t4*outm4(j)+t5*outm5(j) enddo i=npad-4 out2=i*stride do j=1+out2,area+out2 out(j)=t0*out(j-(i-1)*stride)+t1*out(j-i*stride) $ +t2*outm2(j)+t3*outm3(j) $ +t4*outm4(j)+t5*outm5(j) enddo i=npad-5 out2=i*stride do j=1+out2,area+out2 out(j)=t0*out(j-i*stride)+t1*outm1(j) $ +t2*outm2(j)+t3*outm3(j) $ +t4*outm4(j)+t5*outm5(j) enddo do i=npad-6,5,-1 out2=i*stride do j=1+out2,area+out2 out(j)= $ t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) $ +t4*outm4(j)+t5*outm5(j) enddo enddo out2=4*stride do j=1+out2,area+out2 out(j)=out(j+(npad-5)*stride) c & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j)+t4*outm4(j) enddo out2=3*stride do j=1+out2,area+out2 out(j)=out(j+(npad-5)*stride) c & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) enddo out2=2*stride do j=1+out2,area+out2 out(j)=out(j+(npad-5)*stride) c & t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=1*stride do j=1+out2,area+out2 out(j)=out(j+(npad-5)*stride) c & t0*out(j)+t1*outm1(j) enddo out2=0*stride do j=1+out2,area+out2 out(j)=out(j+(npad-5)*stride) c & t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) $ +t4*outp4(j)+t5*outp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsio5(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,tr,tf) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) do i=npad-1,5,-1 work2=i*stride do j=1+work2,area+work2 work(j)= $ t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo enddo work2=4*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j) & +t3*workm3(j)+t4*workm4(j) enddo work2=3*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) enddo work2=2*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j) enddo work2=1*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j) enddo work2=0*stride do j=1+work2,area+work2 work(j)= & t0*work(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)= & t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsio5full(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,filter) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision filter(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=filter(1) t1=filter(2) t2=filter(3) t3=filter(4) t4=filter(5) t5=filter(6) if (t0 /= 0.0) then c symmetric filter i=0 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo do i=5,n-6 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo i=n-5 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j) enddo i=n-4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) enddo i=n-3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j) enddo i=n-2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j) enddo i=n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j)+t0*work(j) enddo c end of symmetric part else c antisymmetric filter i=0 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo do i=5,n-6 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo i=n-5 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j) enddo i=n-4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) enddo i=n-3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j) enddo i=n-2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j) enddo i=n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=-1.*t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j)+t0*work(j) enddo c end of antisymmetric part endif return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsioppercum5(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,tr,tf) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) i=npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-4)*stride)+t1*work(j-(i-3)*stride) $ +t2*work(j-(i-2)*stride)+t3*work(j-(i-1)*stride) $ +t4*work(j-i*stride)+t5*workm5(j) enddo i=npad-2 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-3)*stride)+t1*work(j-(i-2)*stride) $ +t2*work(j-(i-1)*stride)+t3*work(j-i*stride) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-3 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-2)*stride)+t1*work(j-(i-1)*stride) $ +t2*work(j-i*stride)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-4 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-1)*stride)+t1*work(j-i*stride) $ +t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-5 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-i*stride)+t1*workm1(j) $ +t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo do i=npad-6,5,-1 work2=i*stride do j=1+work2,area+work2 work(j)= $ t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo enddo work2=4*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j)+t4*workm4(j) enddo work2=3*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) enddo work2=2*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j) enddo work2=1*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j) enddo work2=0*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)= & t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsiop5(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,tr,tf) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) do i=npad-1,5,-1 work2=i*stride do j=1+work2,area+work2 work(j)= $ t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo enddo work2=4*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j) & +t3*workm3(j)+t4*workm4(j) enddo work2=3*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) enddo work2=2*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j)+t2*workm2(j) enddo work2=1*stride do j=1+work2,area+work2 work(j)= & t0*work(j)+t1*workm1(j) enddo work2=0*stride do j=1+work2,area+work2 work(j)= & t0*work(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+ & t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsiop5full(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,filter) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision filter(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=filter(1) t1=filter(2) t2=filter(3) t3=filter(4) t4=filter(5) t5=filter(6) if (t0 /= 0.0) then c symmetric part i=0 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t0*work(j)+t1*workp1(j)+t2*workp2(j) & +t3*workp3(j)+t4*workp4(j)+t5*workp5(j) enddo i=1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo do i=5,n-6 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo i=n-5 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j) enddo i=n-4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) enddo i=n-3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j) enddo i=n-2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j) & +t0*work(j)+t1*workp1(j) enddo i=n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t5*workm5(j)+t4*workm4(j)+t3*workm3(j) & +t2*workm2(j)+t1*workm1(j)+t0*work(j) enddo c end of symmetric part else c antisymmetric part i=0 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+t0*work(j)+t1*workp1(j)+t2*workp2(j) & +t3*workp3(j)+t4*workp4(j)+t5*workp5(j) enddo i=1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo i=4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo do i=5,n-6 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo i=n-5 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j) enddo i=n-4 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) enddo i=n-3 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j)+t2*workp2(j) enddo i=n-2 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j) & +t0*work(j)+t1*workp1(j) enddo i=n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)-t5*workm5(j)-t4*workm4(j)-t3*workm3(j) & -t2*workm2(j)-t1*workm1(j)+t0*work(j) enddo c end of antisymmetric part end if return end c Back and forth, linear convolution of sets of points using padding c for filters of extent +/- 5 (IN PLACE) c c out: input/output array for in place computation c outm5-outp5,out: pointers in in/output array to c start of data sets c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsiopper5(datout, & work,workm5,workm4,workm3,workm2,workm1, & workp1,workp2,workp3,workp4,workp5, & area,stride,n,npad, & stridel,tr,tf) double precision datout(*) double precision work(*) double precision workm5(*),workm4(*),workm3(*),workm2(*),workm1(*) double precision workp1(*),workp2(*),workp3(*),workp4(*),workp5(*) integer stride,area,n,npad,stridel double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5 integer i,j,work2,out2s,jl c First, make sure pad is zero do i=n,npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) i=npad-1 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-4)*stride)+t1*work(j-(i-3)*stride) $ +t2*work(j-(i-2)*stride)+t3*work(j-(i-1)*stride) $ +t4*work(j-i*stride)+t5*workm5(j) enddo i=npad-2 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-3)*stride)+t1*work(j-(i-2)*stride) $ +t2*work(j-(i-1)*stride)+t3*work(j-i*stride) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-3 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-2)*stride)+t1*work(j-(i-1)*stride) $ +t2*work(j-i*stride)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-4 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-(i-1)*stride)+t1*work(j-i*stride) $ +t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo i=npad-5 work2=i*stride do j=1+work2,area+work2 work(j)=t0*work(j-i*stride)+t1*workm1(j) $ +t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo do i=npad-6,5,-1 work2=i*stride do j=1+work2,area+work2 work(j)= $ t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) $ +t4*workm4(j)+t5*workm5(j) enddo enddo work2=4*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j)+t4*workm4(j) enddo work2=3*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j)+t3*workm3(j) enddo work2=2*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j)+t2*workm2(j) enddo work2=1*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j)+t1*workm1(j) enddo work2=0*stride do j=1+work2,area+work2 work(j)=work(j+(npad-5)*stride) c & t0*work(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) do i=0,n-1 work2=i*stride out2s=i*(stridel-stride) do j=1+work2,area+work2 jl=j+out2s datout(jl)=datout(jl)+ & t0*work(j)+t1*workp1(j)+t2*workp2(j)+t3*workp3(j) & +t4*workp4(j)+t5*workp5(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding, c computing data ONLY on ghostpoints (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine pconvsets3(out,outm3,outm2,outm1,outp1,outp2,outp3, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3 integer i,j,out2 c First, make sure pad is zero for forward phase. c For pretty code the last number would be npad-1, but technically c npad-2 (for extra speed) is all that is needed. do i=n,npad-2 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) do i=n-1,3,-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) enddo enddo out2=2*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j) enddo out2=0 do j=1+out2,area+out2 out(j)= & t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=0,n-1,2 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding, c computing data ONLY on ghostpoints (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine pconvsetsio3full(out,in,inm3,inm2,inm1,inp1,inp2,inp3, & area,stride,n,npad,filt) implicit none double precision out(*) double precision in(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) integer stride,area,n,npad double precision filt(*) double precision t0,t1,t3 integer i,j,out2 c Forward phase (padding now ensures room!) t0=filt(1) t1=filt(2) t3=filt(4) i=0 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inp3(j) enddo i=2 out2=i*stride do j=1+out2,area+out2 out(j)=in(j)+t1*inp1(j)+t3*inp3(j) enddo do i=4,n-6,2 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inm3(j)+t1*inm1(j)+ & in(j)+t1*inp1(j)+t3*inp3(j) enddo enddo i=n-4 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inm3(j)+t1*inm1(j)+in(j) enddo i=n-2 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inm3(j) enddo end c Back and forth, linear convolution of sets of points using padding, c computing data ONLY on ghostpoints (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine pconvsetsio8full(out,in,inm8,inm7,inm6,inm5,inm4, & inm3,inm2,inm1,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, & area,stride,n,npad,filt) implicit none double precision out(*) double precision in(*) double precision inm8(*),inm7(*),inm6(*),inm5(*),inm4(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) double precision inp4(*),inp5(*),inp6(*),inp7(*),inp8(*) integer stride,area,n,npad double precision filt(*) double precision t0,t2,t1,t3,t4,t5,t6,t7,t8 integer i,j,out2 c Forward phase (padding now ensures room!) t0=filt(1) t1=filt(2) t2=filt(3) t3=filt(4) t4=filt(5) t5=filt(6) t6=filt(7) t7=filt(8) t8=filt(9) if (t0 /= 0.0) then c symmetric part i=0 out2=i*stride do j=1+out2,area+out2 out(j)=t8*inp8(j) enddo i=2 out2=i*stride do j=1+out2,area+out2 out(j)=t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=4 out2=i*stride do j=1+out2,area+out2 out(j)=t4*inp4(j)+t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=6 out2=i*stride do j=1+out2,area+out2 out(j)=t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=8 out2=i*stride do j=1+out2,area+out2 out(j)=t0*in(j)+t1*inp1(j)+t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=10 out2=i*stride do j=1+out2,area+out2 out(j)=t2*inm2(j)+t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=12 out2=i*stride do j=1+out2,area+out2 out(j)=t4*inm4(j)+t3*inm3(j)+t2*inm2(j)+t1*inm1(j) & +t0*in(j)+t1*inp1(j)+t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=14 out2=i*stride do j=1+out2,area+out2 out(j)=t6*inm6(j)+t5*inm5(j)+t4*inm4(j)+t3*inm3(j) & +t2*inm2(j)+t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo do i=16,n-18,2 out2=i*stride do j=1+out2,area+out2 out(j)=t8*inm8(j)+t7*inm7(j)+t6*inm6(j)+t5*inm5(j) & +t4*inm4(j)+t3*inm3(j) & +t2*inm2(j)+t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo enddo i=n-16 out2=i*stride do j=1+out2,area+out2 out(j)=t7*inp7(j)+t6*inp6(j)+t5*inp5(j)+t4*inp4(j) & +t3*inp3(j)+t2*inp2(j)+t1*inp1(j)+t0*in(j) & +t1*inm1(j)+t2*inm2(j)+t3*inm3(j)+t4*inm4(j) & +t5*inm5(j)+t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-14 out2=i*stride do j=1+out2,area+out2 out(j)=t5*inp5(j)+t4*inp4(j)+t3*inp3(j)+t2*inp2(j) & +t1*inp1(j)+t0*in(j)+t1*inm1(j)+t2*inm2(j) & +t3*inm3(j)+t4*inm4(j)+t5*inm5(j) & +t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-12 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inp3(j)+t2*inp2(j)+t1*inp1(j)+t0*in(j) & +t1*inm1(j)+t2*inm2(j)+t3*inm3(j)+t4*inm4(j)+t5*inm5(j) & +t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-10 out2=i*stride do j=1+out2,area+out2 out(j)=t1*inp1(j)+t0*in(j)+t1*inm1(j)+t2*inm2(j) & +t3*inm3(j)+t4*inm4(j)+t5*inm5(j) & +t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-8 out2=i*stride do j=1+out2,area+out2 out(j)=t1*inm1(j)+t2*inm2(j)+t3*inm3(j)+t4*inm4(j)+t5*inm5(j) & +t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-6 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inm3(j)+t4*inm4(j)+t5*inm5(j) & +t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-4 out2=i*stride do j=1+out2,area+out2 out(j)=t5*inm5(j)+t6*inm6(j)+t7*inm7(j)+t8*inm8(j) enddo i=n-2 out2=i*stride do j=1+out2,area+out2 out(j)=t7*inm7(j)+t8*inm8(j) enddo c end of symmetric part else c antisymmetric part i=0 out2=i*stride do j=1+out2,area+out2 out(j)=t8*inp8(j) enddo i=2 out2=i*stride do j=1+out2,area+out2 out(j)=t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=4 out2=i*stride do j=1+out2,area+out2 out(j)=t4*inp4(j)+t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=6 out2=i*stride do j=1+out2,area+out2 out(j)=t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=8 out2=i*stride do j=1+out2,area+out2 out(j)=t0*in(j)+t1*inp1(j)+t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=10 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t2*inm2(j)-t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=12 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t4*inm4(j)-t3*inm3(j)-t2*inm2(j)-t1*inm1(j) & +t0*in(j)+t1*inp1(j)+t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo i=14 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t6*inm6(j)-t5*inm5(j)-t4*inm4(j)-t3*inm3(j) & -t2*inm2(j)-t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo do i=16,n-18,2 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t8*inm8(j)-t7*inm7(j)-t6*inm6(j)-t5*inm5(j) & -t4*inm4(j)-t3*inm3(j) & -t2*inm2(j)-t1*inm1(j)+t0*in(j)+t1*inp1(j) & +t2*inp2(j)+t3*inp3(j)+t4*inp4(j) & +t5*inp5(j)+t6*inp6(j)+t7*inp7(j)+t8*inp8(j) enddo enddo i=n-16 out2=i*stride do j=1+out2,area+out2 out(j)=t7*inp7(j)+t6*inp6(j)+t5*inp5(j)+t4*inp4(j) & +t3*inp3(j)+t2*inp2(j)+t1*inp1(j)+t0*in(j) & -t1*inm1(j)-t2*inm2(j)-t3*inm3(j)-t4*inm4(j) & -t5*inm5(j)-t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-14 out2=i*stride do j=1+out2,area+out2 out(j)=t5*inp5(j)+t4*inp4(j)+t3*inp3(j)+t2*inp2(j) & +t1*inp1(j)+t0*in(j)-t1*inm1(j)-t2*inm2(j) & -t3*inm3(j)-t4*inm4(j)-t5*inm5(j) & -t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-12 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inp3(j)+t2*inp2(j)+t1*inp1(j)+t0*in(j) & -t1*inm1(j)-t2*inm2(j)-t3*inm3(j)-t4*inm4(j)-t5*inm5(j) & -t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-10 out2=i*stride do j=1+out2,area+out2 out(j)=t1*inp1(j)+t0*in(j)-t1*inm1(j)-t2*inm2(j) & -t3*inm3(j)-t4*inm4(j)-t5*inm5(j) & -t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-8 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t1*inm1(j)-t2*inm2(j)-t3*inm3(j)-t4*inm4(j) & -t5*inm5(j)-t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-6 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t3*inm3(j)-t4*inm4(j)-t5*inm5(j) & -t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-4 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t5*inm5(j)-t6*inm6(j)-t7*inm7(j)-t8*inm8(j) enddo i=n-2 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t7*inm7(j)-t8*inm8(j) enddo c end of antisymmetric part end if end c Back and forth, linear convolution of sets of points using padding, c computing data ONLY on ghostpoints (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine pconvsets3new(outx,out,outm3,outm2,outm1,outp1, & outp2,outp3, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outx(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3 integer i,j,out2 c First, make sure pad is zero for forward phase. c For pretty code the last number would be npad-1, but technically c npad-2 (for extra speed) is all that is needed. do i=n,npad-2 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) do i=n-1,3,-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) enddo enddo out2=2*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j) enddo out2=0 do j=1+out2,area+out2 out(j)= & t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=0,n-1,2 out2=i*stride do j=1+out2,area+out2 outx(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding, c computing data ONLY on ghostpoints (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine pconvsets8( & out,outm8,outm7,outm6,outm5,outm4, & outm3,outm2,outm1, & outp1,outp2,outp3, & outp4,outp5,outp6,outp7,outp8, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm8(*),outm7(*) double precision outm6(*),outm5(*),outm4(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) double precision outp4(*),outp5(*),outp6(*) double precision outp7(*),outp8(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2 c First, make sure pad is zero for forward phase. c For pretty code the last number would be npad-1, but technically c npad-2 (for extra speed) is all that is needed. do i=n,npad-2 out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo enddo c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) t6=tr(7) t7=tr(8) t8=tr(9) do i=n-1,8,-1 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j)+t3*outm3(j) & +t4*outm4(j)+t5*outm5(j)+t6*outm6(j)+t7*outm7(j) & +t8*outm8(j) enddo enddo out2=7*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) & +t3*outm3(j)+t4*outm4(j)+t5*outm5(j) & +t6*outm6(j)+t7*outm7(j) enddo out2=6*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) & +t3*outm3(j)+t4*outm4(j)+t5*outm5(j) & +t6*outm6(j) enddo out2=5*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) & +t3*outm3(j)+t4*outm4(j)+t5*outm5(j) enddo out2=4*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) & +t3*outm3(j)+t4*outm4(j) enddo out2=3*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outm1(j)+t2*outm2(j) & +t3*outm3(j) enddo out2=2*stride do j=1+out2,area+out2 out(j)= t0*out(j)+t1*outm1(j)+t2*outm2(j) enddo out2=1*stride do j=1+out2,area+out2 out(j)= t0*out(j)+t1*outm1(j) enddo out2=0*stride do j=1+out2,area+out2 out(j)= t0*out(j) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) t6=tf(7) t7=tf(8) t8=tf(9) do i=0,n-1,2 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) & +t4*outp4(j)+t5*outp5(j)+t6*outp6(j)+t7*outp7(j) & +t8*outp8(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c computing data ONLY on ghostpoints, AND folding in input of a c different data set (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c in: input/output array for in place computation c inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of out sets c stride: separation between out data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: separation between in data sets (stride=&in-&inm1) c fptl, lptl: first and last data components actually present in input c tr: reverse-time filter coefficients c tf: forward-time filter coefficients c eflag: mod(,2) of present line (tells if special ghost handling needed) c zcflag: =1 means we can count on zero ghosts and need special handling subroutine pconvsetsi3old( $ out,outm3,outm2,outm1,outp1,outp2,outp3, $ in,inm3,inm2,inm1,inp1,inp2,inp3, $ area,stride,n,npad, $ stridel,fptl,lptl, $ tr,tf,eflag,zcflag) implicit none double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) double precision in(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) integer stride,area,n,npad integer stridel,fptl,lptl double precision tr(*),tf(*) integer eflag,zcflag double precision t0,t1,t2,t3 integer i,j,out2 c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) c First, make sure pad is zero for forward phase. c For pretty code the last number would be npad-1, but technically c npad-2 (for extra speed) is all that is needed. do i=n,npad-2 call zeroline3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) enddo call zeroline3(out,in,inm1,inm2,inm3,n-1, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) if (eflag.eq.0 .and. zcflag.eq.1) then c For Idag (zcflag=1) we are coming from real space, and so can't c depend on the ghostpoints being zero. Thus, on the even planes (eflag), c we must skip over the data in the appropriate lines for the even c points. This is what the line??[eo] routines do. c For the UNUSED lines below, we point to "inm3" to avoid seg faults c because we know "inm3" will always be in range (unless the box is c ridiculously small). Although multiplied by zero, the code still c accesses the dummy data! call lineioe3(out,inm3,inm3,inm3,inm3,n-2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,t3) call lineioo3(out,inm3,inm3,inm2,inm3,n-3, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t2,t3) call lineioe3(out,inm3,inm1,inm2,inm3,n-4, $ stridel,stride,fptl,lptl,area, $ 0.d0,t1,t2,t3) do i=n-5,5,-1 if (mod(i,2).eq.0) then call lineioe3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) else call lineioo3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) endif enddo call lineioe3(out,in,inm1,inm2,in,4, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,0.d0) call lineioo3(out,in,inm1,in,in,3, $ stridel,stride,fptl,lptl,area, $ t0,t1,0.d0,0.d0) call lineioe3(out,in,in,in,in,2, $ stridel,stride,fptl,lptl,area, $ t0,0.d0,0.d0,0.d0) else c For the UNUSED lines below, we point to "inm3" to avoid seg faults c because we know "inm3" will always be in range (unless the box is c ridiculously small). Although multiplied by zero, the code still c accesses the dummy data! call lineio3(out,inm3,inm3,inm3,inm3,n-2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,t3) call lineio3(out,inm3,inm3,inm2,inm3,n-3, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t2,t3) call lineio3(out,inm3,inm1,inm2,inm3,n-4, $ stridel,stride,fptl,lptl,area, $ 0.d0,t1,t2,t3) do i=n-5,5,-1 call lineio3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) enddo c For the UNUSED lines below, we must point to "in" to avoid seg c faults. On this end of the convolution, we know that "in" will c always be in range. call lineio3(out,in,inm1,inm2,in,4, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,0.d0) call lineio3(out,in,inm1,in,in,3, $ stridel,stride,fptl,lptl,area, $ t0,t1,0.d0,0.d0) call lineio3(out,in,in,in,in,2, $ stridel,stride,fptl,lptl,area, $ t0,0.d0,0.d0,0.d0) endif call zeroline3(out,in,inm1,inm2,inm3,1, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) call zeroline3(out,in,inm1,inm2,inm3,0, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=0,n-1,2 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c Back and forth, linear convolution of sets of points using padding c computing data ONLY on ghostpoints, AND folding in input of a c different data set (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c in: input/output array for in place computation c inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of out sets c stride: separation between out data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: separation between in data sets (stride=&in-&inm1) c fptl, lptl: first and last data components actually present in input c tr: reverse-time filter coefficients c tf: forward-time filter coefficients c eflag: mod(,2) of present line (tells if special ghost handling needed) c zcflag: =1 means we can count on zero ghosts and need special handling subroutine pconvsetsi3new(out, $ in,inm3,inm2,inm1,inp1,inp2,inp3, $ area,stride,n,npad, $ stridel,fptl,lptl, $ filt,eflag,zcflag) implicit none double precision out(*) double precision in(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) integer stride,area,n,npad integer stridel,fptl,lptl double precision filt(*) integer eflag,zcflag double precision t0,t1,t2,t3 integer i c Filter coefficients t0=filt(1) t1=filt(2) t2=filt(3) t3=filt(4) do i=0,npad-2 call zeroline3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) enddo c For the UNUSED lines below, we point to "inm3" or "inp3" to avoid c seg faults because we know "inm3" or "inp3" will always be in range c (unless the box is ridiculously small). Although multiplied c by zero, the code still accesses the dummy data! if (eflag.eq.0 .and. zcflag.eq.1) then call lineio3full(out,inp3,inp3,inp3,inp3,inp3,inp3,inp3,0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,t3) call lineio3full(out,inp3,inp3,inp3,in,inp1,inp2,inp3,2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,t0,t1,t2,t3) call lineio3full(out,inp3,inp3,inm1,in,inp1,inp2,inp3,4, $ stridel,stride,fptl,lptl,area, $ 0.d0,t2,t1,t0,t1,t2,t3) do i=6,n-8,2 call lineio3full(out,inm3,inm2,inm1,in,inp1,inp2,inp3,i, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,t0,t1,t2,t3) enddo call lineio3full(out,inm3,inm2,inm1, $ in,inp1,inm3,inm3,n-6, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,t0,t1,0.d0,0.d0) call lineio3full(out,inm3,inm2,inm1, $ inm3,inm3,inm3,inm3,n-4, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,0.d0,0.d0,0.d0,0.d0) call lineio3full(out,inm3,inm3,inm3, $ inm3,inm3,inm3,inm3,n-2, $ stridel,stride,fptl,lptl,area, $ t3,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0) else call lineio3fullacc(out,inp3,inp3,inp3,inp3,inp3,inp3,inp3,0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,t3) call lineio3fullacc(out,inp3,inp3,inp3,in,inp1,inp2,inp3,2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,t0,t1,t2,t3) call lineio3fullacc(out,inp3,inp3,inm1,in,inp1,inp2,inp3,4, $ stridel,stride,fptl,lptl,area, $ 0.d0,t2,t1,t0,t1,t2,t3) do i=6,n-8,2 call lineio3fullacc(out,inm3,inm2,inm1,in,inp1,inp2,inp3,i, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,t0,t1,t2,t3) enddo call lineio3fullacc(out,inm3,inm2,inm1, $ in,inp1,inm3,inm3,n-6, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,t0,t1,0.d0,0.d0) call lineio3fullacc(out,inm3,inm2,inm1, $ inm3,inm3,inm3,inm3,n-4, $ stridel,stride,fptl,lptl,area, $ t3,t2,t1,0.d0,0.d0,0.d0,0.d0) call lineio3fullacc(out,inm3,inm3,inm3, $ inm3,inm3,inm3,inm3,n-2, $ stridel,stride,fptl,lptl,area, $ t3,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0) endif return end c Back and forth, linear convolution of sets of points using padding c computing data ONLY on ghostpoints, AND folding in input of a c different data set (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c in: input/output array for in place computation c inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of out sets c stride: separation between out data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: separation between in data sets (stride=&in-&inm1) c fptl, lptl: first and last data components actually present in input c tr: reverse-time filter coefficients c tf: forward-time filter coefficients c eflag: mod(,2) of present line (tells if special ghost handling needed) c zcflag: =1 means we can count on zero ghosts and need special handling subroutine pconvsetsi8new(out, $ in, $ inm8,inm7,inm6,inm5,inm4,inm3,inm2,inm1, $ inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ area,stride,n,npad, $ stridel,fptl,lptl, $ filt,eflag,zcflag) implicit none double precision out(*) double precision in(*) double precision inm8(*),inm7(*),inm6(*) double precision inm5(*),inm4(*),inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*),inp4(*) double precision inp5(*),inp6(*),inp7(*),inp8(*) integer stride,area,n,npad integer stridel,fptl,lptl double precision filt(*) integer eflag,zcflag double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i c Filter coefficients t0=filt(1) t1=filt(2) t2=filt(3) t3=filt(4) t4=filt(5) t5=filt(6) t6=filt(7) t7=filt(8) t8=filt(9) c For the UNUSED lines below, we point to "inm3" or "inp3" to avoid c seg faults because we know "inm3" or "inp3" will always be in range c (unless the box is ridiculously small). Although multiplied c by zero, the code still accesses the dummy data! if (t0 /= 0.0) then c symmetric filter if (2.GT.3) then c if (eflag.eq.0 .and. zcflag.eq.1) then call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ 0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t3) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp6,inp7,inp8, $ 2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp4,inp5,inp6,inp7,inp8, $ 4, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 6, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 8, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 10, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 12, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 14, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) do i=16,n-18,2 call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ i, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) enddo call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inm8, $ n-16, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inm8,inm8,inm8, $ n-14, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inm8,inm8,inm8,inm8,inm8, $ n-12, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-10, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-8, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-6, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-4, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-2, $ stridel,stride,fptl,lptl,area, $ t8,t7,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) else call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ 0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp6,inp7,inp8, $ 2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp4,inp5,inp6,inp7,inp8, $ 4, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 6, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 8, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 10, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 12, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 14, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) do i=16,n-18,2 call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ i, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) enddo call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inm8, $ n-16, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inm8,inm8,inm8, $ n-14, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inm8,inm8,inm8,inm8,inm8, $ n-12, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-10, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-8, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-6, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-4, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-2, $ stridel,stride,fptl,lptl,area, $ t8,t7,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) endif c end of symmetric part else c antisymmetric filter if (2.GT.3) then c if (eflag.eq.0 .and. zcflag.eq.1) then call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ 0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t3) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp6,inp7,inp8, $ 2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp4,inp5,inp6,inp7,inp8, $ 4, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 6, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 8, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inp8,inp8,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 10, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inp8,inp8,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 12, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8full(out,inp8,inp8,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 14, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) do i=16,n-18,2 call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ i, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) enddo call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inm8, $ n-16, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,t6,t7,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inm8,inm8,inm8, $ n-14, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ t4,t5,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inm8,inm8,inm8,inm8,inm8, $ n-12, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,t2,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-10, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,t0,t1,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-8, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ t2,t1,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm4,inm3,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-6, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,t4,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm6,inm5,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-4, $ stridel,stride,fptl,lptl,area, $ t8,t7,t6,t5,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8full(out,inm8,inm7,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-2, $ stridel,stride,fptl,lptl,area, $ t8,t7,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) c Is any of the above code ever called? Is 2 ever greater than 3? else call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ 0, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp8,inp8,inp6,inp7,inp8, $ 2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp8,inp8,inp4,inp5,inp6,inp7,inp8, $ 4, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,inp8,inp8,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 6, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inp8, $ inp8,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 8, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inp8,inp8,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 10, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inp8,inp8,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 12, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) call lineio8fullacc(out,inp8,inp8,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ 14, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) do i=16,n-18,2 call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, $ i, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,t6,t7,t8) enddo call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inm8, $ n-16, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,t6,t7,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inp4,inp5,inm8,inm8,inm8, $ n-14, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ t4,t5,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inp2,inp3,inm8,inm8,inm8,inm8,inm8, $ n-12, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,t2,t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,in,inp1,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-10, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,t0,t1,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm2, $ inm1,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-8, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ -1.*t2,-1.*t1,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4,inm3,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-6, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,-1.*t4,-1.*t3, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm6,inm5,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-4, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,-1.*t6,-1.*t5,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8fullacc(out,inm8,inm7,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8,inm8, $ n-2, $ stridel,stride,fptl,lptl,area, $ -1.*t8,-1.*t7,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0, $ 0.d0,0.d0,0.d0,0.d0,0.d0) endif c end of antisymmetric part endif return end c Back and forth, linear convolution of sets of points using padding c computing data ONLY on ghostpoints, AND folding in input of a c different data set (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c in: input/output array for in place computation c inm3,inm2,inm1,inp1,inp2,inp3: pointers in input array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of out sets c stride: separation between out data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: separation between in data sets (stride=&in-&inm1) c fptl, lptl: first and last data components actually present in input c tr: reverse-time filter coefficients c tf: forward-time filter coefficients c eflag: mod(,2) of present line (tells if special ghost handling needed) c zcflag: =1 means we can count on zero ghosts and need special handling subroutine pconvsetsi8( $ out,outm8,outm7,outm6,outm5,outm4, $ outm3,outm2,outm1,outp1,outp2,outp3, $ outp4,outp5,outp6,outp7,outp8, $ in,inm8,inm7,inm6,inm5,inm4, $ inm3,inm2,inm1,inp1,inp2,inp3, $ inp4,inp5,inp6,inp7,inp8, $ area,stride,n,npad, $ stridel,fptl,lptl, $ tr,tf,eflag,zcflag) implicit none double precision out(*) double precision outm8(*),outm7(*) double precision outm6(*),outm5(*),outm4(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) double precision outp4(*),outp5(*),outp6(*) double precision outp7(*),outp8(*) double precision in(*) double precision inm8(*),inm7(*) double precision inm6(*),inm5(*),inm4(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) double precision inp4(*),inp5(*),inp6(*) double precision inp7(*),inp8(*) integer stride,area,n,npad integer stridel,fptl,lptl double precision tr(*),tf(*) integer eflag,zcflag double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2 c Reverse phase t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) t6=tr(7) t7=tr(8) t8=tr(9) c First, make sure pad is zero for forward phase. c For pretty code the last number would be npad-1, but technically c npad-2 (for extra speed) is all that is needed. do i=n,npad-2 call zeroline8(out,in,inm1,inm2,inm3,inm4, $ inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) enddo c call zeroline8(out,in,inm1,inm2,inm3, c $ inm4,inm5,inm6,inm7,inm8,n-1, c $ stridel,stride,fptl,lptl,area, c $ t0,t1,t2,t3,t4,t5,t6,t7,t8) c CALL flush(6) c For the UNUSED lines below, we point to "inm8" to avoid seg faults c because we know "inm8" will always be in range (unless the box is c ridiculously small). Although multiplied by zero, the code still c accesses the dummy data! call lineio8(out,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm8,inm8,n-1, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,t8) call lineio8(out,inm8,inm8,inm8,inm8,inm8, $ inm8,inm8,inm7,inm8,n-2, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,t7,t8) call lineio8(out,inm8,inm8,inm8,inm8,inm8, $ inm8,inm6,inm7,inm8,n-3, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,t6,t7,t8) call lineio8(out,inm8,inm8,inm8,inm8,inm8, $ inm5,inm6,inm7,inm8,n-4, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,0.d0,t5,t6,t7,t8) call lineio8(out,inm8,inm8,inm8,inm8,inm4, $ inm5,inm6,inm7,inm8,n-5, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,0.d0,t4,t5,t6,t7,t8) call lineio8(out,inm8,inm8,inm8,inm3,inm4, $ inm5,inm6,inm7,inm8,n-6, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,0.d0,t3,t4,t5,t6,t7,t8) call lineio8(out,inm8,inm8,inm2,inm3,inm4, $ inm5,inm6,inm7,inm8,n-7, $ stridel,stride,fptl,lptl,area, $ 0.d0,0.d0,t2,t3,t4,t5,t6,t7,t8) call lineio8(out,inm8,inm1,inm2,inm3,inm4, $ inm5,inm6,inm7,inm8,n-8, $ stridel,stride,fptl,lptl,area, $ 0.d0,t1,t2,t3,t4,t5,t6,t7,t8) do i=n-9,16,-1 call lineio8(out,in,inm1,inm2,inm3, $ inm4,inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) enddo c For the UNUSED lines below, we must point to "in" to avoid seg c faults. On this end of the convolution, we know that "in" will c always be in range. call lineio8(out,in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,in,15, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,0.d0) call lineio8(out,in,inm1,inm2,inm3,inm4,inm5,inm6,in,in,14, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,0.d0,0.d0) call lineio8(out,in,inm1,inm2,inm3,inm4,inm5,in,in,in,13, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,0.d0,0.d0,0.d0) call lineio8(out,in,inm1,inm2,inm3,inm4,in,in,in,in,12, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,0.d0,0.d0,0.d0,0.d0) call lineio8(out,in,inm1,inm2,inm3,in,in,in,in,in,11, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8(out,in,inm1,inm2,in,in,in,in,in,in,10, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8(out,in,inm1,in,in,in,in,in,in,in,9, $ stridel,stride,fptl,lptl,area, $ t0,t1,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0) call lineio8(out,in,in,in,in,in,in,in,in,in,8, $ stridel,stride,fptl,lptl,area, $ t0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0,0.d0) do i=7,0,-1 call zeroline8(out,in,inm1,inm2,inm3, $ inm4,inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) enddo c Forward phase (padding now ensures room!) t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) t6=tf(7) t7=tf(8) t8=tf(9) do i=0,n-1,2 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) & +t4*outp4(j)+t5*outp5(j)+t6*outp6(j)+t7*outp7(j) & +t8*outp8(j) enddo enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineio3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points do j=1+fptl+out2,1+lptl+out2 out(j+dout2)= & t0*in(j)+t1*inm1(j)+t2*inm2(j)+t3*inm3(j) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineio3full(out,inm3,inm2,inm1,in,inp1,inp2,inp3,i, $ stridel,stride,fptl,lptl,area, $ tl3,tl2,tl1,t0,tr1,tr2,tr3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inp1(*),inp2(*),inp3(*) integer i,stridel,stride,fptl,lptl,area double precision tl3,tl2,tl1,t0,tr1,tr2,tr3 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) do j=1+out2+2,lptl+out2+1,2 c do j=1+out2+1,lptl+out2+1 out(j+dout2)=tr3*inp3(j)+tr1*inp1(j)+ & tl1*inm1(j)+tl3*inm3(j) out(j+dout2+1)=tr3*inp3(j+1)+tr1*inp1(j+1)+t0*in(j+1)+ & tl1*inm1(j+1)+tl3*inm3(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineio8full(out,inm8,inm7,inm6,inm5,inm4, $ inm3,inm2,inm1,in,inp1,inp2,inp3, $ inp4,inp5,inp6,inp7,inp8, $ i, $ stridel,stride,fptl,lptl,area, $ tl8,tl7,tl6,tl5,tl4,tl3,tl2,tl1, $ t0,tr1,tr2,tr3,tr4,tr5,tr6,tr7,tr8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) double precision inp1(*),inp2(*),inp3(*),inp4(*) double precision inp5(*),inp6(*),inp7(*),inp8(*) integer i,stridel,stride,fptl,lptl,area double precision tl8,tl7,tl6,tl5,tl4,tl3,tl2,tl1 double precision t0,tr1,tr2,tr3,tr4,tr5,tr6,tr7,tr8 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) do j=1+out2+2,lptl+out2+1,2 out(j+dout2)=tr8*inp8(j)+tr7*inp7(j)+tr6*inp6(j) & +tr5*inp5(j)+tr4*inp4(j)+tr3*inp3(j)+tr2*inp2(j) & +tr1*inp1(j)+tl1*inm1(j)+tl2*inm2(j)+tl3*inm3(j) & +tl3*inm3(j)+tl4*inm4(j)+tl5*inm4(j)+tl6*inm6(j) & +tl7*inm7(j)+tl8*inm8(j) out(j+dout2+1)=tr8*inp8(j+1)+tr7*inp7(j+1)+tr6*inp6(j+1) & +tr5*inp5(j+1)+tr4*inp4(j+1)+tr3*inp3(j+1)+tr2*inp2(j+1) & +tr1*inp1(j+1)+tl1*inm1(j+1)+tl2*inm2(j+1)+tl3*inm3(j+1) & +tl3*inm3(j+1)+tl4*inm4(j+1)+tl5*inm4(j+1)+tl6*inm6(j+1) & +tl7*inm7(j+1)+tl8*inm8(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineio3fullacc(out,inm3,inm2,inm1,in,inp1,inp2,inp3,i, $ stridel,stride,fptl,lptl,area, $ tl3,tl2,tl1,t0,tr1,tr2,tr3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inp1(*),inp2(*),inp3(*) integer i,stridel,stride,fptl,lptl,area double precision tl3,tl2,tl1,t0,tr1,tr2,tr3 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) do j=1+out2+fptl,lptl+out2+1 c do j=1+out2+1,lptl+out2+1 out(j+dout2)=tr3*inp3(j)+tr1*inp1(j)+t0*in(j)+ & tl1*inm1(j)+tl3*inm3(j) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineio8fullacc(out,inm8,inm7,inm6,inm5,inm4, $ inm3,inm2,inm1,in,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8,i, $ stridel,stride,fptl,lptl,area, $ tl8,tl7,tl6,tl5,tl4,tl3,tl2,tl1, $ t0,tr1,tr2,tr3,tr4,tr5,tr6,tr7,tr8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) double precision inp1(*),inp2(*),inp3(*),inp4(*) double precision inp5(*),inp6(*),inp7(*),inp8(*) integer i,stridel,stride,fptl,lptl,area double precision tl8,tl7,tl6,tl5,tl4,tl3,tl2,tl1 double precision t0,tr1,tr2,tr3,tr4,tr5,tr6,tr7,tr8 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo do j=1+out2+fptl,lptl+out2+1 c do j=1+out2+2,lptl+out2+1 c do j=1+out2+1,lptl+out2+1 out(j+dout2)=tr8*inp8(j)+tr7*inp7(j)+tr6*inp6(j)+tr5*inp5(j) & +tr4*inp4(j)+tr3*inp3(j)+tr2*inp2(j)+tr1*inp1(j)+t0*in(j) & +tl1*inm1(j)+tl2*inm2(j)+tl3*inm3(j)+tl4*inm4(j) & +tl5*inm5(j)+tl6*inm6(j)+tl7*inm7(j)+tl8*inm8(j) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single data set ("line") during reverse c phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients subroutine lineio8(out,in,inm1,inm2,inm3, $ inm4,inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points do j=1+fptl+out2,1+lptl+out2 out(j+dout2)= & t0*in(j)+t1*inm1(j)+t2*inm2(j)+t3*inm3(j) & +t4*inm4(j)+t5*inm5(j)+t6*inm6(j)+t7*inm7(j) & +t8*inm8(j) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Zeros an output line (OUTPUT) c consistent calling sequence with lineio c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine zeroline3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3 integer out2,j out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo return end c Zeros an output line (OUTPUT) c consistent calling sequence with lineio c c out: output set c i: index of line to be done c stride: stride of output c area: full area to be generated subroutine zeroline8(out,in,inm1,inm2,inm3, $ inm4,inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer out2,j out2=i*stride do j=1+out2,area+out2 out(j)=0. enddo return end c Computes convolution for a single "even" data set ("line") during c reverse phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineioe3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points: c Ignore any even input data on even input lines! do j=1+fptl+out2,1+lptl+out2,2 out(j+dout2)= & t1*inm1(j) +t3*inm3(j) out(j+1+dout2)= & t0*in(j+1)+t1*inm1(j+1)+t2*inm2(j+1)+t3*inm3(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single "even" data set ("line") during c reverse phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients subroutine lineioe8(out,in,inm1,inm2,inm3,inm4,inm5,inm6, $ inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points: c Ignore any even input data on even input lines! do j=1+fptl+out2,1+lptl+out2,2 out(j+dout2)= & t1*inm1(j) +t3*inm3(j) out(j+1+dout2)= & t0*in(j+1)+t1*inm1(j+1)+t2*inm2(j+1)+t3*inm3(j+1) & +t4*inm4(j+1)+t5*inm5(j+1)+t6*inm6(j+1)+t7*inm7(j+1) & +t8*inm8(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single "odd" data set ("line") during c reverse phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3: filter coefficients subroutine lineioo3(out,in,inm1,inm2,inm3,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points: c Ignore any even input data on even input lines! do j=1+fptl+out2,1+lptl+out2,2 out(j+dout2)= & t0*in(j) +t2*inm2(j) out(j+1+dout2)= & t0*in(j+1)+t1*inm1(j+1)+t2*inm2(j+1)+t3*inm3(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Computes convolution for a single "odd" data set ("line") during c reverse phases (INPUT/OUTPUT) c c out: output set c in,inm1,inm2,inm3,inm4,inm5,inm6,inm7,inm8: input set c i: index of line to be done c stridel: stride of input c stride: stride of output c fptl: first point of input actually used c lptl: last point of input actually used c area: full area to be generated c t0,t1,t2,t3,t4,t5,t6,t7,t8: filter coefficients subroutine lineioo8(out,in,inm1,inm2,inm3, $ inm4,inm5,inm6,inm7,inm8,i, $ stridel,stride,fptl,lptl,area, $ t0,t1,t2,t3,t4,t5,t6,t7,t8) implicit none double precision out(*),in(*),inm1(*),inm2(*),inm3(*) double precision inm4(*),inm5(*),inm6(*),inm7(*),inm8(*) integer i,stridel,stride,fptl,lptl,area double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer out2,dout2,j out2=i*stridel dout2=i*(stride-stridel) c Zero out data before fptl do j=1+out2,fptl+out2 out(j+dout2)=0. enddo c Do the line for "active" points: c Ignore any even input data on even input lines! do j=1+fptl+out2,1+lptl+out2,2 out(j+dout2)= & t0*in(j) +t2*inm2(j) out(j+1+dout2)= & t0*in(j+1)+t1*inm1(j+1)+t2*inm2(j+1)+t3*inm3(j+1) & +t4*inm4(j+1)+t5*inm5(j+1)+t6*inm6(j+1) & +t8*inm8(j+1) enddo c Zero out data after lptl do j=1+lptl+out2+1,area+out2 out(j+dout2)=0. enddo return end c Back and forth, linear convolution of sets of points for +/- 3 c onvolutions using padding, and ignoring input data on non-ghost c points (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsp3(out,in,inm3,inm2,inm1,inp1,inp2,inp3, & area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision in(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t3 integer i,j,out2 c Forward phase (padding now ensures room!) c No change here relative to convsets t0=tf(1) t1=tf(2) t3=tf(4) do i=3,n-5,2 out2=i*stride do j=1+out2,area+out2 out(j)=t3*inm3(j)+t1*inm1(j)+ & t1*inp1(j)+t3*inp3(j) enddo enddo do i=2,n-5,2 out2=i*stride do j=1+out2,area+out2 out(j)=in(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 8 c onvolutions using padding, and ignoring input data on non-ghost c points (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsp8(out, $ outm8,outm7,outm6,outm5,outm4,outm3,outm2,outm1, $ outp1,outp2,outp3,outp4,outp5,outp6,outp7,outp8, $ area,stride,n,npad,tr,tf) implicit none double precision out(*) double precision outm8(*),outm7(*),outm6(*),outm5(*) double precision outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*) double precision outp5(*),outp6(*),outp7(*),outp8(*) integer stride,area,n,npad double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2 c Reverse direction c Relative to convsets, we wiped out the non-ghost references t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) t6=tr(7) t7=tr(8) t8=tr(9) c Pad region c out2=(npad-1)*stride c do j=1+out2,area+out2 c out(j)=0. c enddo c out2=(npad-2)*stride c do j=1+out2,area+out2 c out(j)=t8*outm8(j) c enddo c out2=(npad-3)*stride c do j=1+out2,area+out2 c out(j)=t7*outm7(j) c enddo c out2=(npad-4)*stride c do j=1+out2,area+out2 c out(j)=t6*outm6(j)+t8*outm8(j) c enddo c out2=(npad-5)*stride c do j=1+out2,area+out2 c out(j)=t5*outm5(j)+t7*outm7(j) c enddo c out2=(npad-6)*stride c do j=1+out2,area+out2 c out(j)=t4*outm4(j)+t6*outm6(j)+t8*outm8(j) c enddo c out2=(npad-7)*stride c do j=1+out2,area+out2 c out(j)=t3*outm3(j)+t5*outm5(j)+t7*outm7(j) c enddo c out2=(npad-8)*stride c do j=1+out2,area+out2 c out(j)=t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) c enddo c Main reverse phase exploiting zeros! do i=n-1,8,-2 out2=i*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j)+t7*outm7(j) enddo out2=(i-1)*stride do j=1+out2,area+out2 out(j)=t0*out(j) $ +t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo enddo c Fence posting c out2=7*stride c do j=1+out2,area+out2 c out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j)+t7*outm7(j) c enddo c out2=6*stride c do j=1+out2,area+out2 c out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j)+t6*outm6(j) c enddo c out2=5*stride c do j=1+out2,area+out2 c out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j) c enddo c out2=4*stride c do j=1+out2,area+out2 c out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j) c enddo c out2=3*stride c do j=1+out2,area+out2 c out(j)=t1*outm1(j)+t3*outm3(j) c enddo c out2=2*stride c do j=1+out2,area+out2 c out(j)=t0*out(j)+t2*outm2(j) c enddo c out2=1*stride c do j=1+out2,area+out2 c out(j)=t1*outm1(j) c enddo c out2=0*stride c do j=1+out2,area+out2 c out(j)=t0*out(j) c enddo c Forward phase (padding now ensures room!) c No change here relative to convsets t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) t6=tf(7) t7=tf(8) t8=tf(9) c Need to verify that this is OK..... c do i=0,n-1 do i=8,n-8 out2=i*stride do j=1+out2,area+out2 out(j)= & t0*out(j) $ +t1*outp1(j)+t2*outp2(j)+t3*outp3(j)+t4*outp4(j) $ +t5*outp5(j)+t6*outp6(j)+t7*outp7(j)+t8*outp8(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 3 c onvolutions using padding, and ignoring input data on non-ghost c points (IN PLACE) c c out: input/output array for in place computation c outm3,outm2,outm1,outp1,outp2,outp3: pointers in output array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetsp8full(out,in,inm8,inm7,inm6,inm5,inm4, & inm3,inm2,inm1,inp1,inp2,inp3,inp4,inp5,inp6,inp7,inp8, & area,stride,n,npad,filter) implicit none double precision out(*) double precision in(*) double precision inm8(*),inm7(*),inm6(*),inm5(*),inm4(*) double precision inm3(*),inm2(*),inm1(*) double precision inp1(*),inp2(*),inp3(*) double precision inp4(*),inp5(*),inp6(*),inp7(*),inp8(*) integer stride,area,n,npad double precision filter(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2 c Forward phase (padding now ensures room!) c No change here relative to convsets t0=filter(1) t1=filter(2) t2=filter(3) t3=filter(4) t4=filter(5) t5=filter(6) t6=filter(7) t7=filter(8) t8=filter(9) if (t0 /= 0.0) then c symmetric filter c thin func do i=9,n-8,2 out2=i*stride do j=1+out2,area+out2 out(j)=t7*inm7(j)+t5*inm5(j)+t3*inm3(j)+t1*inm1(j)+ & t1*inp1(j)+t3*inp3(j)+t5*inp5(j)+t7*inp7(j) enddo enddo c fat func do i=8,n-8,2 out2=i*stride do j=1+out2,area+out2 out(j)=t8*inm8(j)+t6*inm6(j)+t4*inm4(j)+t2*inm2(j)+ & t0*in(j)+t2*inp2(j)+t4*inp4(j)+t6*inp6(j)+t8*inp8(j) enddo enddo c end of symmetric part else c antisymmetric filter c thin func do i=9,n-8,2 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t7*inm7(j)-t5*inm5(j)-t3*inm3(j)-t1*inm1(j)+ & t1*inp1(j)+t3*inp3(j)+t5*inp5(j)+t7*inp7(j) enddo enddo c fat func do i=8,n-8,2 out2=i*stride do j=1+out2,area+out2 out(j)=-1.*t8*inm8(j)-t6*inm6(j)-t4*inm4(j)-t2*inm2(j)+ & t0*in(j)+t2*inp2(j)+t4*inp4(j)+t6*inp6(j)+t8*inp8(j) enddo enddo c end of antisymmetric part endif return end c Back and forth, linear convolution of sets of points using padding, c ignoring input data on non-ghost points (IN/OUT) c c lower: output array c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspo3(lower,out, $ outm3,outm2,outm1,outp1,outp2,outp3, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ tr,tf) implicit none double precision lower(*) double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision tr(*),tf(*) double precision t0,t1,t2,t3 integer i,j,out2,out2s,jl c Reverse direction c No change here relative to convsetsp t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) c Pad region out2=(npad-1)*stride do j=1+out2,area+out2 out(j)=0. enddo out2=(npad-2)*stride do j=1+out2,area+out2 out(j)= t3*outm3(j) enddo out2=(npad-3)*stride do j=1+out2,area+out2 out(j)= t2*outm2(j) enddo c Main reverse phase exploiting zeros! do i=n-1,3,-2 out2=i*stride do j=1+out2,area+out2 out(j)= $ t1*outm1(j) +t3*outm3(j) enddo out2=(i-1)*stride do j=1+out2,area+out2 out(j)= $ t0*out(j) +t2*outm2(j) enddo enddo c Fence posting out2=stride do j=1+out2,area+out2 out(j)= $ t1*outm1(j) enddo out2=0 do j=1+out2,area+out2 out(j)= $ t0*out(j) enddo c Forward phase (padding now ensures room!) c Relative to convsetsp, output now places into lower instead t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=fsetl,lsetl out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s c$$$ out(j)= c$$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) lower(jl)= $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c lower: output array c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspoadd3(lower,addon,out, $ outm3,outm2,outm1,outp1,outp2,outp3, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ tr,tf) implicit none double precision lower(*) double precision addon(*) double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision tr(*),tf(*) double precision t0,t1,t2,t3 integer i,j,out2,out2s,jl c Reverse direction c No change here relative to convsetsp t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) c Pad region out2=(npad-1)*stride do j=1+out2,area+out2 out(j)=0. enddo out2=(npad-2)*stride do j=1+out2,area+out2 out(j)= t3*outm3(j) enddo out2=(npad-3)*stride do j=1+out2,area+out2 out(j)= t2*outm2(j) enddo c Main reverse phase exploiting zeros! do i=n-1,3,-2 out2=i*stride do j=1+out2,area+out2 out(j)= $ t1*outm1(j) +t3*outm3(j) enddo out2=(i-1)*stride do j=1+out2,area+out2 out(j)= $ t0*out(j) +t2*outm2(j) enddo enddo c Fence posting out2=stride do j=1+out2,area+out2 out(j)= $ t1*outm1(j) enddo out2=0 do j=1+out2,area+out2 out(j)= $ t0*out(j) enddo c Forward phase (padding now ensures room!) c Relative to convsetsp, output now places into lower instead t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) do i=fsetl,lsetl out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s c$$$ out(j)= c$$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) lower(jl)=addon(jl)+ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 3 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspop3(lower,out, $ outm3,outm2,outm1,outp1,outp2,outp3, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ filt) implicit none double precision lower(*) double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision filt(*) double precision t0,t1,t3 integer i,j,out2,out2s,jl c Forward phase (padding now ensures room!) c Relative to convsetspo, only change is accumulation into lower t0=filt(1) t1=filt(2) t3=filt(4) do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl)+t3*outm3(j)+t1*outm1(j)+ $ t1*outp1(j)+t3*outp3(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl)+out(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 3 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspom3(lower,out, $ outm3,outm2,outm1,outp1,outp2,outp3, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ filt) implicit none double precision lower(*) double precision out(*) double precision outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision filt(*) double precision t0,t1,t3 integer i,j,out2,out2s,jl c Forward phase (padding now ensures room!) c Relative to convsetspo, only change is accumulation into lower t0=filt(1) t1=filt(2) t3=filt(4) do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl)-t3*outm3(j)-t1*outm1(j)- $ t1*outp1(j)-t3*outp3(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl)-out(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 8 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspop8(lower,out, $ outm8,outm7,outm6,outm5,outm4,outm3,outm2,outm1, $ outp1,outp2,outp3,outp4,outp5,outp6,outp7,outp8, & area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ tr,tf) implicit none double precision lower(*) double precision out(*) double precision outm8(*),outm7(*),outm6(*),outm5(*) double precision outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*) double precision outp5(*),outp6(*),outp7(*),outp8(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2,out2s,jl c Reverse direction c Relative to convsets, we wiped out the non-ghost references t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) t6=tr(7) t7=tr(8) t8=tr(9) c Pad region out2=(npad-1)*stride do j=1+out2,area+out2 out(j)=0. enddo out2=(npad-2)*stride do j=1+out2,area+out2 out(j)=t8*outm8(j) enddo out2=(npad-3)*stride do j=1+out2,area+out2 out(j)=t7*outm7(j) enddo out2=(npad-4)*stride do j=1+out2,area+out2 out(j)=t6*outm6(j)+t8*outm8(j) enddo out2=(npad-5)*stride do j=1+out2,area+out2 out(j)=t5*outm5(j)+t7*outm7(j) enddo out2=(npad-6)*stride do j=1+out2,area+out2 out(j)=t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo out2=(npad-7)*stride do j=1+out2,area+out2 out(j)=t3*outm3(j)+t5*outm5(j)+t7*outm7(j) enddo out2=(npad-8)*stride do j=1+out2,area+out2 out(j)=t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo c Main reverse phase exploiting zeros! do i=n-1,8,-2 out2=i*stride do j=1+out2,area+out2 out(j)=t7*outm7(j)+t5*outm5(j)+t3*outm3(j)+t1*outm1(j) enddo out2=(i-1)*stride do j=1+out2,area+out2 out(j)=t0*out(j) $ +t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo enddo c Fence posting out2=7*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j)+t7*outm7(j) enddo out2=6*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j)+t6*outm6(j) enddo out2=5*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j) enddo out2=4*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j) enddo out2=3*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j) enddo out2=2*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j) enddo out2=1*stride do j=1+out2,area+out2 out(j)=t1*outm1(j) enddo out2=0*stride do j=1+out2,area+out2 out(j)=t0*out(j) enddo c Forward phase (padding now ensures room!) c No change here relative to convsets t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) t6=tf(7) t7=tf(8) t8=tf(9) do i=fsetl,lsetl out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s c$$$ out(j)= c$$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) lower(jl)=lower(jl)+ & t0*out(j) $ +t1*outp1(j)+t2*outp2(j)+t3*outp3(j)+t4*outp4(j) $ +t5*outp5(j)+t6*outp6(j)+t7*outp7(j)+t8*outp8(j) enddo enddo return end c Back and forth, linear convolution of sets of points for +/- 3 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspop8full(lower,out, $ outm8,outm7,outm6,outm5,outm4,outm3,outm2,outm1, $ outp1,outp2,outp3,outp4,outp5,outp6,outp7,outp8, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ filt) implicit none double precision lower(*) double precision out(*) double precision outm8(*),outm7(*),outm6(*),outm5(*) double precision outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*) double precision outp5(*),outp6(*),outp7(*),outp8(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision filt(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2,out2s,jl c Forward phase (padding now ensures room!) c Relative to convsetspo, only change is accumulation into lower t0=filt(1) t1=filt(2) t2=filt(3) t3=filt(4) t4=filt(5) t5=filt(6) t6=filt(7) t7=filt(8) t8=filt(9) if (t0 /= 0.0) then c symmetric filter do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl) $ +t7*outm7(j)+t5*outm5(j)+t3*outm3(j)+t1*outm1(j) $ +t1*outp1(j)+t3*outp3(j)+t5*outp5(j)+t7*outp7(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl) $ +t8*outm8(j)+t6*outm6(j)+t4*outm4(j)+t2*outm2(j) $ +t0*out(j)+t2*outp2(j)+t4*outp4(j)+t6*outp6(j) $ +t8*outp8(j) enddo enddo c end of symmetric part else c antisymmetric filter do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl) $ -t7*outm7(j)-t5*outm5(j)-t3*outm3(j)-t1*outm1(j) $ +t1*outp1(j)+t3*outp3(j)+t5*outp5(j)+t7*outp7(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=lower(jl) $ -t8*outm8(j)-t6*outm6(j)-t4*outm4(j)-t2*outm2(j) $ +t0*out(j)+t2*outp2(j)+t4*outp4(j)+t6*outp6(j) $ +t8*outp8(j) enddo enddo c end of antisymmetric part endif return end c Back and forth, linear convolution of sets of points for +/- 3 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspo8full(lower,out, $ outm8,outm7,outm6,outm5,outm4,outm3,outm2,outm1, $ outp1,outp2,outp3,outp4,outp5,outp6,outp7,outp8, $ area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ filt) implicit none double precision lower(*) double precision out(*) double precision outm8(*),outm7(*),outm6(*),outm5(*) double precision outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*) double precision outp5(*),outp6(*),outp7(*),outp8(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision filt(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2,out2s,jl c Forward phase (padding now ensures room!) c Relative to convsetspo, only change is accumulation into lower t0=filt(1) t1=filt(2) t2=filt(3) t3=filt(4) t4=filt(5) t5=filt(6) t6=filt(7) t7=filt(8) t8=filt(9) if (t0 /= 0.0) then c symmetric filter do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=t7*outm7(j)+t5*outm5(j)+t3*outm3(j)+t1*outm1(j) $ +t1*outp1(j)+t3*outp3(j)+t5*outp5(j)+t7*outp7(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=t8*outm8(j)+t6*outm6(j)+t4*outm4(j)+t2*outm2(j) $ +t0*out(j)+t2*outp2(j)+t4*outp4(j)+t6*outp6(j) $ +t8*outp8(j) enddo enddo c end of symmetric part else c antisymmetric filter do i=fsetl+1,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=-1.*t7*outm7(j)-t5*outm5(j)-t3*outm3(j) $ -t1*outm1(j)+t1*outp1(j)+t3*outp3(j)+t5*outp5(j) $ +t7*outp7(j) enddo enddo do i=fsetl,lsetl,2 out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s lower(jl)=-1.*t8*outm8(j)-t6*outm6(j)-t4*outm4(j) $ -t2*outm2(j)+t0*out(j)+t2*outp2(j)+t4*outp4(j) $ +t6*outp6(j)+t8*outp8(j) enddo enddo c end of antisymmetric part endif return end c Back and forth, linear convolution of sets of points for +/- 8 c filters using padding, ignoring input data on non-ghost points c (IN/+=OUT) c c lower: output array for accumulation c outm3,outm2,outm1,outp1,outp2,outp3: pointers to input work array to c start of data sets, 3 previous, 2 previous, 1 previous, 1 forward, c 2 forward, 3 forward c area: total size of each set c stride: separation between data sets (stride=&out-&outm1) c n: actual number of data sets being convolved c npad: size of padded region needed for "leakage" of foward filters c stridel: stride on output c fptl,lptl: first,last actual points in output set c fsetl,lsetl: first,last actual sets in output c tr: reverse-time filter coefficients c tf: forward-time filter coefficients subroutine convsetspo8(lower,out, $ outm8,outm7,outm6,outm5,outm4,outm3,outm2,outm1, $ outp1,outp2,outp3,outp4,outp5,outp6,outp7,outp8, & area,stride,n,npad, $ stridel,fptl,lptl,fsetl,lsetl, $ tr,tf) implicit none double precision lower(*) double precision out(*) double precision outm8(*),outm7(*),outm6(*),outm5(*) double precision outm4(*),outm3(*),outm2(*),outm1(*) double precision outp1(*),outp2(*),outp3(*),outp4(*) double precision outp5(*),outp6(*),outp7(*),outp8(*) integer stride,area,n,npad integer stridel,fptl,lptl,fsetl,lsetl double precision tr(*),tf(*) double precision t0,t1,t2,t3,t4,t5,t6,t7,t8 integer i,j,out2,out2s,jl c Reverse direction c Relative to convsets, we wiped out the non-ghost references t0=tr(1) t1=tr(2) t2=tr(3) t3=tr(4) t4=tr(5) t5=tr(6) t6=tr(7) t7=tr(8) t8=tr(9) c Pad region out2=(npad-1)*stride do j=1+out2,area+out2 out(j)=0. enddo out2=(npad-2)*stride do j=1+out2,area+out2 out(j)=t8*outm8(j) enddo out2=(npad-3)*stride do j=1+out2,area+out2 out(j)=t7*outm7(j) enddo out2=(npad-4)*stride do j=1+out2,area+out2 out(j)=t6*outm6(j)+t8*outm8(j) enddo out2=(npad-5)*stride do j=1+out2,area+out2 out(j)=t5*outm5(j)+t7*outm7(j) enddo out2=(npad-6)*stride do j=1+out2,area+out2 out(j)=t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo out2=(npad-7)*stride do j=1+out2,area+out2 out(j)=t3*outm3(j)+t5*outm5(j)+t7*outm7(j) enddo out2=(npad-8)*stride do j=1+out2,area+out2 out(j)=t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo c Main reverse phase exploiting zeros! do i=n-1,8,-2 out2=i*stride do j=1+out2,area+out2 out(j)=t7*outm7(j)+t5*outm5(j)+t3*outm3(j)+t1*outm1(j) enddo out2=(i-1)*stride do j=1+out2,area+out2 out(j)=t0*out(j) $ +t2*outm2(j)+t4*outm4(j)+t6*outm6(j)+t8*outm8(j) enddo enddo c Fence posting out2=7*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j)+t7*outm7(j) enddo out2=6*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j)+t6*outm6(j) enddo out2=5*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j)+t5*outm5(j) enddo out2=4*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j)+t4*outm4(j) enddo out2=3*stride do j=1+out2,area+out2 out(j)=t1*outm1(j)+t3*outm3(j) enddo out2=2*stride do j=1+out2,area+out2 out(j)=t0*out(j)+t2*outm2(j) enddo out2=1*stride do j=1+out2,area+out2 out(j)=t1*outm1(j) enddo out2=0*stride do j=1+out2,area+out2 out(j)=t0*out(j) enddo c Forward phase (padding now ensures room!) c No change here relative to convsets t0=tf(1) t1=tf(2) t2=tf(3) t3=tf(4) t4=tf(5) t5=tf(6) t6=tf(7) t7=tf(8) t8=tf(9) do i=fsetl,lsetl out2=i*stride out2s=i*(stridel-stride) do j=1+out2+fptl,1+out2+lptl jl=j+out2s c$$$ out(j)= c$$$ $ t0*out(j)+t1*outp1(j)+t2*outp2(j)+t3*outp3(j) lower(jl)= & t0*out(j) $ +t1*outp1(j)+t2*outp2(j)+t3*outp3(j)+t4*outp4(j) $ +t5*outp5(j)+t6*outp6(j)+t7*outp7(j)+t8*outp8(j) enddo enddo return end