[
Date Prev][
Date Next][
Thread Prev][
Thread Next][
Date Index][
Thread Index]
[patch] Memory allocation cleanup, disable non power-of-2 FFT tests when using SAL.
- To: VSIPL++ Developers List <vsipl++@xxxxxxxxxxxxxxxx>
- Subject: [patch] Memory allocation cleanup, disable non power-of-2 FFT tests when using SAL.
- From: Jules Bergmann <jules@xxxxxxxxxxxxxxxx>
- Date: Fri, 03 Mar 2006 13:14:41 -0500
This patch
* Closes a memory leak (Fft objects were creating a reference counted
Fft_core object, but the initial reference count was 2, making it
impossible for the count to go to 0).
* Changes the FFT/FFTM/FFTM-par tests to disable non-power-of-2
FFT sizes when using SAL, and other cleanup.
* Changes the window test to disable non-power-of-2 Chebychev
test when using SAL.
-- Jules
--
Jules Bergmann
CodeSourcery
jules@xxxxxxxxxxxxxxxx
(650) 331-3385 x705
Index: ChangeLog
===================================================================
RCS file: /home/cvs/Repository/vpp/ChangeLog,v
retrieving revision 1.401
diff -u -r1.401 ChangeLog
--- ChangeLog 3 Mar 2006 14:30:31 -0000 1.401
+++ ChangeLog 3 Mar 2006 18:05:38 -0000
@@ -1,5 +1,25 @@
2006-03-03 Jules Bergmann <jules@xxxxxxxxxxxxxxxx>
+ * configure.ac (--with-mpi-prefix64): New option, similar to
+ --with-mpi-prefix, that searches include64, bin64, and lib64
+ instead. Useful for 64-bit Intel MPI library.
+ * src/vsip/impl/signal-fft.hpp: Allocate ref-counted Fft_core
+ with 'noincrement' to prevent memory leak.
+ * src/vsip/impl/sal/fft.hpp: Fix Wall signed/unsigned warning.
+ * tests/fft.cpp: Cleanup, move reference dft into ref_dft.hpp.
+ Use error_db from error_db.hpp.
+ * tests/fftm-par.cpp: Cleanup, move reference dft into ref_dft.hpp,
+ use error_db from error_db.hpp, disable non-power-of-2 FFTs
+ when using SAL.
+ * tests/fftm.cpp: Likewise.
+ * tests/random.cpp: Deallocate vsip_randstate.
+ * tests/ref_dft.hpp: New file, reference implementation of
+ DFT and DFTM.
+ * tests/window.cpp: Disable test of non-power-of-2 Chebychev
+ window when using SAL FFT.
+
+2006-03-03 Jules Bergmann <jules@xxxxxxxxxxxxxxxx>
+
* autogen.sh: Doctor generated configure file to ignore extra
files generated by GreenHills compiler when detecting object
and executable exentions.
Index: configure.ac
===================================================================
RCS file: /home/cvs/Repository/vpp/configure.ac,v
retrieving revision 1.83
diff -u -r1.83 configure.ac
--- configure.ac 3 Mar 2006 14:30:31 -0000 1.83
+++ configure.ac 3 Mar 2006 18:05:38 -0000
@@ -71,6 +71,12 @@
must be in PATH/include; libraries in PATH/lib.]),
dnl If the user specified --with-mpi-prefix, they mean to use MPI for sure.
[enable_mpi=yes])
+AC_ARG_WITH(mpi_prefix64,
+ AS_HELP_STRING([--with-mpi-prefix64=PATH],
+ [Specify the installation prefix of the MPI library. Headers
+ must be in PATH/include64; libraries in PATH/lib64.]),
+ dnl If the user specified --with-mpi-prefix64, they mean to use MPI for sure.
+ [enable_mpi=yes])
### Mercury Scientific Algorithm (SAL)
AC_ARG_ENABLE([sal],
@@ -732,9 +738,12 @@
if test "$enable_mpi" != "no"; then
vsipl_par_service=1
- if test -n "$with_mpi_prefix"; then
+ if test -n "$with_mpi_prefix" -a "$with_mpi_prefix" != "/usr"; then
MPI_CPPFLAGS="-I$with_mpi_prefix/include"
MPI_LIBS="-L$with_mpi_prefix/lib -L$with_mpi_prefix/lib64"
+ elif test -n "$with_mpi_prefix64"; then
+ MPI_CPPFLAGS="-I$with_mpi_prefix64/include64"
+ MPI_LIBS="-L$with_mpi_prefix64/lib64"
fi
save_CPPFLAGS="$CPPFLAGS"
CPPFLAGS="$CPPFLAGS $MPI_CPPFLAGS"
@@ -785,9 +794,11 @@
mpich)
check_mpicxx="yes"
if test -n "$with_mpi_prefix"; then
- MPICXX="$with_mpi_prefix/bin/mpicxx -show"
+ MPICXX="$with_mpi_prefix/bin/mpicxx -nocompchk -show"
+ elif test -n "$with_mpi_prefix64"; then
+ MPICXX="$with_mpi_prefix64/bin64/mpicxx -nocompchk -show"
else
- MPICXX="mpicxx -show"
+ MPICXX="mpicxx -nocompchk -show"
fi
;;
Index: src/vsip/impl/signal-fft.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/signal-fft.hpp,v
retrieving revision 1.32
diff -u -r1.32 signal-fft.hpp
--- src/vsip/impl/signal-fft.hpp 1 Mar 2006 16:36:47 -0000 1.32
+++ src/vsip/impl/signal-fft.hpp 3 Mar 2006 18:05:39 -0000
@@ -330,7 +330,7 @@
Fft_imp(const vsip::Domain<dim>& dom, impl_scalar_type scale)
VSIP_THROW((std::bad_alloc))
- : core_(new impl_core_type)
+ : core_(new impl_core_type, noincrement)
, input_size_( Impl_arg_size<inT,outT>::size(dom))
, output_size_(Impl_arg_size<outT,inT>::size(dom))
, scale_(scale)
Index: src/vsip/impl/sal/fft.hpp
===================================================================
RCS file: /home/cvs/Repository/vpp/src/vsip/impl/sal/fft.hpp,v
retrieving revision 1.4
diff -u -r1.4 fft.hpp
--- src/vsip/impl/sal/fft.hpp 1 Mar 2006 16:36:47 -0000 1.4
+++ src/vsip/impl/sal/fft.hpp 3 Mar 2006 18:05:39 -0000
@@ -610,7 +610,7 @@
// Unpack the data (see SAL reference for details), and scale back by 1/2.
int const N = (1 << self.size_[1]) + 2;
float scale = 0.5f;
- for (unsigned int i = 0; i != self.runs_; ++i, out += self.dist_ + 2)
+ for (int i = 0; i != self.runs_; ++i, out += self.dist_ + 2)
{
out[(N - 2) * self.stride_] = out[self.stride_];
out[self.stride_] = 0.f;
@@ -633,7 +633,7 @@
// Unpack the data (see SAL reference for details), and scale back by 1/2.
int const N = (1 << self.size_[1]) + 2;
double scale = 0.5f;
- for (unsigned int i = 0; i != self.runs_; ++i, out += self.dist_ + 2)
+ for (int i = 0; i != self.runs_; ++i, out += self.dist_ + 2)
{
out[(N - 2) * self.stride_] = out[self.stride_];
out[self.stride_] = 0.f;
Index: tests/fft.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/fft.cpp,v
retrieving revision 1.13
diff -u -r1.13 fft.cpp
--- tests/fft.cpp 2 Mar 2006 20:07:20 -0000 1.13
+++ tests/fft.cpp 3 Mar 2006 18:05:40 -0000
@@ -21,10 +21,8 @@
#include "test.hpp"
#include "output.hpp"
-
-#ifndef HAVE_SINL
-# define HAVE_SINL 0
-#endif
+#include "error_db.hpp"
+#include "ref_dft.hpp"
@@ -37,152 +35,6 @@
-/// Return sin and cos of phi as complex.
-
-#if HAVE_SINL
-template <typename T>
-vsip::complex<T>
-SinCos(long double phi)
-{
- return vsip::complex<T>(cos(phi), sin(phi));
-}
-#endif
-
-template <typename T>
-vsip::complex<T>
-SinCos(double phi)
-{
- return vsip::complex<T>(cos(phi), sin(phi));
-}
-
-
-
-// Reference 1-D DFT algorithm. Brutes it out, but easy to validate
-// and works for any size.
-
-// Requires:
-// IN to be input Vector.
-// OUT to be output Vector, of same size as IN.
-// IDIR to be sign of exponential.
-// -1 => Forward Fft,
-// +1 => Inverse Fft.
-
-template <typename T,
- typename Block1,
- typename Block2>
-void dft(
- vsip::Vector<vsip::complex<T>, Block1> in,
- vsip::Vector<vsip::complex<T>, Block2> out,
- int idir)
-{
- length_type const size = in.size();
- test_assert(in.size() == out.size());
-#if HAVE_SINL
- typedef long double AccT;
-#else
- typedef double AccT;
-#endif
-
- AccT const phi = idir * 2.0 * M_PI/size;
-
- for (index_type w=0; w<size; ++w)
- {
- vsip::complex<AccT> sum = vsip::complex<AccT>();
- for (index_type k=0; k<size; ++k)
- sum += vsip::complex<AccT>(in(k)) * SinCos<AccT>(phi*k*w);
- out(w) = vsip::complex<T>(sum);
- }
-}
-
-
-
-// Error metric between two vectors.
-
-template <typename T1,
- typename T2,
- typename Block1,
- typename Block2>
-double
-error_db(
- const_Vector<T1, Block1> v1,
- const_Vector<T2, Block2> v2)
-{
- double refmax = 0.0;
- double maxsum = -250;
- double sum;
-
- for (index_type i=0; i<v1.size(); ++i)
- {
- double val = magsq(v1.get(i));
- if (val > refmax)
- refmax = val;
- }
-
-
- for (index_type i=0; i<v1.size(); ++i)
- {
- double val = magsq(v1.get(i) - v2.get(i));
-
- if (val < 1.e-20)
- sum = -201.;
- else
- sum = 10.0 * log10(val/(2.0*refmax));
-
- if (sum > maxsum)
- maxsum = sum;
- }
-
- return maxsum;
-}
-
-// Error metric between two Matrices.
-
-template <typename T1,
- typename T2,
- typename Block1,
- typename Block2>
-double
-error_db(
- const_Matrix<T1, Block1> v1,
- const_Matrix<T2, Block2> v2)
-{
- double maxsum = -250;
- for (unsigned i = 0; i < v1.size(0); ++i)
- {
- double sum = error_db(v1.row(i), v2.row(i));
- if (sum > maxsum)
- maxsum = sum;
- }
- return maxsum;
-}
-
-
-
-// Error metric between two Tensors.
-
-template <typename T1,
- typename T2,
- typename Block1,
- typename Block2>
-double
-error_db(
- const_Tensor<T1, Block1> v1,
- const_Tensor<T2, Block2> v2)
-{
- double maxsum = -250;
- for (unsigned i = 0; i < v1.size(0); ++i)
- {
- vsip::Domain<1> y(v1.size(1));
- vsip::Domain<1> x(v1.size(2));
- double sum = error_db(v1(i,y,x), v2(i,y,x));
- if (sum > maxsum)
- maxsum = sum;
- }
- return maxsum;
-}
-
-
-
// Setup input data for Fft.
template <typename T,
@@ -205,11 +57,11 @@
case 2:
in = T();
in(0) = T(1);
- in(Domain<1>(0, 1, N)) += T(3);
- in(Domain<1>(0, 4, N/4)) += T(-2);
- in(Domain<1>(0, 13, N/13)) += T(7);
- in(Domain<1>(0, 27, N/27)) += T(-15);
- in(Domain<1>(0, 37, N/37)) += T(31);
+ if (N > 1) in(Domain<1>(0, 1, N)) += T(3);
+ if (N > 4) in(Domain<1>(0, 4, N/4)) += T(-2);
+ if (N > 13) in(Domain<1>(0, 13, N/13)) += T(7);
+ if (N > 27) in(Domain<1>(0, 27, N/27)) += T(-15);
+ if (N > 37) in(Domain<1>(0, 37, N/37)) += T(31);
break;
case 3:
in = T(scale);
@@ -251,7 +103,7 @@
setup_data(set, in);
- dft(in, ref, -1);
+ ref::dft(in, ref, -1);
f_fft(in, out);
i_fft(out, inv);
@@ -299,7 +151,7 @@
setup_data(set, in);
- dft(in, ref, -1);
+ ref::dft(in, ref, -1);
out = f_fft(in);
inv = i_fft(out);
Index: tests/fftm-par.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/fftm-par.cpp,v
retrieving revision 1.5
diff -u -r1.5 fftm-par.cpp
--- tests/fftm-par.cpp 20 Dec 2005 12:48:40 -0000 1.5
+++ tests/fftm-par.cpp 3 Mar 2006 18:05:40 -0000
@@ -26,6 +26,19 @@
#include "output.hpp"
#include "util.hpp"
#include "util-par.hpp"
+#include "error_db.hpp"
+#include "ref_dft.hpp"
+
+#define VERBOSE 0
+
+#if VSIP_IMPL_SAL_FFT
+# define TEST_NON_REALCOMPLEX 0
+# define TEST_NON_POWER_OF_2 0
+#else
+# define TEST_NON_REALCOMPLEX 1
+# define TEST_NON_POWER_OF_2 1
+#endif
+
/***********************************************************************
@@ -38,53 +51,6 @@
int number_of_processors;
-/// Return sin and cos of phi as complex.
-
-template <typename T>
-vsip::complex<T>
-sin_cos(T phi)
-{
- return vsip::complex<T>(std::cos(phi), std::sin(phi));
-}
-
-
-template <typename T,
- typename Block1,
- typename Block2>
-void dft_x(
- vsip::Matrix<vsip::complex<T>, Block1> in,
- vsip::Matrix<vsip::complex<T>, Block2> out)
-{
- length_type const xsize = in.size(1);
- test_assert(in.size(0) == out.size(0));
- test_assert(in.size(1) == out.size(1));
-
- typedef vsip::complex<T> CT;
-
- Fft<const_Vector,CT,CT,fft_fwd,by_reference,1>
- fft(Domain<1>(xsize), 1.0);
-
- Block1& in_block = in.block();
- Matrix<CT,Dense<2,CT,tuple<0,1,2>,Local_map> > local_in(
- in_block.get_local_block());
-
- Block2& out_block = out.block();
- Matrix<CT,Dense<2,CT,tuple<0,1,2>,Local_map> > local_out(
- out_block.get_local_block());
-
- typedef const_Matrix<CT> in_mat_type;
- typedef Matrix<CT> out_mat_type;
- const in_mat_type const_in(local_in);
-
- length_type const lysize = const_in.size(0);
-
- for (index_type v=0; v < lysize; ++v)
- {
- const typename in_mat_type::const_row_type inrow = const_in.row(v);
- typename out_mat_type::row_type outrow = local_out.row(v);
- fft.operator()(inrow, outrow);
- }
-}
template <typename Block>
void
@@ -124,136 +90,6 @@
}
-template <typename T,
- typename Block1,
- typename Block2>
-void dft_y(
- vsip::Matrix<vsip::complex<T>, Block1> in,
- vsip::Matrix<vsip::complex<T>, Block2> out)
-{
- length_type const xsize = in.size(1);
- length_type const ysize = in.size(0);
- test_assert(in.size(0) == out.size(0));
- test_assert(in.size(1) == out.size(1));
-
- typedef vsip::complex<T> CT;
-
- Fft<const_Vector,CT,CT,fft_fwd,by_reference,1>
- fft(Domain<1>(ysize), 1.0);
-
- Block1& in_block = in.block();
- Matrix<CT,Dense<2,CT,tuple<0,1,2>,Local_map> > local_in(
- in_block.get_local_block());
-
- Block2& out_block = out.block();
- Matrix<CT,Dense<2,CT,tuple<0,1,2>,Local_map> > local_out(
- out_block.get_local_block());
-
- typedef const_Matrix<CT> in_mat_type;
- typedef Matrix<CT> out_mat_type;
- const in_mat_type const_in(local_in);
-
- for (index_type v=0; v < xsize; ++v)
- {
- if (in_block.subblock() == in_block.map().impl_subblock_from_index(1, v))
- {
- unsigned ix = in_block.map().impl_local_from_global_index(1, v);
- const typename in_mat_type::const_col_type incol = const_in.col(ix);
- typename out_mat_type::col_type outcol = local_out.col(ix);
- fft.operator()(incol, outcol);
- }
- }
-}
-
-
-// Error metric between two vectors.
-
-template <typename T,
- typename Block1,
- typename Block2>
-double
-error_db(
- const_Matrix<T, Block1> v1,
- const_Matrix<T, Block2> v2)
-{
- Matrix<T> lv1(v1.block().get_local_block());
- Matrix<T> lv2(v2.block().get_local_block());
-
- double refmax = 0.0;
- double maxsum = -250;
- double sum;
-
- for (index_type j=0; j < lv1.size(0); ++j)
- for (index_type i=0; i < lv1.size(1); ++i)
- {
- double tmp = std::abs(lv1.get(j,i));
- double val = tmp*tmp;
- if (val > refmax)
- refmax = val;
- }
-
- impl::Communicator comm = impl::default_communicator();
- int rank = comm.rank();
- int size = comm.size();
-
- if (rank != 0)
- {
- comm.buf_send(0, &refmax, 1);
- comm.recv(0, &refmax, 1);
- }
- else
- {
- for (int i = 1; i < size; ++i)
- {
- double otherefmax;
- comm.recv(i, &otherefmax, 1);
- if (refmax < otherefmax)
- refmax = otherefmax;
- }
- for (int i = 1; i < size; ++i)
- comm.buf_send(i, &refmax, 1);
- }
-
-
- for (index_type j=0; j < lv1.size(0); ++j)
- for (index_type i=0; i < lv1.size(1); ++i)
- {
- double tmp = std::abs(lv1.get(j,i) - lv2.get(j,i));
- double val = tmp*tmp;
-
- if (val < 1.e-20)
- sum = -201.;
- else
- sum = 10.0 * log10(val/(2.0*refmax));
-
- if (sum > maxsum)
- maxsum = sum;
- }
-
- if (rank != 0)
- {
- comm.buf_send(0, &maxsum, 1);
- comm.recv(0, &maxsum, 1);
- }
- else
- {
- for (int i = 1; i < size; ++i)
- {
- double othersum;
- comm.recv(i, &othersum, 1);
- if (maxsum < othersum)
- maxsum = othersum;
- }
- for (int i = 1; i < size; ++i)
- comm.buf_send(i, &maxsum, 1);
- return maxsum;
- }
-
- return -200; // only root can test_assert failure.
-
-}
-
-
// Set up input data for Fftm.
@@ -355,7 +191,7 @@
dist_matrix_type inv(inv_block);
setup_data_x(in);
- dft_x(in, ref);
+ ref::dft_x(in, ref, -1);
f_fftm(in, out);
i_fftm(out, inv);
@@ -413,7 +249,7 @@
setup_data_x(in);
- dft_x(in, ref);
+ ref::dft_x(in, ref, -1);
out = f_fftm(in);
inv = i_fftm(out);
@@ -564,28 +400,28 @@
dist_matrix_type inv(inv_block);
setup_data_y(in);
- dft_y(in, ref);
+ ref::dft_y(in, ref, -1);
-#if 0
+#if VERBOSE
cout.precision(3);
cout.setf(ios_base::fixed);
#endif
-#if 0
+#if VERBOSE
dump_matrix(in.block(), N, 1);
dump_matrix(ref.block(), N, 1);
#endif
f_fftm(in, out);
-#if 0
+#if VERBOSE
dump_matrix(in.block(), N, 1);
dump_matrix(out.block(), N, 1);
#endif
i_fftm(out, inv);
-#if 0
+#if VERBOSE
dump_matrix(out.block(), N, 1);
dump_matrix(inv.block(), N, 1);
#endif
@@ -641,7 +477,7 @@
setup_data_y(in);
- dft_y(in, ref);
+ ref::dft_y(in, ref, -1);
out = f_fftm(in);
inv = i_fftm(out);
@@ -709,6 +545,47 @@
#endif
+
+
+template <typename T>
+void
+test()
+{
+ // Test powers-of-2
+ test_by_ref_x<complex<T> >(64);
+ test_by_ref_x<complex<T> >(256);
+
+ test_by_ref_y<complex<T> >(256);
+
+ test_by_val_x<complex<T> >(128);
+ test_by_val_x<complex<T> >(256);
+ test_by_val_x<complex<T> >(512);
+
+ test_by_val_y<complex<T> >(256);
+
+
+#if TEST_NON_REALCOMPLEX
+ test_real<T>(128);
+ test_real<T>(16);
+#endif
+
+#if TEST_NON_POWER_OF_2
+ // Test non-powers-of-2
+ test_by_ref_x<complex<T> >(18);
+ test_by_ref_x<complex<T> >(68);
+ test_by_ref_x<complex<T> >(252);
+
+ test_by_ref_y<complex<T> >(68);
+
+ test_by_val_y<complex<T> >(18);
+
+ // Tests for test r->c, c->r.
+# if TEST_NON_REALCOMPLEX
+ test_real<T>(242);
+# endif
+#endif
+};
+
int
main(int argc, char** argv)
{
@@ -734,53 +611,11 @@
#endif
#if defined(VSIP_IMPL_FFT_USE_FLOAT)
- test_by_ref_x<complex<float> >(18);
- test_by_ref_x<complex<float> >(64);
- test_by_ref_x<complex<float> >(68);
- test_by_ref_x<complex<float> >(256);
- test_by_ref_x<complex<float> >(252);
-
- test_by_ref_y<complex<float> >(68);
- test_by_ref_y<complex<float> >(256);
-
- test_by_val_x<complex<float> >(128);
- test_by_val_x<complex<float> >(256);
- test_by_val_x<complex<float> >(512);
-
- test_by_val_y<complex<float> >(18);
- test_by_val_y<complex<float> >(256);
-
-# if 0
- // Tests for test r->c, c->r.
- test_real<float>(128);
- test_real<float>(242);
- test_real<float>(16);
-# endif
+ test<float>();
#endif
#if defined(VSIP_IMPL_FFT_USE_DOUBLE)
- test_by_ref_x<complex<double> >(18);
- test_by_ref_x<complex<double> >(64);
- test_by_ref_x<complex<double> >(68);
- test_by_ref_x<complex<double> >(256);
- test_by_ref_x<complex<double> >(252);
-
- test_by_ref_y<complex<double> >(68);
- test_by_ref_y<complex<double> >(256);
-
- test_by_val_x<complex<double> >(128);
- test_by_val_x<complex<double> >(256);
- test_by_val_x<complex<double> >(512);
-
- test_by_val_y<complex<double> >(18);
- test_by_val_y<complex<double> >(256);
-
-# if 0
- // Tests for test r->c, c->r.
- test_real<double>(128);
- test_real<double>(242);
- test_real<double>(16);
-# endif
+ test<double>();
#endif
return 0;
Index: tests/fftm.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/fftm.cpp,v
retrieving revision 1.10
diff -u -r1.10 fftm.cpp
--- tests/fftm.cpp 22 Dec 2005 08:21:23 -0000 1.10
+++ tests/fftm.cpp 3 Mar 2006 18:05:40 -0000
@@ -22,6 +22,17 @@
#include "test.hpp"
#include "output.hpp"
+#include "error_db.hpp"
+#include "ref_dft.hpp"
+
+#if VSIP_IMPL_SAL_FFT
+# define TEST_NON_REALCOMPLEX 0
+# define TEST_NON_POWER_OF_2 0
+#else
+# define TEST_NON_REALCOMPLEX 1
+# define TEST_NON_POWER_OF_2 1
+#endif
+
/***********************************************************************
@@ -33,148 +44,6 @@
-/// Return sin and cos of phi as complex.
-
-template <typename T>
-vsip::complex<T>
-sin_cos(T phi)
-{
- return vsip::complex<T>(std::cos(phi), std::sin(phi));
-}
-
-
-// Reference 1-D DFT algorithm. Brutes it out, but easy to validate
-// and works for any size.
-
-// Requires:
-// IN to be input Matrix.
-// OUT to be output Matrix, of same size as IN.
-// IDIR to be sign of exponential.
-// -1 => Forward Fft,
-// +1 => Inverse Fft.
-
-template <typename T,
- typename Block1,
- typename Block2>
-void dft_x(
- vsip::Matrix<vsip::complex<T>, Block1> in,
- vsip::Matrix<vsip::complex<T>, Block2> out,
- int idir)
-{
- length_type const xsize = in.size(1);
- length_type const ysize = in.size(0);
- test_assert(in.size(0) == out.size(0));
- test_assert(in.size(1) == out.size(1));
- typedef long double AccT;
-
- AccT const phi = idir * 2.0 * M_PI/xsize;
-
- for (index_type v=0; v < ysize; ++v)
- for (index_type w=0; w < xsize; ++w)
- {
- vsip::complex<AccT> sum = vsip::complex<AccT>();
- for (index_type k=0; k<xsize; ++k)
- sum += vsip::complex<AccT>(in(v,k)) * sin_cos<AccT>(phi*k*w);
- out(v,w) = vsip::complex<T>(sum);
- }
-}
-
-
-template <typename T,
- typename Block1,
- typename Block2>
-void dft_y(
- vsip::Matrix<vsip::complex<T>, Block1> in,
- vsip::Matrix<vsip::complex<T>, Block2> out,
- int idir)
-{
- length_type const xsize = in.size(1);
- length_type const ysize = in.size(0);
- test_assert(in.size(0) == out.size(0));
- test_assert(in.size(1) == out.size(1));
- typedef long double AccT;
-
- AccT const phi = idir * 2.0 * M_PI/ysize;
-
- for (index_type v=0; v < xsize; ++v)
- for (index_type w=0; w < ysize; ++w)
- {
- vsip::complex<AccT> sum = vsip::complex<AccT>();
- for (index_type k=0; k<ysize; ++k)
- sum += vsip::complex<AccT>(in(k,v)) * sin_cos<AccT>(phi*k*w);
- out(w,v) = vsip::complex<T>(sum);
- }
-}
-
-
-template <typename T,
- typename Block1,
- typename Block2>
-void dft_y_real(
- vsip::Matrix<T, Block1> in,
- vsip::Matrix<vsip::complex<T>, Block2> out)
-{
- length_type const xsize = in.size(1);
- length_type const ysize = in.size(0);
- test_assert(in.size(0)/2 + 1 == out.size(0));
- test_assert(in.size(1) == out.size(1));
- typedef long double AccT;
-
- AccT const phi = -2.0 * M_PI/ysize;
-
- for (index_type v=0; v < xsize; ++v)
- for (index_type w=0; w < ysize / 2 + 1; ++w)
- {
- vsip::complex<AccT> sum = vsip::complex<AccT>();
- for (index_type k=0; k<ysize; ++k)
- sum += vsip::complex<AccT>(in(k,v)) * sin_cos<AccT>(phi*k*w);
- out(w,v) = vsip::complex<T>(sum);
- }
-}
-
-
-// Error metric between two vectors.
-
-template <typename T1,
- typename T2,
- typename Block1,
- typename Block2>
-double
-error_db(
- const_Matrix<T1, Block1> v1,
- const_Matrix<T2, Block2> v2)
-{
- double refmax = 0.0;
- double maxsum = -250;
- double sum;
-
- for (index_type j=0; j < v1.size(0); ++j)
- for (index_type i=0; i < v1.size(1); ++i)
- {
- double val = magsq(v1.get(j,i));
- if (val > refmax)
- refmax = val;
- }
-
- for (index_type j=0; j < v1.size(0); ++j)
- for (index_type i=0; i < v1.size(1); ++i)
- {
- double val = magsq(v1.get(j,i) - v2.get(j,i));
-
- if (val < 1.e-20)
- sum = -201.;
- else
- sum = 10.0 * log10(val/(2.0*refmax));
-
- if (sum > maxsum)
- maxsum = sum;
- }
-
- return maxsum;
-}
-
-
-
// Set up input data for Fftm.
template <typename T,
@@ -240,7 +109,7 @@
setup_data_x(in);
- dft_x(in, ref, -1);
+ ref::dft_x(in, ref, -1);
f_fftm(in, out);
@@ -294,7 +163,7 @@
setup_data_x(in);
- dft_x(in, ref, -1);
+ ref::dft_x(in, ref, -1);
out = f_fftm(in);
inv = i_fftm(out);
@@ -364,7 +233,7 @@
setup_data_y(in);
- dft_y(in, ref, -1);
+ ref::dft_y(in, ref, -1);
#if 0
cout.precision(3);
@@ -429,7 +298,7 @@
Matrix<T> inv(N, 5);
setup_data_y(in);
- dft_y(in, ref, -1);
+ ref::dft_y(in, ref, -1);
out = f_fftm(in);
inv = i_fftm(out);
@@ -471,7 +340,7 @@
Matrix<T> inv(N, 5);
setup_data_y(in);
- dft_y_real(in, ref);
+ ref::dft_y_real(in, ref);
out = f_fftm(in);
inv = i_fftm(out);
@@ -480,83 +349,62 @@
}
-
-int
-main()
+template <typename T>
+void
+test()
{
- vsipl init;
+ // Test powers-of-2
+ test_by_ref_x<complex<T> >(64);
+ test_by_ref_x<complex<T> >(256);
-#if defined(VSIP_IMPL_FFT_USE_FLOAT)
- test_by_ref_x<complex<float> >(18);
- test_by_ref_x<complex<float> >(64);
- test_by_ref_x<complex<float> >(68);
- test_by_ref_x<complex<float> >(256);
- test_by_ref_x<complex<float> >(252);
-
- test_by_ref_y<complex<float> >(68);
- test_by_ref_y<complex<float> >(256);
-
- test_by_val_x<complex<float> >(128);
- test_by_val_x<complex<float> >(256);
- test_by_val_x<complex<float> >(512);
+ test_by_ref_y<complex<T> >(256);
- test_by_val_y<complex<float> >(18);
- test_by_val_y<complex<float> >(256);
+ test_by_val_x<complex<T> >(128);
+ test_by_val_x<complex<T> >(256);
+ test_by_val_x<complex<T> >(512);
- // Tests for test r->c, c->r.
- test_real<float>(128);
- test_real<float>(242);
- test_real<float>(16);
+ test_by_val_y<complex<T> >(256);
+
+
+#if TEST_NON_REALCOMPLEX
+ test_real<T>(128);
+ test_real<T>(16);
#endif
-#if defined(VSIP_IMPL_FFT_USE_DOUBLE)
- test_by_ref_x<complex<double> >(18);
- test_by_ref_x<complex<double> >(64);
- test_by_ref_x<complex<double> >(68);
- test_by_ref_x<complex<double> >(256);
- test_by_ref_x<complex<double> >(252);
-
- test_by_ref_y<complex<double> >(68);
- test_by_ref_y<complex<double> >(256);
-
- test_by_val_x<complex<double> >(128);
- test_by_val_x<complex<double> >(256);
- test_by_val_x<complex<double> >(512);
+#if TEST_NON_POWER_OF_2
+ // Test non-powers-of-2
+ test_by_ref_x<complex<T> >(18);
+ test_by_ref_x<complex<T> >(68);
+ test_by_ref_x<complex<T> >(252);
+
+ test_by_ref_y<complex<T> >(68);
- test_by_val_y<complex<double> >(18);
- test_by_val_y<complex<double> >(256);
+ test_by_val_y<complex<T> >(18);
// Tests for test r->c, c->r.
- test_real<double>(128);
- test_real<double>(242);
- test_real<double>(16);
+ test_real<T>(242);
#endif
+};
+
-#if defined(VSIP_IMPL_FFT_USE_LONG_DOUBLE)
-#if ! defined(VSIP_IMPL_IPP_FFT)
- test_by_ref_x<complex<long double> >(18);
- test_by_ref_x<complex<long double> >(64);
- test_by_ref_x<complex<long double> >(68);
- test_by_ref_x<complex<long double> >(256);
- test_by_ref_x<complex<long double> >(252);
-
- test_by_ref_y<complex<long double> >(68);
- test_by_ref_y<complex<long double> >(256);
-
- test_by_val_x<complex<long double> >(128);
- test_by_val_x<complex<long double> >(256);
- test_by_val_x<complex<long double> >(512);
+int
+main()
+{
+ vsipl init;
- test_by_val_y<complex<long double> >(18);
- test_by_val_y<complex<long double> >(256);
+#if defined(VSIP_IMPL_FFT_USE_FLOAT)
+ test<float>();
+#endif
- // Tests for test r->c, c->r.
- test_real<long double>(128);
- test_real<long double>(242);
- test_real<long double>(16);
-#endif /* VSIP_IMPL_IPP_FFT */
+#if defined(VSIP_IMPL_FFT_USE_DOUBLE)
+ test<double>();
+#endif
+#if defined(VSIP_IMPL_FFT_USE_LONG_DOUBLE)
+# if ! defined(VSIP_IMPL_IPP_FFT)
+ test<long double>();
+# endif /* VSIP_IMPL_IPP_FFT */
#endif
return 0;
Index: tests/random.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/random.cpp,v
retrieving revision 1.2
diff -u -r1.2 random.cpp
--- tests/random.cpp 20 Dec 2005 12:48:41 -0000 1.2
+++ tests/random.cpp 3 Mar 2006 18:05:40 -0000
@@ -139,6 +139,14 @@
}
+int
+vsip_randdestroy(
+ vsip_randstate *state)
+{
+ if (state != NULL) free(state);
+ return 0;
+}
+
vsip_scalar_d vsip_randu_d(
vsip_randstate *state)
@@ -360,6 +368,7 @@
double b = vsip_randn_d( rstate );
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Normal distribution, non-portable
@@ -373,6 +382,7 @@
double b = vsip_randn_d( rstate );
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Uniform distribution, portable
@@ -386,6 +396,7 @@
double b = vsip_randu_d( rstate );
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Uniform distribution, non-portable
@@ -399,6 +410,7 @@
double b = vsip_randu_d( rstate );
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
@@ -416,6 +428,7 @@
complex<double> b(z.r, z.i);
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Normal distribution, non-portable
@@ -430,6 +443,7 @@
complex<double> b(z.r, z.i);
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Uniform distribution, portable
@@ -444,6 +458,7 @@
complex<double> b(z.r, z.i);
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
// Uniform distribution, non-portable
@@ -458,6 +473,7 @@
complex<double> b(z.r, z.i);
test_assert( equal( a, b ) );
}
+ vsip_randdestroy( rstate );
}
Index: tests/ref_dft.hpp
===================================================================
RCS file: tests/ref_dft.hpp
diff -N tests/ref_dft.hpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ tests/ref_dft.hpp 3 Mar 2006 18:05:40 -0000
@@ -0,0 +1,147 @@
+/* Copyright (c) 2006 by CodeSourcery, LLC. All rights reserved. */
+
+/** @file ref_dft.cpp
+ @author Jules Bergmann
+ @date 2006-03-03
+ @brief VSIPL++ Library: Reference implementation of DFT.
+*/
+
+#ifndef VSIP_REF_DFT_HPP
+#define VSIP_REF_DFT_HPP
+
+/***********************************************************************
+ Included Files
+***********************************************************************/
+
+#include <cassert>
+
+#include <vsip/complex.hpp>
+#include <vsip/vector.hpp>
+#include <vsip/math.hpp>
+
+
+
+/***********************************************************************
+ Definitions
+***********************************************************************/
+
+namespace ref
+{
+
+/// Return sin and cos of phi as complex.
+
+template <typename T>
+inline vsip::complex<T>
+sin_cos(double phi)
+{
+ return vsip::complex<T>(cos(phi), sin(phi));
+}
+
+
+
+// Reference 1-D DFT algorithm. Brutes it out, but easy to validate
+// and works for any size.
+
+// Requires:
+// IN to be input Vector.
+// OUT to be output Vector, of same size as IN.
+// IDIR to be sign of exponential.
+// -1 => Forward Fft,
+// +1 => Inverse Fft.
+
+template <typename T1,
+ typename T2,
+ typename Block1,
+ typename Block2>
+void dft(
+ vsip::Vector<T1, Block1> in,
+ vsip::Vector<T2, Block2> out,
+ int idir)
+{
+ using vsip::length_type;
+ using vsip::index_type;
+
+ length_type const size = in.size();
+ assert(in.size() == out.size());
+ typedef double AccT;
+
+ AccT const phi = idir * 2.0 * M_PI/size;
+
+ for (index_type w=0; w<size; ++w)
+ {
+ vsip::complex<AccT> sum = vsip::complex<AccT>();
+ for (index_type k=0; k<size; ++k)
+ sum += vsip::complex<AccT>(in(k)) * sin_cos<AccT>(phi*k*w);
+ out.put(w, T2(sum));
+ }
+}
+
+
+
+// Reference 1-D multi-DFT algorithm on rows of a matrix.
+
+// Requires:
+// IN to be input Matrix.
+// OUT to be output Matrix, of same size as IN.
+// IDIR to be sign of exponential.
+// -1 => Forward Fft,
+// +1 => Inverse Fft.
+
+template <typename T,
+ typename Block1,
+ typename Block2>
+void dft_x(
+ vsip::Matrix<vsip::complex<T>, Block1> in,
+ vsip::Matrix<vsip::complex<T>, Block2> out,
+ int idir)
+{
+ test_assert(in.size(0) == out.size(0));
+ test_assert(in.size(1) == out.size(1));
+ test_assert(in.local().size(0) == out.local().size(0));
+ test_assert(in.local().size(1) == out.local().size(1));
+
+ for (vsip::index_type r=0; r < in.local().size(0); ++r)
+ dft(in.local().row(r), out.local().row(r), idir);
+}
+
+
+
+// Reference 1-D multi-DFT algorithm on columns of a matrix.
+
+template <typename T,
+ typename Block1,
+ typename Block2>
+void dft_y(
+ vsip::Matrix<vsip::complex<T>, Block1> in,
+ vsip::Matrix<vsip::complex<T>, Block2> out,
+ int idir)
+{
+ test_assert(in.size(0) == out.size(0));
+ test_assert(in.size(1) == out.size(1));
+ test_assert(in.local().size(0) == out.local().size(0));
+ test_assert(in.local().size(1) == out.local().size(1));
+
+ for (vsip::index_type c=0; c < in.local().size(1); ++c)
+ dft(in.local().col(c), out.local().col(c), idir);
+}
+
+
+template <typename T,
+ typename Block1,
+ typename Block2>
+void dft_y_real(
+ vsip::Matrix<T, Block1> in,
+ vsip::Matrix<vsip::complex<T>, Block2> out)
+{
+ test_assert(in.size(0) == out.size(0));
+ test_assert(in.size(1) == out.size(1));
+ test_assert(in.local().size(0) == out.local().size(0));
+ test_assert(in.local().size(1) == out.local().size(1));
+
+ for (vsip::index_type c=0; c < in.local().size(1); ++c)
+ dft(in.local().col(c), out.local().col(c), -1);
+}
+
+} // namespace ref
+
+#endif // VSIP_REF_DFT_HPP
Index: tests/window.cpp
===================================================================
RCS file: /home/cvs/Repository/vpp/tests/window.cpp,v
retrieving revision 1.3
diff -u -r1.3 window.cpp
--- tests/window.cpp 20 Dec 2005 12:48:41 -0000 1.3
+++ tests/window.cpp 3 Mar 2006 18:05:40 -0000
@@ -17,6 +17,12 @@
#include <vsip/vector.hpp>
#include "test.hpp"
+#if VSIP_IMPL_SAL_FFT
+# define TEST_NON_POWER_OF_2 0
+#else
+# define TEST_NON_POWER_OF_2 1
+#endif
+
using namespace vsip;
/***********************************************************************
@@ -103,6 +109,7 @@
}
#if defined(VSIP_IMPL_FFT_USE_FLOAT)
+# if TEST_NON_POWER_OF_2
// Chebyshev
{
const length_type N = 24;
@@ -123,6 +130,7 @@
for ( unsigned int n = 0; n < N; ++n )
test_assert( equal( v.get(n), testvec_cheby_odd[n] ) );
}
+# endif
#endif
// Hanning
