fftsg3d.c

2维fft程序

    if (itnull != 0) {
        free(t);
    }
}


/* -------- child routines -------- */


void xdft3da_sub(int n1, int n2, int n3, int icr, int isgn, 
    double ***a, double *t, int *ip, double *w)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    void rdft(int n, int isgn, double *a, int *ip, double *w);
    int i, j, k;
    
    for (i = 0; i < n1; i++) {
        if (icr == 0) {
            for (j = 0; j < n2; j++) {
                cdft(n3, isgn, a[i][j], ip, w);
            }
        } else if (isgn >= 0) {
            for (j = 0; j < n2; j++) {
                rdft(n3, isgn, a[i][j], ip, w);
            }
        }
        if (n3 > 4) {
            for (k = 0; k < n3; k += 8) {
                for (j = 0; j < n2; j++) {
                    t[2 * j] = a[i][j][k];
                    t[2 * j + 1] = a[i][j][k + 1];
                    t[2 * n2 + 2 * j] = a[i][j][k + 2];
                    t[2 * n2 + 2 * j + 1] = a[i][j][k + 3];
                    t[4 * n2 + 2 * j] = a[i][j][k + 4];
                    t[4 * n2 + 2 * j + 1] = a[i][j][k + 5];
                    t[6 * n2 + 2 * j] = a[i][j][k + 6];
                    t[6 * n2 + 2 * j + 1] = a[i][j][k + 7];
                }
                cdft(2 * n2, isgn, t, ip, w);
                cdft(2 * n2, isgn, &t[2 * n2], ip, w);
                cdft(2 * n2, isgn, &t[4 * n2], ip, w);
                cdft(2 * n2, isgn, &t[6 * n2], ip, w);
                for (j = 0; j < n2; j++) {
                    a[i][j][k] = t[2 * j];
                    a[i][j][k + 1] = t[2 * j + 1];
                    a[i][j][k + 2] = t[2 * n2 + 2 * j];
                    a[i][j][k + 3] = t[2 * n2 + 2 * j + 1];
                    a[i][j][k + 4] = t[4 * n2 + 2 * j];
                    a[i][j][k + 5] = t[4 * n2 + 2 * j + 1];
                    a[i][j][k + 6] = t[6 * n2 + 2 * j];
                    a[i][j][k + 7] = t[6 * n2 + 2 * j + 1];
                }
            }
        } else if (n3 == 4) {
            for (j = 0; j < n2; j++) {
                t[2 * j] = a[i][j][0];
                t[2 * j + 1] = a[i][j][1];
                t[2 * n2 + 2 * j] = a[i][j][2];
                t[2 * n2 + 2 * j + 1] = a[i][j][3];
            }
            cdft(2 * n2, isgn, t, ip, w);
            cdft(2 * n2, isgn, &t[2 * n2], ip, w);
            for (j = 0; j < n2; j++) {
                a[i][j][0] = t[2 * j];
                a[i][j][1] = t[2 * j + 1];
                a[i][j][2] = t[2 * n2 + 2 * j];
                a[i][j][3] = t[2 * n2 + 2 * j + 1];
            }
        } else if (n3 == 2) {
            for (j = 0; j < n2; j++) {
                t[2 * j] = a[i][j][0];
                t[2 * j + 1] = a[i][j][1];
            }
            cdft(2 * n2, isgn, t, ip, w);
            for (j = 0; j < n2; j++) {
                a[i][j][0] = t[2 * j];
                a[i][j][1] = t[2 * j + 1];
            }
        }
        if (icr != 0 && isgn < 0) {
            for (j = 0; j < n2; j++) {
                rdft(n3, isgn, a[i][j], ip, w);
            }
        }
    }
}


