Commit 6fe37363 authored by Leo Singer's avatar Leo Singer

New upstream version 3.50.0

parent 01b54810
This diff is collapsed.
...@@ -9,13 +9,9 @@ src_c_utils = \ ...@@ -9,13 +9,9 @@ src_c_utils = \
c_utils/walltime_c.h \ c_utils/walltime_c.h \
c_utils/trig_utils.c \ c_utils/trig_utils.c \
c_utils/trig_utils.h c_utils/trig_utils.h
src_libfftpack = \ src_libpocketfft = \
libfftpack/bluestein.h \ pocketfft/pocketfft.c \
libfftpack/bluestein.c \ pocketfft/pocketfft.h
libfftpack/fftpack.h \
libfftpack/fftpack.c \
libfftpack/ls_fft.h \
libfftpack/ls_fft.c
src_libsharp = \ src_libsharp = \
libsharp/sharp.c \ libsharp/sharp.c \
libsharp/sharp.h \ libsharp/sharp.h \
...@@ -24,6 +20,7 @@ src_libsharp = \ ...@@ -24,6 +20,7 @@ src_libsharp = \
libsharp/sharp_almhelpers.h \ libsharp/sharp_almhelpers.h \
libsharp/sharp_complex_hacks.h \ libsharp/sharp_complex_hacks.h \
libsharp/sharp_core.c \ libsharp/sharp_core.c \
libsharp/sharp_core_avx.c \
libsharp/sharp_core.h \ libsharp/sharp_core.h \
libsharp/sharp_geomhelpers.c \ libsharp/sharp_geomhelpers.c \
libsharp/sharp_geomhelpers.h \ libsharp/sharp_geomhelpers.h \
...@@ -115,11 +112,13 @@ src_healpix_cxx_fits= \ ...@@ -115,11 +112,13 @@ src_healpix_cxx_fits= \
Healpix_cxx/moc_fitsio.cc Healpix_cxx/moc_fitsio.cc
EXTRA_DIST = \ EXTRA_DIST = \
libfftpack/fftpack_inc.c \ libsharp/sharp_core_inc0.c \
libsharp/sharp_core_inc.c \ libsharp/sharp_core_inc.c \
libsharp/sharp_core_inc2.c \ libsharp/sharp_core_inc2.c \
libsharp/sharp_core_inchelper.c \ libsharp/sharp_core_inchelper.c \
cxxsupport/font_data.inc cxxsupport/font_data.inc \
pocketfft/LICENSE.md
include_HEADERS = \ include_HEADERS = \
libsharp/sharp.h \ libsharp/sharp.h \
libsharp/sharp_lowlevel.h \ libsharp/sharp_lowlevel.h \
...@@ -158,9 +157,9 @@ include_HEADERS = \ ...@@ -158,9 +157,9 @@ include_HEADERS = \
Healpix_cxx/moc_query.h \ Healpix_cxx/moc_query.h \
Healpix_cxx/weight_utils.h Healpix_cxx/weight_utils.h
libhealpix_cxx_la_SOURCES = $(src_c_utils) $(src_libfftpack) $(src_libsharp) $(src_cxxsupport) $(src_cxxsupport_fits) $(src_healpix_cxx) $(src_healpix_cxx_fits) $(src_healpix_cxx_mod) libhealpix_cxx_la_SOURCES = $(src_c_utils) $(src_libpocketfft) $(src_libsharp) $(src_cxxsupport) $(src_cxxsupport_fits) $(src_healpix_cxx) $(src_healpix_cxx_fits) $(src_healpix_cxx_mod)
libhealpix_cxx_la_LIBADD = $(CFITSIO_LIBS) libhealpix_cxx_la_LIBADD = $(CFITSIO_LIBS)
libhealpix_cxx_la_LDFLAGS = -version-info 1:0:0 libhealpix_cxx_la_LDFLAGS = -version-info 2:1:0
bin_PROGRAMS = alm2map_cxx anafast_cxx calc_powspec hotspots_cxx map2tga \ bin_PROGRAMS = alm2map_cxx anafast_cxx calc_powspec hotspots_cxx map2tga \
median_filter_cxx mult_alm rotalm_cxx smoothing_cxx syn_alm_cxx udgrade_cxx \ median_filter_cxx mult_alm rotalm_cxx smoothing_cxx syn_alm_cxx udgrade_cxx \
...@@ -199,7 +198,7 @@ hpxtest_LDADD = libhealpix_cxx.la ...@@ -199,7 +198,7 @@ hpxtest_LDADD = libhealpix_cxx.la
TESTS=hpxtest TESTS=hpxtest
AM_CPPFLAGS = $(CFITSIO_CFLAGS) -I$(top_srcdir)/c_utils -I$(top_srcdir)/libfftpack -I$(top_srcdir)/libsharp -I$(top_srcdir)/cxxsupport AM_CPPFLAGS = $(CFITSIO_CFLAGS) -I$(top_srcdir)/c_utils -I$(top_srcdir) -I$(top_srcdir)/libsharp -I$(top_srcdir)/cxxsupport
pkgconfigdir = $(libdir)/pkgconfig pkgconfigdir = $(libdir)/pkgconfig
nodist_pkgconfig_DATA = @PACKAGE_NAME@.pc nodist_pkgconfig_DATA = @PACKAGE_NAME@.pc
......
This diff is collapsed.
Changes since 3.40:
- When compiled with a recent gcc (>=5.0) for an x86_64 target, a special
version for AVX-capable CPUs will always be compiled in addition to the more
portable SSE2 version. This means that whenever AVX is available on a CPU,
SHTs will run up to twice as fast than they did before (unless flags like
"-march=native" were used).
- the FFT library used internally by libsharp was substantially redesigned.
Changes since 3.30: Changes since 3.30:
- IMPORTANT: the syntax for specifying ring weights and pixel windows has - IMPORTANT: the syntax for specifying ring weights and pixel windows has
changed! This affects the facilities anafast_cxx, smoothing_cxx, changed! This affects the facilities anafast_cxx, smoothing_cxx,
......
This diff is collapsed.
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
me=ar-lib me=ar-lib
scriptversion=2012-03-01.08; # UTC scriptversion=2012-03-01.08; # UTC
# Copyright (C) 2010-2017 Free Software Foundation, Inc. # Copyright (C) 2010-2018 Free Software Foundation, Inc.
# Written by Peter Rosin <peda@lysator.liu.se>. # Written by Peter Rosin <peda@lysator.liu.se>.
# #
# This program is free software; you can redistribute it and/or modify # This program is free software; you can redistribute it and/or modify
...@@ -18,7 +18,7 @@ scriptversion=2012-03-01.08; # UTC ...@@ -18,7 +18,7 @@ scriptversion=2012-03-01.08; # UTC
# GNU General Public License for more details. # GNU General Public License for more details.
# #
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
# As a special exception to the GNU General Public License, if you # As a special exception to the GNU General Public License, if you
# distribute this file as part of a program that contains a # distribute this file as part of a program that contains a
......
#! /bin/sh #! /bin/sh
# Wrapper for compilers which do not understand '-c -o'. # Wrapper for compilers which do not understand '-c -o'.
scriptversion=2012-10-14.11; # UTC scriptversion=2018-03-07.03; # UTC
# Copyright (C) 1999-2014 Free Software Foundation, Inc. # Copyright (C) 1999-2018 Free Software Foundation, Inc.
# Written by Tom Tromey <tromey@cygnus.com>. # Written by Tom Tromey <tromey@cygnus.com>.
# #
# This program is free software; you can redistribute it and/or modify # This program is free software; you can redistribute it and/or modify
...@@ -17,7 +17,7 @@ scriptversion=2012-10-14.11; # UTC ...@@ -17,7 +17,7 @@ scriptversion=2012-10-14.11; # UTC
# GNU General Public License for more details. # GNU General Public License for more details.
# #
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
# As a special exception to the GNU General Public License, if you # As a special exception to the GNU General Public License, if you
# distribute this file as part of a program that contains a # distribute this file as part of a program that contains a
...@@ -255,7 +255,8 @@ EOF ...@@ -255,7 +255,8 @@ EOF
echo "compile $scriptversion" echo "compile $scriptversion"
exit $? exit $?
;; ;;
cl | *[/\\]cl | cl.exe | *[/\\]cl.exe ) cl | *[/\\]cl | cl.exe | *[/\\]cl.exe | \
icl | *[/\\]icl | icl.exe | *[/\\]icl.exe )
func_cl_wrapper "$@" # Doesn't return... func_cl_wrapper "$@" # Doesn't return...
;; ;;
esac esac
...@@ -339,9 +340,9 @@ exit $ret ...@@ -339,9 +340,9 @@ exit $ret
# Local Variables: # Local Variables:
# mode: shell-script # mode: shell-script
# sh-indentation: 2 # sh-indentation: 2
# eval: (add-hook 'write-file-hooks 'time-stamp) # eval: (add-hook 'before-save-hook 'time-stamp)
# time-stamp-start: "scriptversion=" # time-stamp-start: "scriptversion="
# time-stamp-format: "%:y-%02m-%02d.%02H" # time-stamp-format: "%:y-%02m-%02d.%02H"
# time-stamp-time-zone: "UTC" # time-stamp-time-zone: "UTC0"
# time-stamp-end: "; # UTC" # time-stamp-end: "; # UTC"
# End: # End:
This diff is collapsed.
AC_INIT([healpix_cxx], [3.40.0]) AC_INIT([healpix_cxx], [3.50.0])
AM_INIT_AUTOMAKE([foreign subdir-objects -Wall -Werror]) AM_INIT_AUTOMAKE([foreign subdir-objects -Wall -Werror])
AM_MAINTAINER_MODE([enable]) AM_MAINTAINER_MODE([enable])
......
...@@ -58,18 +58,12 @@ void openmp_status() ...@@ -58,18 +58,12 @@ void openmp_status()
#endif #endif
} }
extern "C" {
int sharp_veclen(void);
}
void vec_status() void vec_status()
{ {
cout << "Vector math: "; cout << "Supported vector length: " << sharp_veclen() << endl;
#if(defined(__AVX__))
cout << "AVX" << endl;
#elif(defined(__SSE2__))
cout << "SSE2" << endl;
#elif(defined(__SSE__))
cout << "SSE" << endl;
#else
cout << "not supported by this binary" << endl;
#endif
} }
} //unnamed namespace } //unnamed namespace
...@@ -77,7 +71,7 @@ void vec_status() ...@@ -77,7 +71,7 @@ void vec_status()
void announce (const string &name) void announce (const string &name)
{ {
#ifndef VERSION #ifndef VERSION
#define VERSION "?.?" #define VERSION "3.50"
#endif #endif
string version = "v" VERSION; string version = "v" VERSION;
string name2 = name+" "+version; string name2 = name+" "+version;
......
#! /bin/sh #! /bin/sh
# depcomp - compile a program generating dependencies as side-effects # depcomp - compile a program generating dependencies as side-effects
scriptversion=2016-01-11.22; # UTC scriptversion=2018-03-07.03; # UTC
# Copyright (C) 1999-2017 Free Software Foundation, Inc. # Copyright (C) 1999-2018 Free Software Foundation, Inc.
# This program is free software; you can redistribute it and/or modify # This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by # it under the terms of the GNU General Public License as published by
...@@ -16,7 +16,7 @@ scriptversion=2016-01-11.22; # UTC ...@@ -16,7 +16,7 @@ scriptversion=2016-01-11.22; # UTC
# GNU General Public License for more details. # GNU General Public License for more details.
# You should have received a copy of the GNU General Public License # You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
# As a special exception to the GNU General Public License, if you # As a special exception to the GNU General Public License, if you
# distribute this file as part of a program that contains a # distribute this file as part of a program that contains a
...@@ -783,7 +783,7 @@ exit 0 ...@@ -783,7 +783,7 @@ exit 0
# Local Variables: # Local Variables:
# mode: shell-script # mode: shell-script
# sh-indentation: 2 # sh-indentation: 2
# eval: (add-hook 'write-file-hooks 'time-stamp) # eval: (add-hook 'before-save-hook 'time-stamp)
# time-stamp-start: "scriptversion=" # time-stamp-start: "scriptversion="
# time-stamp-format: "%:y-%02m-%02d.%02H" # time-stamp-format: "%:y-%02m-%02d.%02H"
# time-stamp-time-zone: "UTC0" # time-stamp-time-zone: "UTC0"
......
#!/bin/sh #!/bin/sh
# install - install a program, script, or datafile # install - install a program, script, or datafile
scriptversion=2014-09-12.12; # UTC scriptversion=2018-03-11.20; # UTC
# This originates from X11R5 (mit/util/scripts/install.sh), which was # This originates from X11R5 (mit/util/scripts/install.sh), which was
# later released in X11R6 (xc/config/util/install.sh) with the # later released in X11R6 (xc/config/util/install.sh) with the
...@@ -271,15 +271,18 @@ do ...@@ -271,15 +271,18 @@ do
fi fi
dst=$dst_arg dst=$dst_arg
# If destination is a directory, append the input filename; won't work # If destination is a directory, append the input filename.
# if double slashes aren't ignored.
if test -d "$dst"; then if test -d "$dst"; then
if test "$is_target_a_directory" = never; then if test "$is_target_a_directory" = never; then
echo "$0: $dst_arg: Is a directory" >&2 echo "$0: $dst_arg: Is a directory" >&2
exit 1 exit 1
fi fi
dstdir=$dst dstdir=$dst
dst=$dstdir/`basename "$src"` dstbase=`basename "$src"`
case $dst in
*/) dst=$dst$dstbase;;
*) dst=$dst/$dstbase;;
esac
dstdir_status=0 dstdir_status=0
else else
dstdir=`dirname "$dst"` dstdir=`dirname "$dst"`
...@@ -288,6 +291,11 @@ do ...@@ -288,6 +291,11 @@ do
fi fi
fi fi
case $dstdir in
*/) dstdirslash=$dstdir;;
*) dstdirslash=$dstdir/;;
esac
obsolete_mkdir_used=false obsolete_mkdir_used=false
if test $dstdir_status != 0; then if test $dstdir_status != 0; then
...@@ -324,14 +332,16 @@ do ...@@ -324,14 +332,16 @@ do
# is incompatible with FreeBSD 'install' when (umask & 300) != 0. # is incompatible with FreeBSD 'install' when (umask & 300) != 0.
;; ;;
*) *)
# $RANDOM is not portable (e.g. dash); use it when possible to # Note that $RANDOM variable is not portable (e.g. dash); Use it
# lower collision chance # here however when possible just to lower collision chance.
tmpdir=${TMPDIR-/tmp}/ins$RANDOM-$$ tmpdir=${TMPDIR-/tmp}/ins$RANDOM-$$
trap 'ret=$?; rmdir "$tmpdir/a/b" "$tmpdir/a" "$tmpdir" 2>/dev/null; exit $ret' 0 trap 'ret=$?; rmdir "$tmpdir/a/b" "$tmpdir/a" "$tmpdir" 2>/dev/null; exit $ret' 0
# As "mkdir -p" follows symlinks and we work in /tmp possibly; so # Because "mkdir -p" follows existing symlinks and we likely work
# create the $tmpdir first (and fail if unsuccessful) to make sure # directly in world-writeable /tmp, make sure that the '$tmpdir'
# that nobody tries to guess the $tmpdir name. # directory is successfully created first before we actually test
# 'mkdir -p' feature.
if (umask $mkdir_umask && if (umask $mkdir_umask &&
$mkdirprog $mkdir_mode "$tmpdir" && $mkdirprog $mkdir_mode "$tmpdir" &&
exec $mkdirprog $mkdir_mode -p -- "$tmpdir/a/b") >/dev/null 2>&1 exec $mkdirprog $mkdir_mode -p -- "$tmpdir/a/b") >/dev/null 2>&1
...@@ -434,8 +444,8 @@ do ...@@ -434,8 +444,8 @@ do
else else
# Make a couple of temp file names in the proper directory. # Make a couple of temp file names in the proper directory.
dsttmp=$dstdir/_inst.$$_ dsttmp=${dstdirslash}_inst.$$_
rmtmp=$dstdir/_rm.$$_ rmtmp=${dstdirslash}_rm.$$_
# Trap to clean up those temp files at exit. # Trap to clean up those temp files at exit.
trap 'ret=$?; rm -f "$dsttmp" "$rmtmp" && exit $ret' 0 trap 'ret=$?; rm -f "$dsttmp" "$rmtmp" && exit $ret' 0
...@@ -500,9 +510,9 @@ do ...@@ -500,9 +510,9 @@ do
done done
# Local variables: # Local variables:
# eval: (add-hook 'write-file-hooks 'time-stamp) # eval: (add-hook 'before-save-hook 'time-stamp)
# time-stamp-start: "scriptversion=" # time-stamp-start: "scriptversion="
# time-stamp-format: "%:y-%02m-%02d.%02H" # time-stamp-format: "%:y-%02m-%02d.%02H"
# time-stamp-time-zone: "UTC" # time-stamp-time-zone: "UTC0"
# time-stamp-end: "; # UTC" # time-stamp-end: "; # UTC"
# End: # End:
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
* Copyright (C) 2005-2016 Max-Planck-Society
* \author Martin Reinecke
*/
#include <math.h>
#include <stdlib.h>
#include "trig_utils.h"
#include "fftpack.h"
#include "bluestein.h"
/* returns the sum of all prime factors of n */
size_t prime_factor_sum (size_t n)
{
size_t result=0,tmp;
while (((tmp=(n>>1))<<1)==n)
{ result+=2; n=tmp; }
size_t limit=(size_t)sqrt(n+0.01);
for (size_t x=3; x<=limit; x+=2)
while ((tmp=(n/x))*x==n)
{
result+=x;
n=tmp;
limit=(size_t)sqrt(n+0.01);
}
if (n>1) result+=n;
return result;
}
/* returns the smallest composite of 2, 3 and 5 which is >= n */
static size_t good_size(size_t n)
{
if (n<=6) return n;
size_t bestfac=2*n;
for (size_t f2=1; f2<bestfac; f2*=2)
for (size_t f23=f2; f23<bestfac; f23*=3)
for (size_t f235=f23; f235<bestfac; f235*=5)
if (f235>=n) bestfac=f235;
return bestfac;
}
void bluestein_i (size_t n, double **tstorage, size_t *worksize)
{
size_t n2=good_size(n*2-1);
*worksize=2+2*n+8*n2+16;
*tstorage = RALLOC(double,2+2*n+8*n2+16);
((size_t *)(*tstorage))[0]=n2;
double *bk = *tstorage+2;
double *bkf = *tstorage+2+2*n;
double *work= *tstorage+2+2*(n+n2);
/* initialize b_k */
double *tmp=RALLOC(double,4*n);
sincos_2pibyn(2*n,2*n,&tmp[1],&tmp[0],2);
bk[0] = 1;
bk[1] = 0;
size_t coeff=0;
for (size_t m=1; m<n; ++m)
{
coeff+=2*m-1;
if (coeff>=2*n) coeff-=2*n;
bk[2*m ] = tmp[2*coeff ];
bk[2*m+1] = tmp[2*coeff+1];
}
DEALLOC(tmp);
/* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
double xn2 = 1./n2;
bkf[0] = bk[0]*xn2;
bkf[1] = bk[1]*xn2;
for (size_t m=2; m<2*n; m+=2)
{
bkf[m] = bkf[2*n2-m] = bk[m] *xn2;
bkf[m+1] = bkf[2*n2-m+1] = bk[m+1] *xn2;
}
for (size_t m=2*n;m<=(2*n2-2*n+1);++m)
bkf[m]=0.;
cffti (n2,work);
cfftf (n2,bkf,work);
}
void bluestein (size_t n, double *data, double *tstorage, int isign)
{
size_t n2=*((size_t *)tstorage);
double *bk = tstorage+2;
double *bkf = tstorage+2+2*n;
double *work= tstorage+2+2*(n+n2);
double *akf = tstorage+2+2*n+6*n2+16;
/* initialize a_k and FFT it */
if (isign>0)
for (size_t m=0; m<2*n; m+=2)
{
akf[m] = data[m]*bk[m] - data[m+1]*bk[m+1];
akf[m+1] = data[m]*bk[m+1] + data[m+1]*bk[m];
}
else
for (size_t m=0; m<2*n; m+=2)
{
akf[m] = data[m]*bk[m] + data[m+1]*bk[m+1];
akf[m+1] =-data[m]*bk[m+1] + data[m+1]*bk[m];
}
for (size_t m=2*n; m<2*n2; ++m)
akf[m]=0;
cfftf (n2,akf,work);
/* do the convolution */
if (isign>0)
for (size_t m=0; m<2*n2; m+=2)
{
double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1];
akf[m+1] = im;
}
else
for (size_t m=0; m<2*n2; m+=2)
{
double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1];
akf[m+1] = im;
}
/* inverse FFT */
cfftb (n2,akf,work);
/* multiply by b_k* */
if (isign>0)
for (size_t m=0; m<2*n; m+=2)
{
data[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1];
data[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1];
}
else
for (size_t m=0; m<2*n; m+=2)
{
data[m] = bk[m] *akf[m] + bk[m+1]*akf[m+1];
data[m+1] =-bk[m+1]*akf[m] + bk[m] *akf[m+1];
}
}
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*! \file bluestein.h
* Interface for the Bluestein algorithm
*
* Copyright (C) 2005 Max-Planck-Society
* \author Martin Reinecke
*/
#ifndef PLANCK_BLUESTEIN_H
#define PLANCK_BLUESTEIN_H
#include "c_utils.h"
#ifdef __cplusplus
extern "C" {
#endif
size_t prime_factor_sum (size_t n);
void bluestein_i (size_t n, double **tstorage, size_t *worksize);
void bluestein (size_t n, double *data, double *tstorage, int isign);
#ifdef __cplusplus
}
#endif
#endif
This diff is collapsed.
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
fftpack.h : function declarations for fftpack.c
Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
(Version 4, 1985).
Pekka Janhunen 23.2.1995
(reformatted by joerg arndt)
reformatted and slightly enhanced by Martin Reinecke (2004)
*/
#ifndef PLANCK_FFTPACK_H
#define PLANCK_FFTPACK_H
#include "c_utils.h"
#ifdef __cplusplus
extern "C" {
#endif
/*! forward complex transform */
void cfftf(size_t N, double complex_data[], double wrk[]);
/*! backward complex transform */
void cfftb(size_t N, double complex_data[], double wrk[]);
/*! initializer for complex transforms */
void cffti(size_t N, double wrk[]);
/*! forward real transform */
void rfftf(size_t N, double data[], double wrk[]);
/*! backward real transform */
void rfftb(size_t N, double data[], double wrk[]);
/*! initializer for real transforms */
void rffti(size_t N, double wrk[]);
#ifdef __cplusplus
}
#endif
#endif
/*
* This file is part of libfftpack.
*
* libfftpack is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* libfftpack is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with libfftpack; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*
* libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
* and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
* (DLR).
*/
/*
fftpack.c : A set of FFT routines in C.
Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber
(Version 4, 1985).
C port by Martin Reinecke (2010)
*/
#ifdef BACKWARD
#define PSIGN +
#define PMSIGNC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; }
/* a = b*c */
#define MULPMSIGNC(a,b,c) { a.r=b.r*c.r-b.i*c.i; a.i=b.r*c.i+b.i*c.r; }
#else
#define PSIGN -
#define PMSIGNC(a,b,c,d) { a.r=c.r-d.r; a.i=c.i-d.i; b.r=c.r+d.r; b.i=c.i+d.i; }
/* a = conj(b)*c */
#define MULPMSIGNC(a,b,c) { a.r=b.r*c.r+b.i*c.i; a.i=b.r*c.i-b.i*c.r; }
#endif
static void X(2) (size_t ido, size_t l1, const cmplx *cc, cmplx *ch,
const cmplx *wa)
{
const size_t cdim=2;
if (ido==1) // CC(1,2,l1) CH(1,l1,2)
for (size_t k=0;k<l1;++k)
PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k))
else // CC(ido,2,l1) CH(ido,l1,2) WA(1,ido)
for (size_t k=0;k<l1;++k)
for (size_t i=0;i<ido;++i)
{
cmplx t;
PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k))
MULPMSIGNC (CH(i,k,1),WA(0,i),t)