!------------------------------------------------- ! Parallel (OMP) PM code ! 2015, 1997 Anatoly Klypin (aklypin@nmsu.edu) ! Astronomy Department, NMSU !------------------------------------------------- ! ! NGRID= number of cells in one dimension ! AEXPN= expansion parameter = 1/(1+z) Module Structures Real*4 :: AEXPN,AEXP0,AMPLT,ASTEP, box, & EKIN,EKIN1,EKIN2,AEU0, & TINTG,Nspecs,Nseed, & Om0,Oml0,hubble, & extras(100), & ENKIN,ENPOT Integer*4 :: iStep,Ngrid,Nrow Integer*8 :: Nparticles Real*4, Allocatable, Dimension(:,:,:) :: FI Real*4, Allocatable, Dimension(:) :: & XPAR,YPAR,ZPAR, & VX,VY,VZ Real*4, allocatable, Dimension(:) :: Xb,Yb,Zb,VXb,Vyb,Vzb Real*4 :: StepFactor =7.e-2 Real*4 :: da_max =8.e-3 Character*45 :: HEADER end Module Structures !------------------------------------------------- Program PMP use Structures ! Read data and open files OPEN(17,FILE='Result.log',STATUS='UNKNOWN',position='append') WRITE (*,'(A,$)') ' Enter number of steps to go => ' READ (*,*) Nsteps ! Make this number of steps CALL ReadDataPM ! CALL ReadDataGadget if(ISTEP==0) CALL WriteDataPM(1) IF(AEXPN.GE.1.)THEN ! change this if you need to run WRITE (*,*) ' Cannot run over a=1' ! beyond a=1 STOP ENDIF write(*,*) ' Start running:',Nparticles,Ngrid write(*,*) ' Start running:',ISTEP ! Main loop DO i=1,Nsteps Call Timing(0,-1) CALL DENSIT ! Define density CALL POTENTfft5 ! Define potential FactV = 1. ! normal step !if(i==1)FactV =2. ! make 1/2 step in V for the first step If(ASTEP (zyx) XPAR(i+ioff) = Zb(i) +0.5 YPAR(i+ioff) = Yb(i) ZPAR(i+ioff) = Xb(i) end Do read(1) (Xb(i),Yb(i),Zb(i),i=1,Nrecord) !-- read veloc !$OMP PARALLEL DO DEFAULT(SHARED) !--- rescale data Do i =1,Nrecord VX(i+ioff) = Zb(i) VY(i+ioff) = Yb(i) VZ(i+ioff) = Xb(i) end Do Read(1) ii !-- fake read of Id ioff = ioff + Nrecord Deallocate(Xb,Yb,Zb) !-- free buffer end Do !$OMP PARALLEL DO DEFAULT(SHARED) !--- rescale data Do i=1,Nparticles VX(i) = VX(i)/Vscale VY(i) = VY(i)/Vscale ! VZ(i) = VZ(i)/Vscale XPAR(i) = Xscale*XPAR(i)+1.5 ! scale (0-Box) --> (1-ngrid+1) YPAR(i) = Xscale*YPAR(i)+1.5 ! offset by 0.5 ZPAR(i) = Xscale*ZPAR(i)+1.5 If(INT(XPAR(i)).lt.1)XPAR(i) = XPAR(i) + ngrid If(INT(YPAR(i)).lt.1)YPAR(i) = YPAR(i) + ngrid If(INT(ZPAR(i)).lt.1)ZPAR(i) = ZPAR(i) + ngrid If(INT(XPAR(i)).ge.Ngrid+1)XPAR(i) = XPAR(i) - ngrid If(INT(YPAR(i)).ge.Ngrid+1)YPAR(i) = YPAR(i) - ngrid If(INT(ZPAR(i)).ge.Ngrid+1)ZPAR(i) = ZPAR(i) - ngrid EndDo Do i=1,Nparticles If(XPAR(i)<20..and.YPAR(i)<200..and.ZPAR(i)<200.)write(50,'(3ES14.5)')XPAR(i),YPAR(i),ZPAR(i) End Do Allocate (FI(NGRID,NGRID,NGRID)) end SUBROUTINE ReadDataGadget !--------------------------------------- ! Read PMfiles: PMcrd/crs SUBROUTINE ReadDataPM() !--------------------------------------- use Structures Character*80 :: Name Logical :: exst Integer*8 :: iCount,ii,ioff,ip Integer*8 :: Ngal,Nrecord,Jpage,NinPage ! ! Read data and open files Open (4,file ='PMcrd.DAT',form ='UNFORMATTED',status ='UNKNOWN') READ (4) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROW,NGRID,Nspecs,Nseed,Om0,Oml0,hubble, & Nparticles,extras WRITE (*,'(a,/10x,a,f8.4,4(a,i7))') HEADER, & ' a=',AEXPN, ' step= ',ISTEP, & ' Nrow= ', NROW, ' Ngrid=',NGRID WRITE (17,'(a,/10x,a,f8.4,4(a,i7))') HEADER, & ' a=',AEXPN, ' step= ',ISTEP, & ' Nrow= ', NROW, ' Ngrid=',NGRID close(4) Nrecord = 1024**2 Naccess = Nrecord*6 !*4 xR = NGRID +1 boxsize = extras(100) Box = boxsize Npages = (Nparticles-1)/Nrecord+1 ! number of records Nlast = Nparticles - (Npages-1)*Nrecord ! number of particles in last record Jpage = Nrecord write(*,'(a,i10)') ' NROW =',NROW write(*,'(a,i10)') ' Ngal =',Nparticles write(*,'(a,i10)') ' Ngrid =',Ngrid write(*,'(a,f10.1)') ' Box =',Box Allocate (Xb(Nrecord),Yb(Nrecord),Zb(Nrecord)) Allocate (VXb(Nrecord),VYb(Nrecord),VZb(Nrecord)) Allocate (XPAR(Nparticles),YPAR(Nparticles),ZPAR(Nparticles)) Allocate (VX(Nparticles),VY(Nparticles),VZ(Nparticles)) iCount = 0 ifile = 0 jj = 0 write(Name,'(a,i1.1,a)')'PMcrs',ifile,'.DAT' write(*,'(a)') TRIM(Name) INQUIRE(file=TRIM(Name),EXIST=exst) If(.not.exst)Stop ' File PMcrs0.DAT does not exist. Error' OPEN(UNIT=20,FILE=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) Do ii =1,Npages If(ii==Npages)Then NinPage = Nparticles -(ii-1)*JPAGE ! # particles in the current page Else NinPage = JPAGE EndIf jj = jj +1 If(ii<10.or.ii==Npages)write(*,'(3(a,i9))') ' Reading page= ',ii,' record =',jj,' NinPage= ',NinPage 10 Read(20,REC=jj,iostat=ierr) Xb,Yb,Zb,VXb,VYb,VZb If(ierr /=0)Then close(20) ifile = ifile +1 If(ifile<10)Then write(Name,'(a,i1.1,a,i4.4,a)')'PMcrs',ifile,'.DAT' Else write(Name,'(a,i2.2,a,i4.4,a)')'PMcrs',ifile,'.DAT' EndIf jj = 1 INQUIRE(file=TRIM(Name),EXIST=exst) If(.not.exst)Stop ' Error reading PMcrs files: Did not get all' Open(20,file=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) write(*,'(2i7,2a,3x,i9)') ii,ifile,' Open file = ',TRIM(Name),Ninpage go to 10 end If ioff = (ii-1)*JPAGE !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE (ip) Do ip =1,NinPage If(ip+ioff > Nparticles)STOP 'Attempt to read too many particles ' if(INT(Xb(ip))==Ngrid+1)Xb(ip)=Xb(ip)-1.e-3 if(INT(Yb(ip))==Ngrid+1)Yb(ip)=Yb(ip)-1.e-3 if(INT(Zb(ip))==Ngrid+1)Zb(ip)=Zb(ip)-1.e-3 if(INT(Xb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Xb(ip)),Xb(ip) if(INT(Yb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Yb(ip)),Yb(ip) if(INT(Zb(ip))==Ngrid+1)write(*,*)'Error in boundary: ',INT(Zb(ip)),Zb(ip) Xpar(ip+ioff) = Xb(ip) Ypar(ip+ioff) = Yb(ip) Zpar(ip+ioff) = Zb(ip) VX(ip+ioff) = VXb(ip) VY(ip+ioff) = VYb(ip) VZ(ip+ioff) = VZb(ip) end Do end DO close (20) xx = MAXVAL(Xpar) xm = MINVAL(Xpar) write(*,*)' x min/max= ',xm,xx xx = MAXVAL(Ypar) xm = MINVAL(Ypar) write(*,*)' y min/max= ',xm,xx xx = MAXVAL(Zpar) xm = MINVAL(Zpar) write(*,*)' z min/max= ',xm,xx Do ii =1,Nparticles If(Xpar(ii).lt.1.0.or.Ypar(ii).lt.1.0.or.Zpar(ii).lt.1.0)& write(*,'(a,i10,1p,3g14.5)')' Error coord: ',ii,Xpar(ii),Ypar(ii),Zpar(ii) EndDo DEALLOCATE (Xb,Yb,Zb,VXb,VYb,VZb) Allocate (FI(NGRID,NGRID,NGRID)) end SUBROUTINE ReadDataPM !--------------------------------------- ! Write PMfiles: PMcrd/crs SUBROUTINE WriteDataPM(iFlag) !--------------------------------------- use Structures Character*80 :: Name Logical :: exst Integer*8 :: iCount,ii,Jpage,Nrecord Integer*8 :: Ngal,moment,i,jfirst,jlast,j0 Call Timing(4,-1) Npage = 1024**2 ! number of particles per record Nrecord = Npage Jpage = Npage Naccess = Nrecord*6 !!*4 !!! Nrecpage = 512 ! max number of records per file moment = ISTEP Npages = (Nparticles-1)/JPAGE+1 Nfiles = (Npages-1)/Nrecpage +1 ! number of files write(*,*) ' Npages = ',Npages write(*,*) ' Write files = ',iFlag,moment ! open files If(iFlag == 0)Then Open (4,file ='PMcrd.DAT',form ='UNFORMATTED',status ='UNKNOWN') write(Name,'(a,i1,a)')'PMcrs',0,'.DAT' OPEN(UNIT=20,FILE=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) Else write(Name,'(a,i4.4,a)')'PMcrd.',moment,'.DAT' Open (4,file =TRIM(Name),form ='UNFORMATTED',status ='UNKNOWN') write(Name,'(a,i1,a,i4.4,a)')'PMcrs',0,'.',moment,'.DAT' write(*,'(a)') TRIM(Name) OPEN(UNIT=20,FILE=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) end If PARTW = float(NGRID)**3/FLOAT(Nparticles) AU0 = 0. write (4) HEADER, & AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, & TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0, & NROW,NGRID,Nspecs,Nseed,Om0,Oml0,hubble, & Nparticles,extras close(4) Allocate (Xb(Nrecord),Yb(Nrecord),Zb(Nrecord)) Allocate (VXb(Nrecord),VYb(Nrecord),VZb(Nrecord)) jj = 1 Do i=1,Npages !-------- dump into files If(i==Npages)Then NinPage = Nparticles -(i-1)*JPAGE ! # particles in the current page Else NinPage = JPAGE EndIf If(mod(i-1,Nrecpage)==0.and.i>1)Then ! close old and open new file close(20) jj = 1 ifile = (i-1)/Nrecpage ! construct file name If(iFlag==0)Then If(ifile<10)Then write(Name,'(a,i1.1,a)')'PMcrs',ifile,'.DAT' Else write(Name,'(a,i2.2,a)')'PMcrs',ifile,'.DAT' EndIf Else If(ifile<10)Then write(Name,'(a,i1.1,a,i4.4,a)')'PMcrs',ifile,'.',moment,'.DAT' Else write(Name,'(a,i2.2,a,i4.4,a)')'PMcrs',ifile,'.',moment,'.DAT' EndIf end If Open(20,file=TRIM(Name),ACCESS='DIRECT', & FORM='unformatted',STATUS='UNKNOWN',RECL=NACCESS) write(*,'(2i7,2a,3x,i9)') i,ifile,' Open file = ',TRIM(Name),Ninpage end If jfirst = (i-1)*JPAGE +1 ! first and last particles in current record jlast = jfirst + NinPage-1 If(mod(i,10)==0.or.i==Npages) & write(*,'(10x,a,i5,a,4i11)')'Write page=',i,' Particles=',NinPage,jfirst,jlast Do j0 = jfirst,jlast Xb(j0-jfirst+1) = Xpar(j0) Yb(j0-jfirst+1) = Ypar(j0) Zb(j0-jfirst+1) = Zpar(j0) Vxb(j0-jfirst+1) = VX(j0) VYb(j0-jfirst+1) = VY(j0) VZb(j0-jfirst+1) = VZ(j0) EndDo !write(*,'(i10,1p,6g13.5)') (k,XPAR(k),YPAR(k),ZPAR(k),VX(k),VY(k),VZ(k),k=1,1024) WRITE (20,REC=jj) Xb,Yb,Zb,VXb,VYb,VZb jj = jj +1 EndDo ! end lspecies loop DEALLOCATE (Xb,Yb,Zb,VXb,VYb,VZb) Call Timing(4,1) end SUBROUTINE WriteDataPM !-------------------------------------------------- SUBROUTINE SetTest !-------------------------------------------------- use Structures Integer*8 :: ii ISTEP =0 AEXPN = 0.01 ASTEP = 0.001 Om0 = 1. Oml = 0. write(*,*) ' Ngrid =',Ngrid ii = 0 DO M3=1,NGRID,2 DO M2=1,NGRID,2 DO M1=1,NGRID,2 ii = ii +1 If(ii>Nparticles)Stop ' Number of particles is to big in SetTest' XPAR(ii) = M1+0.05 YPAR(ii) = M2+0.05 ZPAR(ii) = M3+0.05 VX(ii) = 0. VY(ii) = 0. VZ(ii) = 0. END DO END DO END DO If(ii/= Nparticles)Stop ' Wrong number of particles in SetTest' end SUBROUTINE SetTest !-------------------------------------------------- ! Advance Aexpn, Istep, tIntg,... ! AEXPN = currnet expansion parameter ! ISTEP = current step ! ASTEP = step in the expansion parameter SUBROUTINE ADDTIME !-------------------------------------------------- use Structures ISTEP = ISTEP + 1 AEXPN = AEXPN + ASTEP ! Energy conservation IF(ISTEP.EQ.1)THEN EKIN1 = EKIN EKIN2 = 0. EKIN = ENKIN AU0 = AEXP0*ENPOT AEU0 = AEXP0*ENPOT + AEXP0*(EKIN+EKIN1)/2. TINTG = 0. WRITE (*,40) ISTEP,AEXPN,EKIN,ENPOT,AU0,AEU0 WRITE (17,40) ISTEP,AEXPN,EKIN,ENPOT,AU0,AEU0 40 FORMAT('**** STEP=',I3,' A=',F10.4,' E KIN=',E12.4, & ' E POT=',E12.4,/' AU0,AEU0=',2E12.4) ELSE EKIN2 = EKIN1 EKIN1 = EKIN EKIN = ENKIN TINTG = TINTG +ASTEP*(EKIN1 +(EKIN -2.*EKIN1 +EKIN2)/24.) ERROR = ((AEXPN-ASTEP)*((EKIN+EKIN1)/2.+ENPOT)-AEU0+TINTG)/ & ((AEXPN-ASTEP)*ENPOT)*100. WRITE (*,50) ISTEP,AEXPN,ERROR,EKIN,ENPOT,TINTG WRITE (17,50) ISTEP,AEXPN,ERROR,EKIN,ENPOT,TINTG ENDIF 50 FORMAT('Step = ',I4,' A=',F7.4,' Error(%)=',f7.2, & ' Ekin=',E11.3,' Epot=',E11.3,' Intg=',E11.3) END SUBROUTINE ADDTIME !------------------------------------------------------------ ! Advance each particle: dA AEXPN dA ! by one step I______._______I_______.______I -> A ! i-1 . i . i+1 step ! . {Fi} . . ! { vx } { x } . . ! { vy } { y } ^ . ! ._______________^ . ! . ^ ! .______________^ ! ! 0.5 ! dP = - A * Grad(Fi) * dA ; A =AEXPN ! i i ! 3/2 ! dX = P(new)/A * dA ; A =AEXPN+dA/2 ! i+1/2 i+1/2 ! !------------------------------------------------ SUBROUTINE MOVE(FactV) !------------------------------------------------ use structures ! PCONST = factor to change velocities ! XCONST = factor to change coordinates ! Note: 0.5 is needed in Pconst because ! Fi(i+1)-Fi(i-1) is used as gradient ! FactV = 1 - normal constant step ! = 1.2 - increase step by 1.5 real*8 :: SVEL,SPHI, PCONST, XCONST, XN, YN, & D1,D2,D3,T1,T2,T3, T2T1,T2D1,D2T1,D2D1, & GX,GY,GZ,FP,VVx,VVY,VVZ, & GX111,GX211,GX121,GX221, & GX112,GX212,GX122,GX222, & GY111,GY211,GY121,GY221, & GY112,GY212,GY122,GY222, & GZ111,GZ211,GZ121,GZ221, & GZ112,GZ212,GZ122,GZ222, & X,Y,Z integer*8 :: IN Call Timing(2,-1) PCONST = - SQRT(AEXPN/(Om0+Oml0*AEXPN**3))*ASTEP*0.5/FactV Ahalf = AEXPN+ASTEP/2. XCONST = ASTEP/SQRT(Ahalf*(Om0+Oml0*Ahalf**3))/Ahalf SVEL = 0. ! counter for \Sum(v_i**2) SPHI = 0. ! counter for \Sum(phi_i) XN = FLOAT(NGRID)+1.-1.E-8 ! N+1 YN = FLOAT(NGRID) ! N Wpar = YN**3/FLOAT(Nparticles) !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (X,Y,Z,GX,GY,GZ,FP,VVx,VVY,VVZ,I,J,K) & !$OMP PRIVATE (I1,J1,K1,K2,K3,I2,J2,I0,J0) & !$OMP PRIVATE (D1,D2,D3,T1,T2,T3, T2T1,T2D1,D2T1,D2D1) & !$OMP PRIVATE (F111,F211,F121,F221,F112,F212,F122,F222) & !$OMP PRIVATE (F113,F213,F123,F223,F110,F210,F120,F220) & !$OMP PRIVATE (F311,F321,F131,F231,F312,F322,F132,F232) & !$OMP PRIVATE (F011,F021,F101,F201,F012,F022,F102,F202) & !$OMP PRIVATE (GX111,GX211,GX121,GX221,GX112,GX212,GX122,GX222) & !$OMP PRIVATE (GY111,GY211,GY121,GY221,GY112,GY212,GY122,GY222) & !$OMP PRIVATE (GZ111,GZ211,GZ121,GZ221,GZ112,GZ212,GZ122,GZ222) & !$OMP REDUCTION(+:SVEL,SPHI) DO IN=1,Nparticles ! Loop over particles X=XPAR(IN) Y=YPAR(IN) Z=ZPAR(IN) VVX=VX(IN) VVY=VY(IN) VVZ=VZ(IN) I=INT(X) J=INT(Y) K=INT(Z) D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2T1 =T2*T1 T2D1 =T2*D1 D2T1 =D2*T1 D2D1 =D2*D1 I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 K2=K+2 IF(K2.GT.NGRID)K2=K2-NGRID K3=K-1 IF(K3.LT.1 )K3=NGRID F111 =FI(I ,J ,K ) ! Read potential to Fij vars F211 =FI(I1,J ,K ) F121 =FI(I ,J1,K ) F221 =FI(I1,J1,K ) F112 =FI(I ,J ,K1) F212 =FI(I1,J ,K1) F122 =FI(I ,J1,K1) F222 =FI(I1,J1,K1) F113 =FI(I ,J ,K2) F213 =FI(I1,J ,K2) F123 =FI(I ,J1,K2) F223 =FI(I1,J1,K2) F110 =FI(I ,J ,K3) F210 =FI(I1,J ,K3) F120 =FI(I ,J1,K3) F220 =FI(I1,J1,K3) I2=I+2 IF(I2.GT.NGRID)I2=I2-NGRID J2=J+2 IF(J2.GT.NGRID)J2=J2-NGRID F311 =FI(I2,J ,K ) F321 =FI(I2,J1,K ) F131 =FI(I ,J2,K ) F231 =FI(I1,J2,K ) F312 =FI(I2,J ,K1) F322 =FI(I2,J1,K1) F132 =FI(I ,J2,K1) F232 =FI(I1,J2,K1) I0=I-1 IF(I0.LT.1)I0=NGRID J0=J-1 IF(J0.LT.1)J0=NGRID F011 =FI(I0,J ,K ) F021 =FI(I0,J1,K ) F101 =FI(I ,J0,K ) F201 =FI(I1,J0,K ) F012 =FI(I0,J ,K1) F022 =FI(I0,J1,K1) F102 =FI(I ,J0,K1) F202 =FI(I1,J0,K1) ! Find {2*gradient} in nods GX111 =F211 -F011 GX211 =F311 -F111 GX121 =F221 -F021 GX221 =F321 -F121 GX112 =F212 -F012 GX212 =F312 -F112 GX122 =F222 -F022 GX222 =F322 -F122 GY111 =F121 -F101 GY211 =F221 -F201 GY121 =F131 -F111 GY221 =F231 -F211 GY112 =F122 -F102 GY212 =F222 -F202 GY122 =F132 -F112 GY222 =F232 -F212 GZ111 =F112 -F110 GZ211 =F212 -F210 GZ121 =F122 -F120 GZ221 =F222 -F220 GZ112 =F113 -F111 GZ212 =F213 -F211 GZ122 =F123 -F121 GZ222 =F223 -F221 ! Interpolate to the point GX=PCONST*(T3*(T2T1*GX111+T2D1*GX211 +D2T1*GX121+D2D1*GX221 )+ & D3*(T2T1*GX112+T2D1*GX212 +D2T1*GX122+D2D1*GX222 )) GY=PCONST*(T3*(T2T1*GY111+T2D1*GY211 +D2T1*GY121+D2D1*GY221 )+ & D3*(T2T1*GY112+T2D1*GY212 +D2T1*GY122+D2D1*GY222 )) GZ=PCONST*(T3*(T2T1*GZ111+T2D1*GZ211 +D2T1*GZ121+D2D1*GZ221 )+ & D3*(T2T1*GZ112+T2D1*GZ212 +D2T1*GZ122+D2D1*GZ222 )) ! Find potential of the point FP= T3*(T2T1*F111+T2D1*F211 +D2T1*F121+D2D1*F221 )+ & D3*(T2T1*F112+T2D1*F212 +D2T1*F122+D2D1*F222 ) SPHI = SPHI + FP*WPAR VVX =VVX+GX ! Move points VVY =VVY+GY VVZ =VVZ+GZ X =X +VVX*XCONST Y =Y +VVY*XCONST Z =Z +VVZ*XCONST IF(X.LT.1.d0)X=X+YN ! Periodical conditions IF(X.GE.XN)X=X-YN IF(Y.LT.1.d0)Y=Y+YN IF(Y.GE.XN)Y=Y-YN IF(Z.LT.1.d0)Z=Z+YN IF(Z.GE.XN)Z=Z-YN SVEL=SVEL+(VVX**2+VVY**2+VVZ**2)*WPAR XPAR(IN)=X ! Write new coordinates YPAR(IN)=Y ZPAR(IN)=Z VX(IN)=VVX VY(IN)=VVY VZ(IN)=VVZ if(INT(Xpar(IN))==Ngrid+1)Xpar(IN)=Xpar(IN)-1.e-3 if(INT(Ypar(IN))==Ngrid+1)Ypar(IN)=Ypar(IN)-1.e-3 if(INT(Zpar(IN))==Ngrid+1)Zpar(IN)=Zpar(IN)-1.e-3 ENDDO ! Set energies: ! Kin energy now at A(i+1/2) ! Pot energy at A(i) ENKIN = SVEL / 2. / (AEXPN+ASTEP/2.)**2 ENPOT = SPHI /2. Call Timing(2,1) !Do i =1,10000,100 ! write(*,'(a,i5,3f10.5,4x,1p,16G14.5)') ' Move : ',i,Xpar(i),Ypar(i),Zpar(i) !end Do end SUBROUTINE MOVE !------------------------------------------------- ! Find potential on Grid FI: DENSITY -> POTENTIAL ! ! O 1 ^ - Fourier component ! | ! 1 |-4 1 ^ ^ 2Pi ! O-----O-----O Fi = Rho / (2cos(--- (i-1))+ ! | i,j i,j Ngrid ! | ! O 1 2Pi ! 2cos(--- (j-1))-4) ! ^ Ngrid ! Fi = 1 (?) ! 11 ! 2 ! NABLA Fi = 3/2 /A * (Rho - ) ; ! X ! = 1 ! SUBROUTINE POTENTfft5 !--------------------------------------------- use Structures integer*4, parameter :: Nlensav = 8192 integer*4, parameter :: Nlenwrk = 8192 real*8, parameter :: P16 = 6.28318530718 real*4, parameter :: sq2 = 1. !1.41421356237 real*8, save :: wsave(1:Nlensav) real*8, save :: work(1:Nlenwrk) REAL*8 :: XX,D1,D2,A1,A2,A3,wi,wj,wk Integer*4 :: OMP_GET_MAX_THREADS,OMP_GET_THREAD_NUM Integer*4 :: Ng,ier,lensav,lenwrk,lenr,inc Real*8 :: GREENf(Nlenwrk) real*8 :: r(Nlenwrk) !$OMP THREADPRIVATE(work,wsave) Call Timing(1,-1) If(Ngrid>Nlenwrk)Stop ' Incresase Nlenwrk in POTENTfft5' Ng = Ngrid lensav = Ngrid+int(log(real(Ngrid,kind = 4))/log(2.0E+00))+4 lenwrk = Ngrid trfi = 1.5*Om0/aexpn write(*,*) ' Om0 =',Om0,trfi ! Set Green function components GREENf(1) = 2. GREENf(Ngrid) = -2. DO i=1,Ngrid/2-1 XX = 2.*COS(P16*i/Ngrid) GREENf(2*i) = XX GREENf(2*i+1) = XX End DO call rfft1i ( Ng, wsave, lensav, ier ) ! Initialize FFT inc = 1 lenr = Ngrid !write(*,'(a,3i8,3x,1p,10G15.7)')' PotentFFT5: ',Ngrid,lensav,lenwrk,(GREENf(i),i=1,10) write(*,*) ' time =',seconds(), ' XY fft' !$OMP PARALLEL DO DEFAULT(SHARED) copyin(wsave,work) & !$OMP PRIVATE ( k,j,i ,r,ier) Do k=1,NGRID ! fft for xy planes in x-dir Do j=1,NGRID Do i=1,NGRID r(i) = FI(i,j,k) EndDo call rfft1f ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do i=1,NGRID FI(i,j,k) = r(i) EndDo EndDo Do i=1,NGRID ! fft xy planes in y-dir Do j=1,NGRID r(j) = FI(i,j,k) EndDo call rfft1f ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do j=1,NGRID FI(i,j,k) = r(j) EndDo EndDo EndDo write(*,*) ' time =',seconds(), ' transposition' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,aa) DO J=1,Ngrid DO K=1,Ngrid-1 DO I=K+1,Ngrid aa = FI(I,J,K) FI(I,J,K) =FI(K,J,I) FI(K,J,I) =aa ENDDO ENDDO ENDDO write(*,*) ' time =',seconds(), ' Z forw/backw' !$OMP PARALLEL DO DEFAULT(SHARED) copyin(wsave,work) & !$OMP PRIVATE ( k,j,i ,r, ier,A1,A2,A3) Do j=1,NGRID ! ------ z-direction Do i=1,NGRID Do k=1,NGRID r(k) = FI(k,j,i) EndDo call rfft1f ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do k=1,NGRID FI(k,j,i) = r(k) EndDo EndDo end Do ww = (P16/Ngrid)**2 !$OMP PARALLEL DO DEFAULT(SHARED) copyin(wsave,work) & !$OMP PRIVATE ( k,j,i ,r, ier,A1,A2,A3,wi,wj,wk) Do j=1,NGRID ! ------ z-direction A3 = GREENf(J) -6. !wj = j/2 Do i=1,NGRID A2 = GREENf(I) + A3 !wi = i/2 Do k=1,NGRID A1 =A2 +GREENf(K) !--- use this for descrete Poisson solver ! wk = k/2 !A1 = -ww*(wi**2+wj**2+wk**2) !--- use this for k**2 Green functions !if(ww<5.)write(50,'(3i4,3x,3f7.3,f9.4,2g14.5)') i,j,k,wi,wj,wk,ww,A1,FI(k,j,i) IF(ABS(A1).LT.1.d-7) A1=1. r(k) = FI(k,j,i)*trfi/A1 EndDo call rfft1b ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do k=1,NGRID FI(k,j,i) = r(k) EndDo EndDo end Do write(*,*) ' time =',seconds(),' transpose ' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE ( k,j,i,aa) DO J=1,Ngrid DO K=1,Ngrid-1 DO I=K+1,Ngrid aa = FI(I,J,K) FI(I,J,K) =FI(K,J,I) FI(K,J,I) =aa ENDDO ENDDO ENDDO write(*,*) ' time =',seconds(), ' XY fft' !$OMP PARALLEL DO DEFAULT(SHARED) copyin(wsave,work) & !$OMP PRIVATE ( k,j,i ,r,ier) Do k=1,NGRID ! fft for xy planes in x-dir Do j=1,NGRID Do i=1,NGRID r(i) = FI(i,j,k) EndDo call rfft1b ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do i=1,NGRID FI(i,j,k) = r(i) EndDo EndDo Do i=1,NGRID ! fft xy planes in y-dir Do j=1,NGRID r(j) = FI(i,j,k) EndDo call rfft1b ( Ng, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) Do j=1,NGRID FI(i,j,k) = r(j) EndDo EndDo EndDo write(*,*) ' time =',seconds(), ' Finished Potent' Call Timing(1,1) end SUBROUTINE POTENTfft5 !------------------------------------------------------ SUBROUTINE DENSIT !------------------------------------------------------ use Structures real*8 :: XN,YN, & X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W integer*8 :: IN Call Timing(3,-1) XN =FLOAT(NGRID)+1.-1.E-7 YN =FLOAT(NGRID) Wpar = YN**3/FLOAT(Nparticles) ! Subtract mean density write(*,*) ' Init density' !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID FI(M1,M2,M3) = -1. END DO END DO END DO !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (IN,X,Y,Z,D1,D2,D3,T1,T2,T3,T2W,D2W) & !$OMP PRIVATE (I,J,K,I1,J1,K1) DO IN=1,Nparticles ! loop over particles X=XPAR(IN) Y=YPAR(IN) Z=ZPAR(IN) I=INT(X) J=INT(Y) K=INT(Z) D1=X-FLOAT(I) D2=Y-FLOAT(J) D3=Z-FLOAT(K) T1=1.-D1 T2=1.-D2 T3=1.-D3 T2W =T2*WPAR D2W =D2*WPAR I1=I+1 IF(I1.GT.NGRID)I1=1 J1=J+1 IF(J1.GT.NGRID)J1=1 K1=K+1 IF(K1.GT.NGRID)K1=1 !$OMP ATOMIC FI(I ,J ,K ) =FI(I ,J ,K ) +T3*T1*T2W !$OMP ATOMIC FI(I1,J ,K ) =FI(I1,J ,K ) +T3*D1*T2W !$OMP ATOMIC FI(I ,J1,K ) =FI(I ,J1,K ) +T3*T1*D2W !$OMP ATOMIC FI(I1,J1,K ) =FI(I1,J1,K ) +T3*D1*D2W !$OMP ATOMIC FI(I ,J ,K1) =FI(I ,J ,K1) +D3*T1*T2W !$OMP ATOMIC FI(I1,J ,K1) =FI(I1,J ,K1) +D3*D1*T2W !$OMP ATOMIC FI(I ,J1,K1) =FI(I ,J1,K1) +D3*T1*D2W !$OMP ATOMIC FI(I1,J1,K1) =FI(I1,J1,K1) +D3*D1*D2W ENDDO Call Timing(3,1) D1 = 0. ; D2 =0. !$OMP PARALLEL DO DEFAULT(SHARED) & !$OMP PRIVATE (M1,M2,M3) REDUCTION(+:D1,D2) DO M3=1,NGRID DO M2=1,NGRID DO M1=1,NGRID D1 = D1 + FI(M1,M2,M3) D2 = D2 + FI(M1,M2,M3)**2 END DO END DO END DO D1 = D1/(float(NGRID))**3 D2 = sqrt(D2/(float(NGRID))**3) write(*,*) ' Finished density: aver/rms=',D1,D2 END SUBROUTINE DENSIT ! ! ---------------------------------------------------- function seconds () ! ---------------------------------------------------- ! ! purpose: returns elapsed time in seconds Integer*8, SAVE :: first=0,rate=0,i0=0 Integer*8 :: i If(first==0)Then CALL SYSTEM_CLOCK(i,rate) first =1 i0 = i seconds = 0. Else CALL SYSTEM_CLOCK(i) seconds = float(i-i0)/float(rate) EndIf end function seconds !-------------------------------------------- subroutine Timing ( ielement , isign ) !-------------------------------------------- ! 0 - total ! 1 - force, 2- move ! 3 - density, 4- IO use Structures Real*8, SAVE :: CPU(0:5) =0. Character*80 :: FName='timing.log' Integer OMP_GET_THREAD_NUM, OMP_GET_MAX_THREADS If(isign == 0)Then Open(30,file=TRIM(FName),position='append') If(iStep==1)Then write(30,'(3(a,i11))') ' Ngrid = ',Ngrid, & ' Npart = ',Nparticles, & ' Nthreads = ',OMP_GET_MAX_THREADS() write(30,'(T3,a,T10,a,T23,a,T35,a,T49,a,T62,a)')& 'Step','Tot/min','Force','Move','Density','IO' End If write(30,'(i5,1p,10G13.4)'),iStep,CPU(0:4)/60. close(30) CPU(:) = 0 return EndIf CPU(ielement) = CPU(ielement) + float(isign) * seconds() end subroutine Timing !---------------------------------------------------------------------- ! FFT5 routines ! !---------------------------- subroutine rfft1b ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) !*****************************************************************************80 ! !! RFFT1B: real single precision backward fast Fourier transform, 1D. ! ! Discussion: ! ! RFFT1B computes the one-dimensional Fourier transform of a periodic ! sequence within a real array. This is referred to as the backward ! transform or Fourier synthesis, transforming the sequence from ! spectral to physical space. This transform is normalized since a ! call to RFFT1B followed by a call to RFFT1F (or vice-versa) reproduces ! the original array within roundoff error. ! ! License: ! ! Licensed under the GNU General Public License (GPL). ! Copyright (C) 1995-2004, Scientific Computing Division, ! University Corporation for Atmospheric Research ! ! Modified: ! ! 25 March 2005 ! ! Author: ! ! Paul Swarztrauber ! Richard Valent ! ! Reference: ! ! Paul Swarztrauber, ! Vectorizing the Fast Fourier Transforms, ! in Parallel Computations, ! edited by G. Rodrigue, ! Academic Press, 1982. ! ! Paul Swarztrauber, ! Fast Fourier Transform Algorithms for Vector Computers, ! Parallel Computing, pages 45-63, 1984. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the length of the sequence to be ! transformed. The transform is most efficient when N is a product of ! small primes. ! ! Input, integer ( kind = 4 ) INC, the increment between the locations, ! in array R, of two consecutive elements within the sequence. ! ! Input/output, real ( kind = 4 ) R(LENR), on input, the data to be ! transformed, and on output, the transformed data. ! ! Input, integer ( kind = 4 ) LENR, the dimension of the R array. ! LENR must be at least INC*(N-1) + 1. ! ! Input, real ( kind = 4 ) WSAVE(LENSAV). WSAVE's contents must be ! initialized with a call to RFFT1I before the first call to routine ! RFFT1F or RFFT1B for a given transform length N. ! ! Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array. ! LENSAV must be at least N + INT(LOG(REAL(N))) + 4. ! ! Workspace, real ( kind = 4 ) WORK(LENWRK). ! ! Input, integer ( kind = 4 ) LENWRK, the dimension of the WORK array. ! LENWRK must be at least N. ! ! Output, integer ( kind = 4 ) IER, error flag. ! 0, successful exit; ! 1, input parameter LENR not big enough; ! 2, input parameter LENSAV not big enough; ! 3, input parameter LENWRK not big enough. ! implicit none integer ( kind = 4 ) lenr integer ( kind = 4 ) lensav integer ( kind = 4 ) lenwrk integer ( kind = 4 ) ier integer ( kind = 4 ) inc integer ( kind = 4 ) n real ( kind = 8 ) r(lenr) real ( kind = 8 ) work(lenwrk) real ( kind = 8 ) wsave(lensav) ier = 0 if ( lenr < inc * ( n - 1 ) + 1 ) then ier = 1 call xerfft ( 'rfft1b ', 6 ) return end if if ( lensav < n + int ( log ( real ( n, kind = 4 ) ) ) + 4 ) then ier = 2 call xerfft ( 'rfft1b ', 8 ) return end if if ( lenwrk < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RFFT1B - Fatal error!' write ( *, '(a)' ) ' LENWRK < N:' write ( *, '(a,i6)' ) ' LENWRK = ', lenwrk write ( *, '(a,i6)' ) ' N = ', n ier = 3 call xerfft ( 'rfft1b ', 10 ) return end if if ( n == 1 ) then return end if call rfftb1 ( n, inc, r, work, wsave, wsave(n+1) ) return end subroutine rfft1b subroutine rfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier ) !*****************************************************************************80 ! !! RFFT1F: real single precision forward fast Fourier transform, 1D. ! ! Discussion: ! ! RFFT1F computes the one-dimensional Fourier transform of a periodic ! sequence within a real array. This is referred to as the forward ! transform or Fourier analysis, transforming the sequence from physical ! to spectral space. This transform is normalized since a call to ! RFFT1F followed by a call to RFFT1B (or vice-versa) reproduces the ! original array within roundoff error. ! ! License: ! ! Licensed under the GNU General Public License (GPL). ! Copyright (C) 1995-2004, Scientific Computing Division, ! University Corporation for Atmospheric Research ! ! Modified: ! ! 25 March 2005 ! ! Author: ! ! Paul Swarztrauber ! Richard Valent ! ! Reference: ! ! Paul Swarztrauber, ! Vectorizing the Fast Fourier Transforms, ! in Parallel Computations, ! edited by G. Rodrigue, ! Academic Press, 1982. ! ! Paul Swarztrauber, ! Fast Fourier Transform Algorithms for Vector Computers, ! Parallel Computing, pages 45-63, 1984. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the length of the sequence to be ! transformed. The transform is most efficient when N is a product of ! small primes. ! ! Input, integer ( kind = 4 ) INC, the increment between the locations, in ! array R, of two consecutive elements within the sequence. ! ! Input/output, real ( kind = 8 ) R(LENR), on input, contains the sequence ! to be transformed, and on output, the transformed data. ! ! Input, integer ( kind = 4 ) LENR, the dimension of the R array. ! LENR must be at least INC*(N-1) + 1. ! ! Input, real ( kind = 8 ) WSAVE(LENSAV). WSAVE's contents must be ! initialized with a call to RFFT1I before the first call to routine RFFT1F ! or RFFT1B for a given transform length N. ! ! Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array. ! LENSAV must be at least N + INT(LOG(REAL(N))) + 4. ! ! Workspace, real ( kind = 8 ) WORK(LENWRK). ! ! Input, integer ( kind = 4 ) LENWRK, the dimension of the WORK array. ! LENWRK must be at least N. ! ! Output, integer ( kind = 4 ) IER, error flag. ! 0, successful exit; ! 1, input parameter LENR not big enough: ! 2, input parameter LENSAV not big enough; ! 3, input parameter LENWRK not big enough. ! implicit none integer ( kind = 4 ) lenr integer ( kind = 4 ) lensav integer ( kind = 4 ) lenwrk integer ( kind = 4 ) ier integer ( kind = 4 ) inc integer ( kind = 4 ) n real ( kind = 8 ) work(lenwrk) real ( kind = 8 ) wsave(lensav) real ( kind = 8 ) r(lenr) ier = 0 if ( lenr < inc * ( n - 1 ) + 1 ) then ier = 1 call xerfft ( 'rfft1f ', 6 ) return end if if ( lensav < n + int ( log ( real ( n, kind = 4 ) ) ) + 4 ) then ier = 2 call xerfft ( 'rfft1f ', 8 ) return end if if ( lenwrk < n ) then ier = 3 call xerfft ( 'rfft1f ', 10 ) return end if if ( n == 1 ) then return end if call rfftf1 ( n, inc, r, work, wsave, wsave(n+1) ) return end subroutine rfft1f subroutine rfft1i ( n, wsave, lensav, ier ) !*****************************************************************************80 ! !! RFFT1I: initialization for RFFT1B and RFFT1F. ! ! Discussion: ! ! RFFT1I initializes array WSAVE for use in its companion routines ! RFFT1B and RFFT1F. The prime factorization of N together with a ! tabulation of the trigonometric functions are computed and stored ! in array WSAVE. Separate WSAVE arrays are required for different ! values of N. ! ! License: ! ! Licensed under the GNU General Public License (GPL). ! Copyright (C) 1995-2004, Scientific Computing Division, ! University Corporation for Atmospheric Research ! ! Modified: ! ! 25 March 2005 ! ! Author: ! ! Paul Swarztrauber ! Richard Valent ! ! Reference: ! ! Paul Swarztrauber, ! Vectorizing the Fast Fourier Transforms, ! in Parallel Computations, ! edited by G. Rodrigue, ! Academic Press, 1982. ! ! Paul Swarztrauber, ! Fast Fourier Transform Algorithms for Vector Computers, ! Parallel Computing, pages 45-63, 1984. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the length of the sequence to be ! transformed. The transform is most efficient when N is a product of ! small primes. ! ! Output, real ( kind = 8 ) WSAVE(LENSAV), containing the prime factors of ! N and also containing certain trigonometric values which will be used in ! routines RFFT1B or RFFT1F. ! ! Input, integer ( kind = 4 ) LENSAV, the dimension of the WSAVE array. ! LENSAV must be at least N + INT(LOG(REAL(N))) + 4. ! ! Output, integer ( kind = 4 ) IER, error flag. ! 0, successful exit; ! 2, input parameter LENSAV not big enough. ! implicit none integer ( kind = 4 ) lensav integer ( kind = 4 ) ier integer ( kind = 4 ) n real ( kind = 8 ) wsave(lensav) ier = 0 if ( lensav < n + int ( log ( real ( n, kind = 4 ) ) ) + 4 ) then ier = 2 call xerfft ( 'rfft1i ', 3 ) return end if if ( n == 1 ) then return end if call rffti1 ( n, wsave(1), wsave(n+1) ) return end subroutine rfft1i subroutine xerfft ( srname, info ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) info character ( len = * ) srname write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'XERFFT - Fatal error!' if ( 1 <= info ) then write ( *, '(a,a,a,i3,a)') ' On entry to ', trim ( srname ), & ' parameter number ', info, ' had an illegal value.' else if ( info == -1 ) then write( *, '(a,a,a,a)') ' On entry to ', trim ( srname ), & ' parameters LOT, JUMP, N and INC are inconsistent.' else if ( info == -2 ) then write( *, '(a,a,a,a)') ' On entry to ', trim ( srname ), & ' parameter L is greater than LDIM.' else if ( info == -3 ) then write( *, '(a,a,a,a)') ' On entry to ', trim ( srname ), & ' parameter M is greater than MDIM.' else if ( info == -5 ) then write( *, '(a,a,a,a)') ' Within ', trim ( srname ), & ' input error returned by lower level routine.' else if ( info == -6 ) then write( *, '(a,a,a,a)') ' On entry to ', trim ( srname ), & ' parameter LDIM is less than 2*(L/2+1).' end if stop end subroutine xerfft subroutine rfftb1 ( n, in, c, ch, wa, fac ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) in integer ( kind = 4 ) n real ( kind = 8 ) c(in,*) real ( kind = 8 ) ch(*) real ( kind = 8 ) fac(15) real ( kind = 8 ) half real ( kind = 8 ) halfm integer ( kind = 4 ) idl1 integer ( kind = 4 ) ido integer ( kind = 4 ) ip integer ( kind = 4 ) iw integer ( kind = 4 ) ix2 integer ( kind = 4 ) ix3 integer ( kind = 4 ) ix4 integer ( kind = 4 ) j integer ( kind = 4 ) k1 integer ( kind = 4 ) l1 integer ( kind = 4 ) l2 integer ( kind = 4 ) modn integer ( kind = 4 ) na integer ( kind = 4 ) nf integer ( kind = 4 ) nl real ( kind = 8 ) wa(n) nf = int ( fac(2) ) na = 0 do k1 = 1, nf ip = int ( fac(k1+2) ) na = 1 - na if ( 5 < ip ) then if ( k1 /= nf ) then na = 1 - na end if end if end do half = 0.5E+00 halfm = -0.5E+00 modn = mod ( n, 2 ) nl = n - 2 if ( modn /= 0 ) then nl = n - 1 end if if ( na == 0 ) then do j = 2, nl, 2 c(1,j) = half * c(1,j) c(1,j+1) = halfm * c(1,j+1) end do else ch(1) = c(1,1) ch(n) = c(1,n) do j = 2, nl, 2 ch(j) = half*c(1,j) ch(j+1) = halfm*c(1,j+1) end do end if l1 = 1 iw = 1 do k1 = 1, nf ip = int ( fac(k1+2) ) l2 = ip * l1 ido = n / l2 idl1 = ido * l1 if ( ip == 4 ) then ix2 = iw + ido ix3 = ix2 + ido if ( na == 0 ) then call r1f4kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3) ) else call r1f4kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3) ) end if na = 1 - na else if ( ip == 2 ) then if ( na == 0 ) then call r1f2kb ( ido, l1, c, in, ch, 1, wa(iw) ) else call r1f2kb ( ido, l1, ch, 1, c, in, wa(iw) ) end if na = 1 - na else if ( ip == 3 ) then ix2 = iw + ido if ( na == 0 ) then call r1f3kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2) ) else call r1f3kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2) ) end if na = 1 - na else if ( ip == 5 ) then ix2 = iw + ido ix3 = ix2 + ido ix4 = ix3 + ido if ( na == 0 ) then call r1f5kb ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3), wa(ix4) ) else call r1f5kb ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3), wa(ix4) ) end if na = 1 - na else if ( na == 0 ) then call r1fgkb ( ido, ip, l1, idl1, c, c, c, in, ch, ch, 1, wa(iw) ) else call r1fgkb ( ido, ip, l1, idl1, ch, ch, ch, 1, c, c, in, wa(iw) ) end if if ( ido == 1 ) then na = 1 - na end if end if l1 = l2 iw = iw + ( ip - 1 ) * ido end do return end subroutine rfftb1 ! !--------------------------------------------- subroutine r1fgkb ( ido, ip, l1, idl1, cc, c1, c2, in1, ch, ch2, in2, wa ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) idl1 integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) ip integer ( kind = 4 ) l1 real ( kind = 8 ) ai1 real ( kind = 8 ) ai2 real ( kind = 8 ) ar1 real ( kind = 8 ) ar1h real ( kind = 8 ) ar2 real ( kind = 8 ) ar2h real ( kind = 8 ) arg real ( kind = 8 ) c1(in1,ido,l1,ip) real ( kind = 8 ) c2(in1,idl1,ip) real ( kind = 8 ) cc(in1,ido,ip,l1) real ( kind = 8 ) ch(in2,ido,l1,ip) real ( kind = 8 ) ch2(in2,idl1,ip) real ( kind = 8 ) dc2 real ( kind = 8 ) dcp real ( kind = 8 ) ds2 real ( kind = 8 ) dsp integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idij integer ( kind = 4 ) idp2 integer ( kind = 4 ) ik integer ( kind = 4 ) ipp2 integer ( kind = 4 ) ipph integer ( kind = 4 ) is integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) jc integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) lc integer ( kind = 4 ) nbd real ( kind = 8 ) tpi real ( kind = 8 ) wa(ido) tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) arg = tpi / real ( ip, kind = 4 ) dcp = cos ( arg ) dsp = sin ( arg ) idp2 = ido + 2 nbd = ( ido - 1 ) / 2 ipp2 = ip + 2 ipph = ( ip + 1 ) / 2 if ( ido < l1 ) then do i = 1, ido do k = 1, l1 ch(1,i,k,1) = cc(1,i,1,k) end do end do else do k = 1, l1 do i = 1, ido ch(1,i,k,1) = cc(1,i,1,k) end do end do end if do j = 2, ipph jc = ipp2 - j j2 = j + j do k = 1, l1 ch(1,1,k,j) = cc(1,ido,j2-2,k)+cc(1,ido,j2-2,k) ch(1,1,k,jc) = cc(1,1,j2-1,k)+cc(1,1,j2-1,k) end do end do if ( ido == 1 ) then else if ( nbd < l1 ) then do j = 2, ipph jc = ipp2 - j do i = 3, ido, 2 ic = idp2 - i do k = 1, l1 ch(1,i-1,k,j) = cc(1,i-1,2*j-1,k)+cc(1,ic-1,2*j-2,k) ch(1,i-1,k,jc) = cc(1,i-1,2*j-1,k)-cc(1,ic-1,2*j-2,k) ch(1,i,k,j) = cc(1,i,2*j-1,k)-cc(1,ic,2*j-2,k) ch(1,i,k,jc) = cc(1,i,2*j-1,k)+cc(1,ic,2*j-2,k) end do end do end do else do j = 2, ipph jc = ipp2 - j do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,k,j) = cc(1,i-1,2*j-1,k)+cc(1,ic-1,2*j-2,k) ch(1,i-1,k,jc) = cc(1,i-1,2*j-1,k)-cc(1,ic-1,2*j-2,k) ch(1,i,k,j) = cc(1,i,2*j-1,k)-cc(1,ic,2*j-2,k) ch(1,i,k,jc) = cc(1,i,2*j-1,k)+cc(1,ic,2*j-2,k) end do end do end do end if ar1 = 1.0E+00 ai1 = 0.0E+00 do l = 2, ipph lc = ipp2 - l ar1h = dcp * ar1 - dsp * ai1 ai1 = dcp * ai1 + dsp * ar1 ar1 = ar1h do ik = 1, idl1 c2(1,ik,l) = ch2(1,ik,1)+ar1*ch2(1,ik,2) c2(1,ik,lc) = ai1*ch2(1,ik,ip) end do dc2 = ar1 ds2 = ai1 ar2 = ar1 ai2 = ai1 do j = 3, ipph jc = ipp2 - j ar2h = dc2*ar2-ds2*ai2 ai2 = dc2*ai2+ds2*ar2 ar2 = ar2h do ik = 1, idl1 c2(1,ik,l) = c2(1,ik,l)+ar2*ch2(1,ik,j) c2(1,ik,lc) = c2(1,ik,lc)+ai2*ch2(1,ik,jc) end do end do end do do j = 2, ipph do ik = 1, idl1 ch2(1,ik,1) = ch2(1,ik,1)+ch2(1,ik,j) end do end do do j = 2, ipph jc = ipp2 - j do k = 1, l1 ch(1,1,k,j) = c1(1,1,k,j)-c1(1,1,k,jc) ch(1,1,k,jc) = c1(1,1,k,j)+c1(1,1,k,jc) end do end do if ( ido == 1 ) then else if ( nbd < l1 ) then do j = 2, ipph jc = ipp2 - j do i = 3, ido, 2 do k = 1, l1 ch(1,i-1,k,j) = c1(1,i-1,k,j) - c1(1,i,k,jc) ch(1,i-1,k,jc) = c1(1,i-1,k,j) + c1(1,i,k,jc) ch(1,i,k,j) = c1(1,i,k,j) + c1(1,i-1,k,jc) ch(1,i,k,jc) = c1(1,i,k,j) - c1(1,i-1,k,jc) end do end do end do else do j = 2, ipph jc = ipp2 - j do k = 1, l1 do i = 3, ido, 2 ch(1,i-1,k,j) = c1(1,i-1,k,j)-c1(1,i,k,jc) ch(1,i-1,k,jc) = c1(1,i-1,k,j)+c1(1,i,k,jc) ch(1,i,k,j) = c1(1,i,k,j)+c1(1,i-1,k,jc) ch(1,i,k,jc) = c1(1,i,k,j)-c1(1,i-1,k,jc) end do end do end do end if if ( ido == 1 ) then return end if do ik = 1, idl1 c2(1,ik,1) = ch2(1,ik,1) end do do j = 2, ip do k = 1, l1 c1(1,1,k,j) = ch(1,1,k,j) end do end do if ( l1 < nbd ) then is = -ido do j = 2, ip is = is + ido do k = 1, l1 idij = is do i = 3, ido, 2 idij = idij + 2 c1(1,i-1,k,j) = wa(idij-1)*ch(1,i-1,k,j)-wa(idij)* ch(1,i,k,j) c1(1,i,k,j) = wa(idij-1)*ch(1,i,k,j)+wa(idij)* ch(1,i-1,k,j) end do end do end do else is = -ido do j = 2, ip is = is + ido idij = is do i = 3, ido, 2 idij = idij + 2 do k = 1, l1 c1(1,i-1,k,j) = wa(idij-1) * ch(1,i-1,k,j) - wa(idij) * ch(1,i,k,j) c1(1,i,k,j) = wa(idij-1) * ch(1,i,k,j) + wa(idij) * ch(1,i-1,k,j) end do end do end do end if return end subroutine r1fgkb ! !--------------------------------------------- subroutine r1fgkf ( ido, ip, l1, idl1, cc, c1, c2, in1, ch, ch2, in2, wa ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) idl1 integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) ip integer ( kind = 4 ) l1 real ( kind = 8 ) ai1 real ( kind = 8 ) ai2 real ( kind = 8 ) ar1 real ( kind = 8 ) ar1h real ( kind = 8 ) ar2 real ( kind = 8 ) ar2h real ( kind = 8 ) arg real ( kind = 8 ) c1(in1,ido,l1,ip) real ( kind = 8 ) c2(in1,idl1,ip) real ( kind = 8 ) cc(in1,ido,ip,l1) real ( kind = 8 ) ch(in2,ido,l1,ip) real ( kind = 8 ) ch2(in2,idl1,ip) real ( kind = 8 ) dc2 real ( kind = 8 ) dcp real ( kind = 8 ) ds2 real ( kind = 8 ) dsp integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idij integer ( kind = 4 ) idp2 integer ( kind = 4 ) ik integer ( kind = 4 ) ipp2 integer ( kind = 4 ) ipph integer ( kind = 4 ) is integer ( kind = 4 ) j integer ( kind = 4 ) j2 integer ( kind = 4 ) jc integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) lc integer ( kind = 4 ) nbd real ( kind = 8 ) tpi real ( kind = 8 ) wa(ido) tpi = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) arg = tpi / real ( ip, kind = 4 ) dcp = cos ( arg ) dsp = sin ( arg ) ipph = ( ip + 1 ) / 2 ipp2 = ip + 2 idp2 = ido + 2 nbd = ( ido - 1 ) / 2 if ( ido == 1 ) then do ik = 1, idl1 c2(1,ik,1) = ch2(1,ik,1) end do else do ik = 1, idl1 ch2(1,ik,1) = c2(1,ik,1) end do do j = 2, ip do k = 1, l1 ch(1,1,k,j) = c1(1,1,k,j) end do end do if ( l1 < nbd ) then is = -ido do j = 2, ip is = is + ido do k = 1, l1 idij = is do i = 3, ido, 2 idij = idij + 2 ch(1,i-1,k,j) = wa(idij-1)*c1(1,i-1,k,j)+wa(idij) *c1(1,i,k,j) ch(1,i,k,j) = wa(idij-1)*c1(1,i,k,j)-wa(idij) *c1(1,i-1,k,j) end do end do end do else is = -ido do j = 2, ip is = is + ido idij = is do i = 3, ido, 2 idij = idij + 2 do k = 1, l1 ch(1,i-1,k,j) = wa(idij-1)*c1(1,i-1,k,j)+wa(idij) *c1(1,i,k,j) ch(1,i,k,j) = wa(idij-1)*c1(1,i,k,j)-wa(idij) *c1(1,i-1,k,j) end do end do end do end if if ( nbd < l1 ) then do j = 2, ipph jc = ipp2 - j do i = 3, ido, 2 do k = 1, l1 c1(1,i-1,k,j) = ch(1,i-1,k,j)+ch(1,i-1,k,jc) c1(1,i-1,k,jc) = ch(1,i,k,j)-ch(1,i,k,jc) c1(1,i,k,j) = ch(1,i,k,j)+ch(1,i,k,jc) c1(1,i,k,jc) = ch(1,i-1,k,jc)-ch(1,i-1,k,j) end do end do end do else do j = 2, ipph jc = ipp2 - j do k = 1, l1 do i = 3, ido, 2 c1(1,i-1,k,j) = ch(1,i-1,k,j)+ch(1,i-1,k,jc) c1(1,i-1,k,jc) = ch(1,i,k,j)-ch(1,i,k,jc) c1(1,i,k,j) = ch(1,i,k,j)+ch(1,i,k,jc) c1(1,i,k,jc) = ch(1,i-1,k,jc)-ch(1,i-1,k,j) end do end do end do end if end if do j = 2, ipph jc = ipp2 - j do k = 1, l1 c1(1,1,k,j) = ch(1,1,k,j)+ch(1,1,k,jc) c1(1,1,k,jc) = ch(1,1,k,jc)-ch(1,1,k,j) end do end do ar1 = 1.0E+00 ai1 = 0.0E+00 do l = 2, ipph lc = ipp2 - l ar1h = dcp * ar1 - dsp * ai1 ai1 = dcp * ai1 + dsp * ar1 ar1 = ar1h do ik = 1, idl1 ch2(1,ik,l) = c2(1,ik,1)+ar1*c2(1,ik,2) ch2(1,ik,lc) = ai1*c2(1,ik,ip) end do dc2 = ar1 ds2 = ai1 ar2 = ar1 ai2 = ai1 do j = 3, ipph jc = ipp2 - j ar2h = dc2 * ar2 - ds2 * ai2 ai2 = dc2 * ai2 + ds2 * ar2 ar2 = ar2h do ik = 1, idl1 ch2(1,ik,l) = ch2(1,ik,l)+ar2*c2(1,ik,j) ch2(1,ik,lc) = ch2(1,ik,lc)+ai2*c2(1,ik,jc) end do end do end do do j = 2, ipph do ik = 1, idl1 ch2(1,ik,1) = ch2(1,ik,1)+c2(1,ik,j) end do end do if ( ido < l1 ) then do i = 1, ido do k = 1, l1 cc(1,i,1,k) = ch(1,i,k,1) end do end do else do k = 1, l1 do i = 1, ido cc(1,i,1,k) = ch(1,i,k,1) end do end do end if do j = 2, ipph jc = ipp2 - j j2 = j+j do k = 1, l1 cc(1,ido,j2-2,k) = ch(1,1,k,j) cc(1,1,j2-1,k) = ch(1,1,k,jc) end do end do if ( ido == 1 ) then return end if if ( nbd < l1 ) then do j = 2, ipph jc = ipp2 - j j2 = j + j do i = 3, ido, 2 ic = idp2 - i do k = 1, l1 cc(1,i-1,j2-1,k) = ch(1,i-1,k,j)+ch(1,i-1,k,jc) cc(1,ic-1,j2-2,k) = ch(1,i-1,k,j)-ch(1,i-1,k,jc) cc(1,i,j2-1,k) = ch(1,i,k,j)+ch(1,i,k,jc) cc(1,ic,j2-2,k) = ch(1,i,k,jc)-ch(1,i,k,j) end do end do end do else do j = 2, ipph jc = ipp2 - j j2 = j + j do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i cc(1,i-1,j2-1,k) = ch(1,i-1,k,j) + ch(1,i-1,k,jc) cc(1,ic-1,j2-2,k) = ch(1,i-1,k,j) - ch(1,i-1,k,jc) cc(1,i,j2-1,k) = ch(1,i,k,j) + ch(1,i,k,jc) cc(1,ic,j2-2,k) = ch(1,i,k,jc) - ch(1,i,k,j) end do end do end do end if return end subroutine r1fgkf ! !------------------------------------------------------- subroutine r1f2kb ( ido, l1, cc, in1, ch, in2, wa1 ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) cc(in1,ido,2,l1) real ( kind = 8 ) ch(in2,ido,l1,2) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) wa1(ido) do k = 1, l1 ch(1,1,k,1) = cc(1,1,1,k) + cc(1,ido,2,k) ch(1,1,k,2) = cc(1,1,1,k) - cc(1,ido,2,k) end do if ( ido < 2 ) then return end if if ( 2 < ido ) then idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,k,1) = cc(1,i-1,1,k) + cc(1,ic-1,2,k) ch(1,i,k,1) = cc(1,i,1,k) - cc(1,ic,2,k) ch(1,i-1,k,2) = wa1(i-2) * ( cc(1,i-1,1,k) - cc(1,ic-1,2,k) ) & - wa1(i-1) * ( cc(1,i,1,k) + cc(1,ic,2,k) ) ch(1,i,k,2) = wa1(i-2) * ( cc(1,i,1,k) + cc(1,ic,2,k) ) & + wa1(i-1) * ( cc(1,i-1,1,k) - cc(1,ic-1,2,k) ) end do end do if ( mod ( ido, 2 ) == 1 ) then return end if end if do k = 1, l1 ch(1,ido,k,1) = cc(1,ido,1,k) + cc(1,ido,1,k) ch(1,ido,k,2) = - ( cc(1,1,2,k) + cc(1,1,2,k) ) end do return end subroutine r1f2kb ! !------------------------------------------------------- subroutine r1f2kf ( ido, l1, cc, in1, ch, in2, wa1 ) !*****************************************************************************80 implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) ch(in2,ido,2,l1) real ( kind = 8 ) cc(in1,ido,l1,2) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) wa1(ido) do k = 1, l1 ch(1,1,1,k) = cc(1,1,k,1) + cc(1,1,k,2) ch(1,ido,2,k) = cc(1,1,k,1) - cc(1,1,k,2) end do if ( ido < 2 ) then return end if if ( 2 < ido ) then idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i,1,k) = cc(1,i,k,1) + ( wa1(i-2) * cc(1,i,k,2) & - wa1(i-1) * cc(1,i-1,k,2) ) ch(1,ic,2,k) = -cc(1,i,k,1) + ( wa1(i-2) * cc(1,i,k,2) & - wa1(i-1) * cc(1,i-1,k,2) ) ch(1,i-1,1,k) = cc(1,i-1,k,1) + ( wa1(i-2) * cc(1,i-1,k,2) & + wa1(i-1) * cc(1,i,k,2)) ch(1,ic-1,2,k) = cc(1,i-1,k,1) - ( wa1(i-2) * cc(1,i-1,k,2) & + wa1(i-1) * cc(1,i,k,2)) end do end do if ( mod ( ido, 2 ) == 1 ) then return end if end if do k = 1, l1 ch(1,1,2,k) = -cc(1,ido,k,2) ch(1,ido,1,k) = cc(1,ido,k,1) end do return end subroutine r1f2kf ! !------------------------------------------------------- subroutine r1f3kb ( ido, l1, cc, in1, ch, in2, wa1, wa2 ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) arg real ( kind = 8 ) cc(in1,ido,3,l1) real ( kind = 8 ) ch(in2,ido,l1,3) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) taui real ( kind = 8 ) taur real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 3.0E+00 taur = cos ( arg ) taui = sin ( arg ) do k = 1, l1 ch(1,1,k,1) = cc(1,1,1,k) + 2.0E+00 * cc(1,ido,2,k) ch(1,1,k,2) = cc(1,1,1,k) + 2.0E+00 * taur * cc(1,ido,2,k) & - 2.0E+00 * taui * cc(1,1,3,k) ch(1,1,k,3) = cc(1,1,1,k) + 2.0E+00 * taur * cc(1,ido,2,k) & + 2.0E+00 * taui * cc(1,1,3,k) end do if ( ido == 1 ) then return end if idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,k,1) = cc(1,i-1,1,k)+(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) ch(1,i,k,1) = cc(1,i,1,k)+(cc(1,i,3,k)-cc(1,ic,2,k)) ch(1,i-1,k,2) = wa1(i-2)* & ((cc(1,i-1,1,k)+taur*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)))- & (taui*(cc(1,i,3,k)+cc(1,ic,2,k)))) -wa1(i-1)* & ((cc(1,i,1,k)+taur*(cc(1,i,3,k)-cc(1,ic,2,k)))+ & (taui*(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))) ch(1,i,k,2) = wa1(i-2)* & ((cc(1,i,1,k)+taur*(cc(1,i,3,k)-cc(1,ic,2,k)))+ & (taui*(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))) +wa1(i-1)* & ((cc(1,i-1,1,k)+taur*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)))- & (taui*(cc(1,i,3,k)+cc(1,ic,2,k)))) ch(1,i-1,k,3) = wa2(i-2)* & ((cc(1,i-1,1,k)+taur*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)))+ & (taui*(cc(1,i,3,k)+cc(1,ic,2,k)))) -wa2(i-1)* & ((cc(1,i,1,k)+taur*(cc(1,i,3,k)-cc(1,ic,2,k)))- & (taui*(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))) ch(1,i,k,3) = wa2(i-2)* & ((cc(1,i,1,k)+taur*(cc(1,i,3,k)-cc(1,ic,2,k)))- & (taui*(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))) +wa2(i-1)* & ((cc(1,i-1,1,k)+taur*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)))+ & (taui*(cc(1,i,3,k)+cc(1,ic,2,k)))) end do end do return end subroutine r1f3kb ! !------------------------------------------------------- subroutine r1f3kf ( ido, l1, cc, in1, ch, in2, wa1, wa2 ) !*****************************************************************************80 implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) arg real ( kind = 8 ) cc(in1,ido,l1,3) real ( kind = 8 ) ch(in2,ido,3,l1) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) taui real ( kind = 8 ) taur real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 3.0E+00 taur = cos ( arg ) taui = sin ( arg ) do k = 1, l1 ch(1,1,1,k) = cc(1,1,k,1) + ( cc(1,1,k,2) + cc(1,1,k,3) ) ch(1,1,3,k) = taui * ( cc(1,1,k,3) - cc(1,1,k,2) ) ch(1,ido,2,k) = cc(1,1,k,1) + taur * ( cc(1,1,k,2) + cc(1,1,k,3) ) end do if ( ido == 1 ) then return end if idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,1,k) = cc(1,i-1,k,1)+((wa1(i-2)*cc(1,i-1,k,2)+ & wa1(i-1)*cc(1,i,k,2))+(wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))) ch(1,i,1,k) = cc(1,i,k,1)+((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))) ch(1,i-1,3,k) = (cc(1,i-1,k,1)+taur*((wa1(i-2)* & cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2))+(wa2(i-2)* & cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3))))+(taui*((wa1(i-2)* & cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2))-(wa2(i-2)* & cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3)))) ch(1,ic-1,2,k) = (cc(1,i-1,k,1)+taur*((wa1(i-2)* & cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2))+(wa2(i-2)* & cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3))))-(taui*((wa1(i-2)* & cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2))-(wa2(i-2)* & cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3)))) ch(1,i,3,k) = (cc(1,i,k,1)+taur*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))))+(taui*((wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))) ch(1,ic,2,k) = (taui*((wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2))))-(cc(1,i,k,1)+taur*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3)))) end do end do return end subroutine r1f3kf ! !------------------------------------------------------- subroutine r1f4kb ( ido, l1, cc, in1, ch, in2, wa1, wa2, wa3 ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) cc(in1,ido,4,l1) real ( kind = 8 ) ch(in2,ido,l1,4) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) sqrt2 real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) real ( kind = 8 ) wa3(ido) sqrt2 = sqrt ( 2.0E+00 ) do k = 1, l1 ch(1,1,k,3) = ( cc(1,1,1,k) + cc(1,ido,4,k) ) & - ( cc(1,ido,2,k) + cc(1,ido,2,k) ) ch(1,1,k,1) = ( cc(1,1,1,k) + cc(1,ido,4,k) ) & + ( cc(1,ido,2,k) + cc(1,ido,2,k) ) ch(1,1,k,4) = ( cc(1,1,1,k) - cc(1,ido,4,k) ) & + ( cc(1,1,3,k) + cc(1,1,3,k) ) ch(1,1,k,2) = ( cc(1,1,1,k) - cc(1,ido,4,k) ) & - ( cc(1,1,3,k) + cc(1,1,3,k) ) end do if ( ido < 2 ) then return end if if ( 2 < ido ) then idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,k,1) = (cc(1,i-1,1,k)+cc(1,ic-1,4,k)) & +(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) ch(1,i,k,1) = (cc(1,i,1,k)-cc(1,ic,4,k)) & +(cc(1,i,3,k)-cc(1,ic,2,k)) ch(1,i-1,k,2) = wa1(i-2)*((cc(1,i-1,1,k)-cc(1,ic-1,4,k)) & -(cc(1,i,3,k)+cc(1,ic,2,k)))-wa1(i-1) & *((cc(1,i,1,k)+cc(1,ic,4,k))+(cc(1,i-1,3,k)-cc(1,ic-1,2,k))) ch(1,i,k,2) = wa1(i-2)*((cc(1,i,1,k)+cc(1,ic,4,k)) & +(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))+wa1(i-1) & *((cc(1,i-1,1,k)-cc(1,ic-1,4,k))-(cc(1,i,3,k)+cc(1,ic,2,k))) ch(1,i-1,k,3) = wa2(i-2)*((cc(1,i-1,1,k)+cc(1,ic-1,4,k)) & -(cc(1,i-1,3,k)+cc(1,ic-1,2,k)))-wa2(i-1) & *((cc(1,i,1,k)-cc(1,ic,4,k))-(cc(1,i,3,k)-cc(1,ic,2,k))) ch(1,i,k,3) = wa2(i-2)*((cc(1,i,1,k)-cc(1,ic,4,k)) & -(cc(1,i,3,k)-cc(1,ic,2,k)))+wa2(i-1) & *((cc(1,i-1,1,k)+cc(1,ic-1,4,k))-(cc(1,i-1,3,k) & +cc(1,ic-1,2,k))) ch(1,i-1,k,4) = wa3(i-2)*((cc(1,i-1,1,k)-cc(1,ic-1,4,k)) & +(cc(1,i,3,k)+cc(1,ic,2,k)))-wa3(i-1) & *((cc(1,i,1,k)+cc(1,ic,4,k))-(cc(1,i-1,3,k)-cc(1,ic-1,2,k))) ch(1,i,k,4) = wa3(i-2)*((cc(1,i,1,k)+cc(1,ic,4,k)) & -(cc(1,i-1,3,k)-cc(1,ic-1,2,k)))+wa3(i-1) & *((cc(1,i-1,1,k)-cc(1,ic-1,4,k))+(cc(1,i,3,k)+cc(1,ic,2,k))) end do end do if ( mod ( ido, 2 ) == 1 ) then return end if end if do k = 1, l1 ch(1,ido,k,1) = ( cc(1,ido,1,k) + cc(1,ido,3,k) ) & + ( cc(1,ido,1,k) + cc(1,ido,3,k)) ch(1,ido,k,2) = sqrt2 * ( ( cc(1,ido,1,k) - cc(1,ido,3,k) ) & - ( cc(1,1,2,k) + cc(1,1,4,k) ) ) ch(1,ido,k,3) = ( cc(1,1,4,k) - cc(1,1,2,k) ) & + ( cc(1,1,4,k) - cc(1,1,2,k) ) ch(1,ido,k,4) = -sqrt2 * ( ( cc(1,ido,1,k) - cc(1,ido,3,k) ) & + ( cc(1,1,2,k) + cc(1,1,4,k) ) ) end do return end subroutine r1f4kb ! !------------------------------------------------------- subroutine r1f4kf ( ido, l1, cc, in1, ch, in2, wa1, wa2, wa3 ) !*****************************************************************************80 ! Parameters: ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) cc(in1,ido,l1,4) real ( kind = 8 ) ch(in2,ido,4,l1) real ( kind = 8 ) hsqt2 integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) real ( kind = 8 ) wa3(ido) hsqt2 = sqrt ( 2.0E+00 ) / 2.0E+00 do k = 1, l1 ch(1,1,1,k) = ( cc(1,1,k,2) + cc(1,1,k,4) ) & + ( cc(1,1,k,1) + cc(1,1,k,3) ) ch(1,ido,4,k) = ( cc(1,1,k,1) + cc(1,1,k,3) ) & - ( cc(1,1,k,2) + cc(1,1,k,4) ) ch(1,ido,2,k) = cc(1,1,k,1) - cc(1,1,k,3) ch(1,1,3,k) = cc(1,1,k,4) - cc(1,1,k,2) end do if ( ido < 2 ) then return end if if ( 2 < ido ) then idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,1,k) = ((wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2))+(wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4)))+(cc(1,i-1,k,1)+(wa2(i-2)*cc(1,i-1,k,3)+ & wa2(i-1)*cc(1,i,k,3))) ch(1,ic-1,4,k) = (cc(1,i-1,k,1)+(wa2(i-2)*cc(1,i-1,k,3)+ & wa2(i-1)*cc(1,i,k,3)))-((wa1(i-2)*cc(1,i-1,k,2)+ & wa1(i-1)*cc(1,i,k,2))+(wa3(i-2)*cc(1,i-1,k,4)+ & wa3(i-1)*cc(1,i,k,4))) ch(1,i,1,k) = ((wa1(i-2)*cc(1,i,k,2)-wa1(i-1)* & cc(1,i-1,k,2))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4)))+(cc(1,i,k,1)+(wa2(i-2)*cc(1,i,k,3)- & wa2(i-1)*cc(1,i-1,k,3))) ch(1,ic,4,k) = ((wa1(i-2)*cc(1,i,k,2)-wa1(i-1)* & cc(1,i-1,k,2))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4)))-(cc(1,i,k,1)+(wa2(i-2)*cc(1,i,k,3)- & wa2(i-1)*cc(1,i-1,k,3))) ch(1,i-1,3,k) = ((wa1(i-2)*cc(1,i,k,2)-wa1(i-1)* & cc(1,i-1,k,2))-(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4)))+(cc(1,i-1,k,1)-(wa2(i-2)*cc(1,i-1,k,3)+ & wa2(i-1)*cc(1,i,k,3))) ch(1,ic-1,2,k) = (cc(1,i-1,k,1)-(wa2(i-2)*cc(1,i-1,k,3)+ & wa2(i-1)*cc(1,i,k,3)))-((wa1(i-2)*cc(1,i,k,2)-wa1(i-1)* & cc(1,i-1,k,2))-(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4))) ch(1,i,3,k) = ((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))+(cc(1,i,k,1)-(wa2(i-2)*cc(1,i,k,3)- & wa2(i-1)*cc(1,i-1,k,3))) ch(1,ic,2,k) = ((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))-(cc(1,i,k,1)-(wa2(i-2)*cc(1,i,k,3)- & wa2(i-1)*cc(1,i-1,k,3))) end do end do if ( mod ( ido, 2 ) == 1 ) then return end if end if do k = 1, l1 ch(1,ido,1,k) = (hsqt2*(cc(1,ido,k,2)-cc(1,ido,k,4)))+ cc(1,ido,k,1) ch(1,ido,3,k) = cc(1,ido,k,1)-(hsqt2*(cc(1,ido,k,2)- cc(1,ido,k,4))) ch(1,1,2,k) = (-hsqt2*(cc(1,ido,k,2)+cc(1,ido,k,4)))- cc(1,ido,k,3) ch(1,1,4,k) = (-hsqt2*(cc(1,ido,k,2)+cc(1,ido,k,4)))+ cc(1,ido,k,3) end do return end subroutine r1f4kf ! !------------------------------------------------------- subroutine r1f5kb ( ido, l1, cc, in1, ch, in2, wa1, wa2, wa3, wa4 ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) arg real ( kind = 8 ) cc(in1,ido,5,l1) real ( kind = 8 ) ch(in2,ido,l1,5) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) ti11 real ( kind = 8 ) ti12 real ( kind = 8 ) tr11 real ( kind = 8 ) tr12 real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) real ( kind = 8 ) wa3(ido) real ( kind = 8 ) wa4(ido) arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 5.0E+00 tr11 = cos ( arg ) ti11 = sin ( arg ) tr12 = cos ( 2.0E+00 * arg ) ti12 = sin ( 2.0E+00 * arg ) do k = 1, l1 ch(1,1,k,1) = cc(1,1,1,k) + 2.0E+00 * cc(1,ido,2,k) & + 2.0E+00 * cc(1,ido,4,k) ch(1,1,k,2) = ( cc(1,1,1,k) & + tr11 * 2.0E+00 * cc(1,ido,2,k) + tr12 * 2.0E+00 * cc(1,ido,4,k) ) & - ( ti11 * 2.0E+00 * cc(1,1,3,k) + ti12 * 2.0E+00 * cc(1,1,5,k)) ch(1,1,k,3) = ( cc(1,1,1,k) & + tr12 * 2.0E+00 * cc(1,ido,2,k) + tr11 * 2.0E+00 * cc(1,ido,4,k) ) & - ( ti12 * 2.0E+00 * cc(1,1,3,k) - ti11 * 2.0E+00 * cc(1,1,5,k)) ch(1,1,k,4) = ( cc(1,1,1,k) & + tr12 * 2.0E+00 * cc(1,ido,2,k) + tr11 * 2.0E+00 * cc(1,ido,4,k) ) & + ( ti12 * 2.0E+00 * cc(1,1,3,k) - ti11 * 2.0E+00 * cc(1,1,5,k)) ch(1,1,k,5) = ( cc(1,1,1,k) & + tr11 * 2.0E+00 * cc(1,ido,2,k) + tr12 * 2.0E+00 * cc(1,ido,4,k) ) & + ( ti11 * 2.0E+00 * cc(1,1,3,k) + ti12 * 2.0E+00 * cc(1,1,5,k) ) end do if ( ido == 1 ) then return end if idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,k,1) = cc(1,i-1,1,k)+(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +(cc(1,i-1,5,k)+cc(1,ic-1,4,k)) ch(1,i,k,1) = cc(1,i,1,k)+(cc(1,i,3,k)-cc(1,ic,2,k)) & +(cc(1,i,5,k)-cc(1,ic,4,k)) ch(1,i-1,k,2) = wa1(i-2)*((cc(1,i-1,1,k)+tr11* & (cc(1,i-1,3,k)+cc(1,ic-1,2,k))+tr12 & *(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))-(ti11*(cc(1,i,3,k) & +cc(1,ic,2,k))+ti12*(cc(1,i,5,k)+cc(1,ic,4,k)))) & -wa1(i-1)*((cc(1,i,1,k)+tr11*(cc(1,i,3,k)-cc(1,ic,2,k)) & +tr12*(cc(1,i,5,k)-cc(1,ic,4,k)))+(ti11*(cc(1,i-1,3,k) & -cc(1,ic-1,2,k))+ti12*(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) ch(1,i,k,2) = wa1(i-2)*((cc(1,i,1,k)+tr11*(cc(1,i,3,k) & -cc(1,ic,2,k))+tr12*(cc(1,i,5,k)-cc(1,ic,4,k))) & +(ti11*(cc(1,i-1,3,k)-cc(1,ic-1,2,k))+ti12 & *(cc(1,i-1,5,k)-cc(1,ic-1,4,k))))+wa1(i-1) & *((cc(1,i-1,1,k)+tr11*(cc(1,i-1,3,k) & +cc(1,ic-1,2,k))+tr12*(cc(1,i-1,5,k)+cc(1,ic-1,4,k))) & -(ti11*(cc(1,i,3,k)+cc(1,ic,2,k))+ti12 & *(cc(1,i,5,k)+cc(1,ic,4,k)))) ch(1,i-1,k,3) = wa2(i-2) & *((cc(1,i-1,1,k)+tr12*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr11*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))-(ti12*(cc(1,i,3,k) & +cc(1,ic,2,k))-ti11*(cc(1,i,5,k)+cc(1,ic,4,k)))) & -wa2(i-1) & *((cc(1,i,1,k)+tr12*(cc(1,i,3,k)- & cc(1,ic,2,k))+tr11*(cc(1,i,5,k)-cc(1,ic,4,k))) & +(ti12*(cc(1,i-1,3,k)-cc(1,ic-1,2,k))-ti11 & *(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) ch(1,i,k,3) = wa2(i-2) & *((cc(1,i,1,k)+tr12*(cc(1,i,3,k)- & cc(1,ic,2,k))+tr11*(cc(1,i,5,k)-cc(1,ic,4,k))) & +(ti12*(cc(1,i-1,3,k)-cc(1,ic-1,2,k))-ti11 & *(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) & +wa2(i-1) & *((cc(1,i-1,1,k)+tr12*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr11*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))-(ti12*(cc(1,i,3,k) & +cc(1,ic,2,k))-ti11*(cc(1,i,5,k)+cc(1,ic,4,k)))) ch(1,i-1,k,4) = wa3(i-2) & *((cc(1,i-1,1,k)+tr12*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr11*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))+(ti12*(cc(1,i,3,k) & +cc(1,ic,2,k))-ti11*(cc(1,i,5,k)+cc(1,ic,4,k)))) & -wa3(i-1) & *((cc(1,i,1,k)+tr12*(cc(1,i,3,k)- & cc(1,ic,2,k))+tr11*(cc(1,i,5,k)-cc(1,ic,4,k))) & -(ti12*(cc(1,i-1,3,k)-cc(1,ic-1,2,k))-ti11 & *(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) ch(1,i,k,4) = wa3(i-2) & *((cc(1,i,1,k)+tr12*(cc(1,i,3,k)- & cc(1,ic,2,k))+tr11*(cc(1,i,5,k)-cc(1,ic,4,k))) & -(ti12*(cc(1,i-1,3,k)-cc(1,ic-1,2,k))-ti11 & *(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) & +wa3(i-1) & *((cc(1,i-1,1,k)+tr12*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr11*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))+(ti12*(cc(1,i,3,k) & +cc(1,ic,2,k))-ti11*(cc(1,i,5,k)+cc(1,ic,4,k)))) ch(1,i-1,k,5) = wa4(i-2) & *((cc(1,i-1,1,k)+tr11*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr12*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))+(ti11*(cc(1,i,3,k) & +cc(1,ic,2,k))+ti12*(cc(1,i,5,k)+cc(1,ic,4,k)))) & -wa4(i-1) & *((cc(1,i,1,k)+tr11*(cc(1,i,3,k)-cc(1,ic,2,k)) & +tr12*(cc(1,i,5,k)-cc(1,ic,4,k)))-(ti11*(cc(1,i-1,3,k) & -cc(1,ic-1,2,k))+ti12*(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) ch(1,i,k,5) = wa4(i-2) & *((cc(1,i,1,k)+tr11*(cc(1,i,3,k)-cc(1,ic,2,k)) & +tr12*(cc(1,i,5,k)-cc(1,ic,4,k)))-(ti11*(cc(1,i-1,3,k) & -cc(1,ic-1,2,k))+ti12*(cc(1,i-1,5,k)-cc(1,ic-1,4,k)))) & +wa4(i-1) & *((cc(1,i-1,1,k)+tr11*(cc(1,i-1,3,k)+cc(1,ic-1,2,k)) & +tr12*(cc(1,i-1,5,k)+cc(1,ic-1,4,k)))+(ti11*(cc(1,i,3,k) & +cc(1,ic,2,k))+ti12*(cc(1,i,5,k)+cc(1,ic,4,k)))) end do end do return end subroutine r1f5kb ! !------------------------------------------------------- subroutine r1f5kf ( ido, l1, cc, in1, ch, in2, wa1, wa2, wa3, wa4 ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) ido integer ( kind = 4 ) in1 integer ( kind = 4 ) in2 integer ( kind = 4 ) l1 real ( kind = 8 ) arg real ( kind = 8 ) cc(in1,ido,l1,5) real ( kind = 8 ) ch(in2,ido,5,l1) integer ( kind = 4 ) i integer ( kind = 4 ) ic integer ( kind = 4 ) idp2 integer ( kind = 4 ) k real ( kind = 8 ) ti11 real ( kind = 8 ) ti12 real ( kind = 8 ) tr11 real ( kind = 8 ) tr12 real ( kind = 8 ) wa1(ido) real ( kind = 8 ) wa2(ido) real ( kind = 8 ) wa3(ido) real ( kind = 8 ) wa4(ido) arg = 2.0E+00 * 4.0E+00 * atan ( 1.0E+00 ) / 5.0E+00 tr11 = cos ( arg ) ti11 = sin ( arg ) tr12 = cos ( 2.0E+00 * arg ) ti12 = sin ( 2.0E+00 * arg ) do k = 1, l1 ch(1,1,1,k) = cc(1,1,k,1) + ( cc(1,1,k,5) + cc(1,1,k,2) ) & + ( cc(1,1,k,4) + cc(1,1,k,3) ) ch(1,ido,2,k) = cc(1,1,k,1) + tr11 * ( cc(1,1,k,5) + cc(1,1,k,2) ) & + tr12 * ( cc(1,1,k,4) + cc(1,1,k,3) ) ch(1,1,3,k) = ti11 * ( cc(1,1,k,5) - cc(1,1,k,2) ) & + ti12 * ( cc(1,1,k,4) - cc(1,1,k,3) ) ch(1,ido,4,k) = cc(1,1,k,1) + tr12 * ( cc(1,1,k,5) + cc(1,1,k,2) ) & + tr11 * ( cc(1,1,k,4) + cc(1,1,k,3) ) ch(1,1,5,k) = ti12 * ( cc(1,1,k,5) - cc(1,1,k,2) ) & - ti11 * ( cc(1,1,k,4) - cc(1,1,k,3) ) end do if ( ido == 1 ) then return end if idp2 = ido + 2 do k = 1, l1 do i = 3, ido, 2 ic = idp2 - i ch(1,i-1,1,k) = cc(1,i-1,k,1)+((wa1(i-2)*cc(1,i-1,k,2)+ & wa1(i-1)*cc(1,i,k,2))+(wa4(i-2)*cc(1,i-1,k,5)+wa4(i-1)* & cc(1,i,k,5)))+((wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))+(wa3(i-2)*cc(1,i-1,k,4)+ & wa3(i-1)*cc(1,i,k,4))) ch(1,i,1,k) = cc(1,i,k,1)+((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)* & cc(1,i-1,k,5)))+((wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4))) ch(1,i-1,3,k) = cc(1,i-1,k,1)+tr11* & ( wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2) & +wa4(i-2)*cc(1,i-1,k,5)+wa4(i-1)*cc(1,i,k,5))+tr12* & ( wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3) & +wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)*cc(1,i,k,4))+ti11* & ( wa1(i-2)*cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2) & -(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)*cc(1,i-1,k,5)))+ti12* & ( wa2(i-2)*cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3) & -(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)*cc(1,i-1,k,4))) ch(1,ic-1,2,k) = cc(1,i-1,k,1)+tr11* & ( wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2) & +wa4(i-2)*cc(1,i-1,k,5)+wa4(i-1)*cc(1,i,k,5))+tr12* & ( wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3) & +wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)*cc(1,i,k,4))-(ti11* & ( wa1(i-2)*cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2) & -(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)*cc(1,i-1,k,5)))+ti12* & ( wa2(i-2)*cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3) & -(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)*cc(1,i-1,k,4)))) ch(1,i,3,k) = (cc(1,i,k,1)+tr11*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)* & cc(1,i-1,k,5)))+tr12*((wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4))))+(ti11*((wa4(i-2)*cc(1,i-1,k,5)+ & wa4(i-1)*cc(1,i,k,5))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))+ti12*((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3)))) ch(1,ic,2,k) = (ti11*((wa4(i-2)*cc(1,i-1,k,5)+wa4(i-1)* & cc(1,i,k,5))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))+ti12*((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))))-(cc(1,i,k,1)+tr11*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)* & cc(1,i-1,k,5)))+tr12*((wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4)))) ch(1,i-1,5,k) = (cc(1,i-1,k,1)+tr12*((wa1(i-2)* & cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2))+(wa4(i-2)* & cc(1,i-1,k,5)+wa4(i-1)*cc(1,i,k,5)))+tr11*((wa2(i-2)* & cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3))+(wa3(i-2)* & cc(1,i-1,k,4)+wa3(i-1)*cc(1,i,k,4))))+(ti12*((wa1(i-2)* & cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2))-(wa4(i-2)* & cc(1,i,k,5)-wa4(i-1)*cc(1,i-1,k,5)))-ti11*((wa2(i-2)* & cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3))-(wa3(i-2)* & cc(1,i,k,4)-wa3(i-1)*cc(1,i-1,k,4)))) ch(1,ic-1,4,k) = (cc(1,i-1,k,1)+tr12*((wa1(i-2)* & cc(1,i-1,k,2)+wa1(i-1)*cc(1,i,k,2))+(wa4(i-2)* & cc(1,i-1,k,5)+wa4(i-1)*cc(1,i,k,5)))+tr11*((wa2(i-2)* & cc(1,i-1,k,3)+wa2(i-1)*cc(1,i,k,3))+(wa3(i-2)* & cc(1,i-1,k,4)+wa3(i-1)*cc(1,i,k,4))))-(ti12*((wa1(i-2)* & cc(1,i,k,2)-wa1(i-1)*cc(1,i-1,k,2))-(wa4(i-2)* & cc(1,i,k,5)-wa4(i-1)*cc(1,i-1,k,5)))-ti11*((wa2(i-2)* & cc(1,i,k,3)-wa2(i-1)*cc(1,i-1,k,3))-(wa3(i-2)* & cc(1,i,k,4)-wa3(i-1)*cc(1,i-1,k,4)))) ch(1,i,5,k) = (cc(1,i,k,1)+tr12*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)* & cc(1,i-1,k,5)))+tr11*((wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4))))+(ti12*((wa4(i-2)*cc(1,i-1,k,5)+ & wa4(i-1)*cc(1,i,k,5))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))-ti11*((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3)))) ch(1,ic,4,k) = (ti12*((wa4(i-2)*cc(1,i-1,k,5)+wa4(i-1)* & cc(1,i,k,5))-(wa1(i-2)*cc(1,i-1,k,2)+wa1(i-1)* & cc(1,i,k,2)))-ti11*((wa3(i-2)*cc(1,i-1,k,4)+wa3(i-1)* & cc(1,i,k,4))-(wa2(i-2)*cc(1,i-1,k,3)+wa2(i-1)* & cc(1,i,k,3))))-(cc(1,i,k,1)+tr12*((wa1(i-2)*cc(1,i,k,2)- & wa1(i-1)*cc(1,i-1,k,2))+(wa4(i-2)*cc(1,i,k,5)-wa4(i-1)* & cc(1,i-1,k,5)))+tr11*((wa2(i-2)*cc(1,i,k,3)-wa2(i-1)* & cc(1,i-1,k,3))+(wa3(i-2)*cc(1,i,k,4)-wa3(i-1)* & cc(1,i-1,k,4)))) end do end do return end subroutine r1f5kf ! !------------------------------------------------------- subroutine r4_factor ( n, nf, fac ) !*****************************************************************************80 ! implicit none real ( kind = 8 ) fac(*) integer ( kind = 4 ) j integer ( kind = 4 ) n integer ( kind = 4 ) nf integer ( kind = 4 ) nl integer ( kind = 4 ) nq integer ( kind = 4 ) nr integer ( kind = 4 ) ntry nl = n nf = 0 j = 0 do while ( 1 < nl ) j = j + 1 if ( j == 1 ) then ntry = 4 else if ( j == 2 ) then ntry = 2 else if ( j == 3 ) then ntry = 3 else if ( j == 4 ) then ntry = 5 else ntry = ntry + 2 end if do nq = nl / ntry nr = nl - ntry * nq if ( nr /= 0 ) then exit end if nf = nf + 1 fac(nf) = real ( ntry, kind = 4 ) nl = nq end do end do return end subroutine r4_factor ! !----------------------------------------------------- subroutine rfftf1 ( n, in, c, ch, wa, fac ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) in integer ( kind = 4 ) n real ( kind = 8 ) c(in,*) real ( kind = 8 ) ch(*) real ( kind = 8 ) fac(15) integer ( kind = 4 ) idl1 integer ( kind = 4 ) ido integer ( kind = 4 ) ip integer ( kind = 4 ) iw integer ( kind = 4 ) ix2 integer ( kind = 4 ) ix3 integer ( kind = 4 ) ix4 integer ( kind = 4 ) j integer ( kind = 4 ) k1 integer ( kind = 4 ) kh integer ( kind = 4 ) l1 integer ( kind = 4 ) l2 integer ( kind = 4 ) modn integer ( kind = 4 ) na integer ( kind = 4 ) nf integer ( kind = 4 ) nl real ( kind = 8 ) sn real ( kind = 8 ) tsn real ( kind = 8 ) tsnm real ( kind = 8 ) wa(n) nf = int ( fac(2) ) na = 1 l2 = n iw = n do k1 = 1, nf kh = nf - k1 ip = int ( fac(kh+3) ) l1 = l2 / ip ido = n / l2 idl1 = ido * l1 iw = iw - ( ip - 1 ) * ido na = 1 - na if ( ip == 4 ) then ix2 = iw + ido ix3 = ix2 + ido if ( na == 0 ) then call r1f4kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3) ) else call r1f4kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3) ) end if else if ( ip == 2 ) then if ( na == 0 ) then call r1f2kf ( ido, l1, c, in, ch, 1, wa(iw) ) else call r1f2kf ( ido, l1, ch, 1, c, in, wa(iw) ) end if else if ( ip == 3 ) then ix2 = iw + ido if ( na == 0 ) then call r1f3kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2) ) else call r1f3kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2) ) end if else if ( ip == 5 ) then ix2 = iw + ido ix3 = ix2 + ido ix4 = ix3 + ido if ( na == 0 ) then call r1f5kf ( ido, l1, c, in, ch, 1, wa(iw), wa(ix2), wa(ix3), wa(ix4) ) else call r1f5kf ( ido, l1, ch, 1, c, in, wa(iw), wa(ix2), wa(ix3), wa(ix4) ) end if else if ( ido == 1 ) then na = 1 - na end if if ( na == 0 ) then call r1fgkf ( ido, ip, l1, idl1, c, c, c, in, ch, ch, 1, wa(iw) ) na = 1 else call r1fgkf ( ido, ip, l1, idl1, ch, ch, ch, 1, c, c, in, wa(iw) ) na = 0 end if end if l2 = l1 end do sn = 1.0E+00 / real ( n, kind = 4 ) tsn = 2.0E+00 / real ( n, kind = 4 ) tsnm = -tsn modn = mod ( n, 2 ) nl = n - 2 if ( modn /= 0 ) then nl = n - 1 end if if ( na == 0 ) then c(1,1) = sn * ch(1) do j = 2, nl, 2 c(1,j) = tsn * ch(j) c(1,j+1) = tsnm * ch(j+1) end do if ( modn == 0 ) then c(1,n) = sn * ch(n) end if else c(1,1) = sn * c(1,1) do j = 2, nl, 2 c(1,j) = tsn * c(1,j) c(1,j+1) = tsnm * c(1,j+1) end do if ( modn == 0 ) then c(1,n) = sn * c(1,n) end if end if return end subroutine rfftf1 ! !------------------------------------------------------ subroutine rffti1 ( n, wa, fac ) !*****************************************************************************80 ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) arg real ( kind = 8 ) argh real ( kind = 8 ) argld real ( kind = 8 ) fac(15) real ( kind = 8 ) fi integer ( kind = 4 ) i integer ( kind = 4 ) ib integer ( kind = 4 ) ido integer ( kind = 4 ) ii integer ( kind = 4 ) ip integer ( kind = 4 ) ipm integer ( kind = 4 ) is integer ( kind = 4 ) j integer ( kind = 4 ) k1 integer ( kind = 4 ) l1 integer ( kind = 4 ) l2 integer ( kind = 4 ) ld integer ( kind = 4 ) nf integer ( kind = 4 ) nfm1 integer ( kind = 4 ) nl integer ( kind = 4 ) nq integer ( kind = 4 ) nr integer ( kind = 4 ) ntry real ( kind = 8 ) tpi real ( kind = 8 ) wa(n) nl = n nf = 0 j = 0 do while ( 1 < nl ) j = j + 1 if ( j == 1 ) then ntry = 4 else if ( j == 2 ) then ntry = 2 else if ( j == 3 ) then ntry = 3 else if ( j == 4 ) then ntry = 5 else ntry = ntry + 2 end if do nq = nl / ntry nr = nl - ntry * nq if ( nr /= 0 ) then exit end if nf = nf + 1 fac(nf+2) = real ( ntry, kind = 4 ) nl = nq ! ! If 2 is a factor, make sure it appears first in the list of factors. ! if ( ntry == 2 ) then if ( nf /= 1 ) then do i = 2, nf ib = nf - i + 2 fac(ib+2) = fac(ib+1) end do fac(3) = 2.0E+00 end if end if end do end do fac(1) = real ( n, kind = 4 ) fac(2) = real ( nf, kind = 4 ) tpi = 8.0D+00 * atan ( 1.0D+00 ) argh = tpi / real ( n, kind = 4 ) is = 0 nfm1 = nf-1 l1 = 1 do k1 = 1, nfm1 ip = int ( fac(k1+2) ) ld = 0 l2 = l1 * ip ido = n / l2 ipm = ip - 1 do j = 1, ipm ld = ld + l1 i = is argld = real ( ld, kind = 4 ) * argh fi = 0.0E+00 do ii = 3, ido, 2 i = i + 2 fi = fi + 1.0E+00 arg = fi * argld wa(i-1) = real ( cos ( arg ), kind = 8 ) wa(i) = real ( sin ( arg ), kind = 8 ) end do is = is + ido end do l1 = l2 end do return end subroutine rffti1