void cdft3db_sub(int n1, int n2, int n3, int isgn, double ***a, 
    double *t, int *ip, double *w)
{
    void cdft(int n, int isgn, double *a, int *ip, double *w);
    int i, j, k;
    
    if (n3 > 4) {
        for (j = 0; j < n2; j++) {
            for (k = 0; k < n3; k += 8) {
                for (i = 0; i < n1; i++) {
                    t[2 * i] = a[i][j][k];
                    t[2 * i + 1] = a[i][j][k + 1];
                    t[2 * n1 + 2 * i] = a[i][j][k + 2];
                    t[2 * n1 + 2 * i + 1] = a[i][j][k + 3];
                    t[4 * n1 + 2 * i] = a[i][j][k + 4];
                    t[4 * n1 + 2 * i + 1] = a[i][j][k + 5];
                    t[6 * n1 + 2 * i] = a[i][j][k + 6];
                    t[6 * n1 + 2 * i + 1] = a[i][j][k + 7];
                }
                cdft(2 * n1, isgn, t, ip, w);
                cdft(2 * n1, isgn, &t[2 * n1], ip, w);
                cdft(2 * n1, isgn, &t[4 * n1], ip, w);
                cdft(2 * n1, isgn, &t[6 * n1], ip, w);
                for (i = 0; i < n1; i++) {
                    a[i][j][k] = t[2 * i];
                    a[i][j][k + 1] = t[2 * i + 1];
                    a[i][j][k + 2] = t[2 * n1 + 2 * i];
                    a[i][j][k + 3] = t[2 * n1 + 2 * i + 1];
                    a[i][j][k + 4] = t[4 * n1 + 2 * i];
                    a[i][j][k + 5] = t[4 * n1 + 2 * i + 1];
                    a[i][j][k + 6] = t[6 * n1 + 2 * i];
                    a[i][j][k + 7] = t[6 * n1 + 2 * i + 1];
                }
            }
        }
    } else if (n3 == 4) {
        for (j = 0; j < n2; j++) {
            for (i = 0; i < n1; i++) {
                t[2 * i] = a[i][j][0];
                t[2 * i + 1] = a[i][j][1];
                t[2 * n1 + 2 * i] = a[i][j][2];
                t[2 * n1 + 2 * i + 1] = a[i][j][3];
            }
            cdft(2 * n1, isgn, t, ip, w);
            cdft(2 * n1, isgn, &t[2 * n1], ip, w);
            for (i = 0; i < n1; i++) {
                a[i][j][0] = t[2 * i];
                a[i][j][1] = t[2 * i + 1];
                a[i][j][2] = t[2 * n1 + 2 * i];
                a[i][j][3] = t[2 * n1 + 2 * i + 1];
            }
        }
    } else if (n3 == 2) {
        for (j = 0; j < n2; j++) {
            for (i = 0; i < n1; i++) {
                t[2 * i] = a[i][j][0];
                t[2 * i + 1] = a[i][j][1];
            }
            cdft(2 * n1, isgn, t, ip, w);
            for (i = 0; i < n1; i++) {
                a[i][j][0] = t[2 * i];
                a[i][j][1] = t[2 * i + 1];
            }
        }
    }
}


