use fftautocorr

OneDQuantum/1DSchrodinger.c

@@ -17,7 +17,7 @@
 #endif 
 #include "science.h"
 #include "band.h"
-#include "pocketfft/pocketfft.h"
+#include "fftautocorr/fftautocorr.h"
 
 #define Y_EPS 0.1 /**< Starting point for ode solver. Default value = 0.1 */
 
@@ -560,7 +560,6 @@ double LOphononScatter(double step, numpyint N, double kl,
     int i;
     double powerUnit = -kl*step*ANG;
     double *psiij;
-    rfft_plan p;
     int start, end;
     for(start = 0; start < N && 
             (fabs(psi_i[start]) < MINPSI || fabs(psi_j[start]) < MINPSI);
@@ -572,33 +571,21 @@ double LOphononScatter(double step, numpyint N, double kl,
         return 0.0;
     end += 1;
     /* end - start is at least 1 */ 
-    psiij = (double *)malloc(sizeof(double)*(end-start)*2);
-    p = make_rfft_plan((end-start)*2);
+    psiij = (double *)malloc(sizeof(double)*(end-start));
     for(i=start; i<end; i++) {
         psiij[i-start] = psi_i[i] * psi_j[i];
     }
-    for(i=0; i<end-start; i++) {
-        psiij[end+i-start] = 0;
-    }
-    end = (end - start);
-    if(rfft_forward(p, psiij, 1.0) != 0) {
+    end = end - start;
+    if(autocorr(psiij, end) != 0) {
         #ifdef _DEBUG
         printf("FFT failed!\n");
         #endif
         return -1.0;
     }
-    psiij[0] = sq(psiij[0]);
-    for(i=1; i<end; i++) {
-        psiij[2*i-1] = sq(psiij[2*i]) + sq(psiij[2*i-1]);
-        psiij[2*i] = 0;
-    }
-    psiij[2*end-1] = sq(psiij[2*end-1]);
-    rfft_backward(p, psiij, 0.5/end);
-    Iij += psiij[0]/2;
+    Iij += psiij[0] / 2;
     for(i=1; i<end; i++) {
         Iij += psiij[i]*exp(powerUnit*i);
     }
-    destroy_rfft_plan(p);
     free(psiij);
     return 2*Iij * sq(step);
 }
@@ -627,20 +614,20 @@ double LOtotal(double step, numpyint N, const double *kls,
     double Iij = 0;
     int j,n;
     int starti, endi;
-    double *psiij;
-    rfft_plan p;
+    autocorr_plan plan;
     for(starti = 0; starti < N && fabs(psi_i[starti]) < MINPSI; starti++);
     for(endi = N-1; endi >= starti && fabs(psi_i[endi]) < MINPSI; endi--);
     if(starti == endi)
         return 0.0;
     endi = endi - starti + 1;
     psi_i += starti;
-    p = make_rfft_plan(endi*2);
+    plan = make_autocorr_plan(endi);
 #ifdef __MP
 #pragma omp parallel
 #endif
     {
-        psiij = (double *)malloc(sizeof(double) * endi * 2);
+        double *psiij = (double *)malloc(sizeof(double) * endi * 2);
+        double *mempool = (double *) malloc(sizeof(double) * mem_len(plan));
         #ifdef __MP
         #pragma omp for reduction(+:Iij)
         #endif
@@ -650,10 +637,7 @@ double LOtotal(double step, numpyint N, const double *kls,
             for(j=0; j<endi; j++) {
                 psiij[j] = psi_i[j] * psi_j[j];
             }
-            for(j=0; j<endi; j++) {
-                psiij[endi+j] = 0;
-            }
-            if(rfft_forward(p, psiij, 1.0) != 0) {
+            if(autocorr_mem(plan, psiij, mempool) != 0) {
                 #ifdef _DEBUG
                 printf("FFT failed!\n");
                 #endif
@@ -664,23 +648,17 @@ double LOtotal(double step, numpyint N, const double *kls,
                 return -1.0;
                 #endif
             }
-            psiij[0] = sq(psiij[0]);
-            for(j=1; j<endi; j++) {
-                psiij[2*j-1] = sq(psiij[2*j]) + sq(psiij[2*j-1]);
-                psiij[2*j] = 0;
-            }
-            psiij[2*endi-1] = sq(psiij[2*endi-1]);
-            rfft_backward(p, psiij, 0.5/endi);
             Iij += psiij[0]/2;
             for(j=1; j<endi; j++) {
                 Iij += psiij[j]*exp(powerUnit*j);
             }
         }
     free(psiij);
+    free(mempool);
     }
     if(Iij < 0)
         return -1;
-    destroy_rfft_plan(p);
+    destroy_autocorr_plan(plan);
     return 2*Iij * sq(step);
 }
 

OneDQuantum/Makefile

@@ -24,27 +24,30 @@ endif
 all: 1DSchrodinger.$(EXT) 1DThermal.$(EXT) 1DMaxwell.$(EXT)
 MP:  1DSchrodinger_MP.$(EXT)
 
-1DSchrodinger_MP.$(EXT) : 1DSchrodinger_MP.o band.o pocketfft.o
+1DSchrodinger_MP.$(EXT) : 1DSchrodinger_MP.o band.o fftautocorr.o
 	$(CC) -shared -fPIC -fopenmp $^ -o $@ -lm
 
-1DSchrodinger.$(EXT) : 1DSchrodinger.o band.o pocketfft.o
+1DSchrodinger.$(EXT) : 1DSchrodinger.o band.o fftautocorr.o
 	$(CC) -shared -fPIC $^ -o $@ -lm
 
 %.$(EXT) : %.o
 	$(CC) -shared -fPIC $^ -o $@ -lm
 
-1DSchrodinger.o : 1DSchrodinger.c science.h band.h
+1DSchrodinger.o : 1DSchrodinger.c science.h band.h fftautocorr/fftautocorr.h
 	$(CC) $(CFLAGS) -c $< -o $@
 
-1DSchrodinger_MP.o : 1DSchrodinger.c science.h band.h
+1DSchrodinger_MP.o : 1DSchrodinger.c science.h band.h fftautocorr/fftautocorr.h
 	$(CC) $(CFLAGS) -fopenmp -D __MP -c $< -o $@
 
-pocketfft.o : pocketfft/pocketfft.c pocketfft/pocketfft.h
+fftautocorr.o : fftautocorr/fftautocorr.c fftautocorr/fftautocorr.h fftautocorr/factortable.h
 	$(CC) $(CFLAGS) -c $< -o $@
 
+fftautocorr/% : 
+	git submodule update --init
+
 %.o : %.c science.h 
 	$(CC) $(CFLAGS) -c $< -o $@
 
 .PHONY : clean
 clean :
-	@$(RM) {1DSchrodinger,1DThermal,1DMaxwell,1DSchrodinger_MP}.{so,o,dll,dylib} band.o
+	@$(RM) {1DSchrodinger,1DThermal,1DMaxwell,1DSchrodinger_MP}.{so,o,dll,dylib} band.o fftautocorr.o