void rdft3d_sub(int n1, int n2, int n3, int isgn, double ***a)
{
    int n1h, n2h, i, j, k, l;
    double xi;
    
    n1h = n1 >> 1;
    n2h = n2 >> 1;
    if (isgn < 0) {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            xi = a[i][0][0] - a[j][0][0];
            a[i][0][0] += a[j][0][0];
            a[j][0][0] = xi;
            xi = a[j][0][1] - a[i][0][1];
            a[i][0][1] += a[j][0][1];
            a[j][0][1] = xi;
            xi = a[i][n2h][0] - a[j][n2h][0];
            a[i][n2h][0] += a[j][n2h][0];
            a[j][n2h][0] = xi;
            xi = a[j][n2h][1] - a[i][n2h][1];
            a[i][n2h][1] += a[j][n2h][1];
            a[j][n2h][1] = xi;
            for (k = 1; k < n2h; k++) {
                l = n2 - k;
                xi = a[i][k][0] - a[j][l][0];
                a[i][k][0] += a[j][l][0];
                a[j][l][0] = xi;
                xi = a[j][l][1] - a[i][k][1];
                a[i][k][1] += a[j][l][1];
                a[j][l][1] = xi;
                xi = a[j][k][0] - a[i][l][0];
                a[j][k][0] += a[i][l][0];
                a[i][l][0] = xi;
                xi = a[i][l][1] - a[j][k][1];
                a[j][k][1] += a[i][l][1];
                a[i][l][1] = xi;
            }
        }
        for (k = 1; k < n2h; k++) {
            l = n2 - k;
            xi = a[0][k][0] - a[0][l][0];
            a[0][k][0] += a[0][l][0];
            a[0][l][0] = xi;
            xi = a[0][l][1] - a[0][k][1];
            a[0][k][1] += a[0][l][1];
            a[0][l][1] = xi;
            xi = a[n1h][k][0] - a[n1h][l][0];
            a[n1h][k][0] += a[n1h][l][0];
            a[n1h][l][0] = xi;
            xi = a[n1h][l][1] - a[n1h][k][1];
            a[n1h][k][1] += a[n1h][l][1];
            a[n1h][l][1] = xi;
        }
    } else {
        for (i = 1; i < n1h; i++) {
            j = n1 - i;
            a[j][0][0] = 0.5 * (a[i][0][0] - a[j][0][0]);
            a[i][0][0] -= a[j][0][0];
            a[j][0][1] = 0.5 * (a[i][0][1] + a[j][0][1]);
            a[i][0][1] -= a[j][0][1];
            a[j][n2h][0] = 0.5 * (a[i][n2h][0] - a[j][n2h][0]);
            a[i][n2h][0] -= a[j][n2h][0];
            a[j][n2h][1] = 0.5 * (a[i][n2h][1] + a[j][n2h][1]);
            a[i][n2h][1] -= a[j][n2h][1];
            for (k = 1; k < n2h; k++) {
                l = n2 - k;
                a[j][l][0] = 0.5 * (a[i][k][0] - a[j][l][0]);
                a[i][k][0] -= a[j][l][0];
                a[j][l][1] = 0.5 * (a[i][k][1] + a[j][l][1]);
                a[i][k][1] -= a[j][l][1];
                a[i][l][0] = 0.5 * (a[j][k][0] - a[i][l][0]);
                a[j][k][0] -= a[i][l][0];
                a[i][l][1] = 0.5 * (a[j][k][1] + a[i][l][1]);
                a[j][k][1] -= a[i][l][1];
            }
        }
        for (k = 1; k < n2h; k++) {
            l = n2 - k;
            a[0][l][0] = 0.5 * (a[0][k][0] - a[0][l][0]);
            a[0][k][0] -= a[0][l][0];
            a[0][l][1] = 0.5 * (a[0][k][1] + a[0][l][1]);
            a[0][k][1] -= a[0][l][1];
            a[n1h][l][0] = 0.5 * (a[n1h][k][0] - a[n1h][l][0]);
            a[n1h][k][0] -= a[n1h][l][0];
            a[n1h][l][1] = 0.5 * (a[n1h][k][1] + a[n1h][l][1]);
            a[n1h][k][1] -= a[n1h][l][1];
        }
    }
}


void ddxt3da_sub(int n1, int n2, int n3, int ics, int isgn, 
    double ***a, double *t, int *ip, double *w)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int i, j, k;
    
    for (i = 0; i < n1; i++) {
        if (ics == 0) {
            for (j = 0; j < n2; j++) {
                ddct(n3, isgn, a[i][j], ip, w);
            }
        } else {
            for (j = 0; j < n2; j++) {
                ddst(n3, isgn, a[i][j], ip, w);
            }
        }
        if (n3 > 2) {
            for (k = 0; k < n3; k += 4) {
                for (j = 0; j < n2; j++) {
                    t[j] = a[i][j][k];
                    t[n2 + j] = a[i][j][k + 1];
                    t[2 * n2 + j] = a[i][j][k + 2];
                    t[3 * n2 + j] = a[i][j][k + 3];
                }
                if (ics == 0) {
                    ddct(n2, isgn, t, ip, w);
                    ddct(n2, isgn, &t[n2], ip, w);
                    ddct(n2, isgn, &t[2 * n2], ip, w);
                    ddct(n2, isgn, &t[3 * n2], ip, w);
                } else {
                    ddst(n2, isgn, t, ip, w);
                    ddst(n2, isgn, &t[n2], ip, w);
                    ddst(n2, isgn, &t[2 * n2], ip, w);
                    ddst(n2, isgn, &t[3 * n2], ip, w);
                }
                for (j = 0; j < n2; j++) {
                    a[i][j][k] = t[j];
                    a[i][j][k + 1] = t[n2 + j];
                    a[i][j][k + 2] = t[2 * n2 + j];
                    a[i][j][k + 3] = t[3 * n2 + j];
                }
            }
        } else if (n3 == 2) {
            for (j = 0; j < n2; j++) {
                t[j] = a[i][j][0];
                t[n2 + j] = a[i][j][1];
            }
            if (ics == 0) {
                ddct(n2, isgn, t, ip, w);
                ddct(n2, isgn, &t[n2], ip, w);
            } else {
                ddst(n2, isgn, t, ip, w);
                ddst(n2, isgn, &t[n2], ip, w);
            }
            for (j = 0; j < n2; j++) {
                a[i][j][0] = t[j];
                a[i][j][1] = t[n2 + j];
            }
        }
    }
}


void ddxt3db_sub(int n1, int n2, int n3, int ics, int isgn, 
    double ***a, double *t, int *ip, double *w)
{
    void ddct(int n, int isgn, double *a, int *ip, double *w);
    void ddst(int n, int isgn, double *a, int *ip, double *w);
    int i, j, k;
    
    if (n3 > 2) {
        for (j = 0; j < n2; j++) {
            for (k = 0; k < n3; k += 4) {
                for (i = 0; i < n1; i++) {
                    t[i] = a[i][j][k];
                    t[n1 + i] = a[i][j][k + 1];
                    t[2 * n1 + i] = a[i][j][k + 2];
                    t[3 * n1 + i] = a[i][j][k + 3];
                }
                if (ics == 0) {
                    ddct(n1, isgn, t, ip, w);
                    ddct(n1, isgn, &t[n1], ip, w);
                    ddct(n1, isgn, &t[2 * n1], ip, w);
                    ddct(n1, isgn, &t[3 * n1], ip, w);
                } else {
                    ddst(n1, isgn, t, ip, w);
                    ddst(n1, isgn, &t[n1], ip, w);
                    ddst(n1, isgn, &t[2 * n1], ip, w);
                    ddst(n1, isgn, &t[3 * n1], ip, w);
                }
                for (i = 0; i < n1; i++) {
                    a[i][j][k] = t[i];
                    a[i][j][k + 1] = t[n1 + i];
                    a[i][j][k + 2] = t[2 * n1 + i];
                    a[i][j][k + 3] = t[3 * n1 + i];
                }
            }
        }
    } else if (n3 == 2) {
        for (j = 0; j < n2; j++) {
            for (i = 0; i < n1; i++) {
                t[i] = a[i][j][0];
                t[n1 + i] = a[i][j][1];
            }
            if (ics == 0) {
                ddct(n1, isgn, t, ip, w);
                ddct(n1, isgn, &t[n1], ip, w);
            } else {
                ddst(n1, isgn, t, ip, w);
                ddst(n1, isgn, &t[n1], ip, w);
            }
            for (i = 0; i < n1; i++) {
                a[i][j][0] = t[i];
                a[i][j][1] = t[n1 + i];
            }
        }
    }
}


#ifdef USE_FFT3D_THREADS
struct fft3d_arg_st {
    int nthread;
    int n0;
    int n1;
    int n2;
    int n3;
    int ic;
    int isgn;
    double ***a;
    double *t;
    int *ip;
    double *w;
};
typedef struct fft3d_arg_st fft3d_arg_t;


void xdft3da_subth(int n1, int n2, int n3, int icr, int isgn, 
    double ***a, double *t, int *ip, double *w)
{
    void *xdft3da_th(void *p);
    fft3d_thread_t th[FFT3D_MAX_THREADS];
    fft3d_arg_t ag[FFT3D_MAX_THREADS];
    int nthread, nt, i;
    
    nthread = FFT3D_MAX_THREADS;
    if (nthread > n1) {
        nthread = n1;
    }
    nt = 8 * n2;
    if (n3 == 4) {
        nt >>= 1;
    } else if (n3 < 4) {
        nt >>= 2;
    }
    for (i = 0; i < nthread; i++) {
        ag[i].nthread = nthread;
        ag[i].n0 = i;
        ag[i].n1 = n1;
        ag[i].n2 = n2;
        ag[i].n3 = n3;
        ag[i].ic = icr;
        ag[i].isgn = isgn;
        ag[i].a = a;
        ag[i].t = &t[nt * i];
        ag[i].ip = ip;
        ag[i].w = w;
        fft3d_thread_create(&th[i], xdft3da_th, &ag[i]);
    }
    for (i = 0; i < nthread; i++) {
        fft3d_thread_wait(th[i]);
    }
}


void cdft3db_subth(int n1, int n2, int n3, int isgn, double ***a,