[Bug c++/12886] New: program compiles incorrectly without optimization on

["gcc at schabel-family dot org" <gcc-bugzilla at gcc dot gnu dot org>]

PLEASE REPLY TO gcc-bugzilla@gcc.gnu.org ONLY, *NOT* gcc-bugs@gcc.gnu.org.

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=12886

           Summary: program compiles incorrectly without optimization on
           Product: gcc
           Version: 3.3
            Status: UNCONFIRMED
          Severity: normal
          Priority: P2
         Component: c++
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: gcc at schabel-family dot org
                CC: gcc-bugs at gcc dot gnu dot org

The code below, when compiled, produces all zeros unless optimization is turned on.  
With any optimization level (-O1 through -O4), the code runs as expected.

running on Apple PowerBook, Mac OS 10.3

g++ --version
g++ (GCC) 3.3 20030304 (Apple Computer, Inc. build 1495)
Copyright (C) 2002 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

command line arguments to compiler:

g++ test_fft.cpp -I../.. -lm -lfftw3               -> produces all zeroes
g++ test_fft.cpp -I../.. -O1 -lm -lfftw3        -> produces expected output
g++ test_fft.cpp -I../.. -O2 -lm -lfftw3        -> "
g++ test_fft.cpp -I../.. -O3 -lm -lfftw3        -> "
g++ test_fft.cpp -I../.. -O4 -lm -lfftw3        -> "

I am linking with the 1.30 release of boost (www.boost.org) and the 3.0.1 release of FFTW 
(www.fftw.org) - the latter was built in the standard way : ./configure, make, make install.

// code below
#include <bitset>
#include <complex>
#include <iostream>

#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>

#include <fftw3.h>

namespace fft {

struct Types
{
	typedef unsigned int				flag_type;
	typedef int							integer_type;
	typedef fftw_r2r_kind_do_not_use_me	r2r_type;
	typedef int							size_type;
};

// complex transforms
struct Forward 
{
	static const Types::integer_type 	kind()	{ return FFTW_FORWARD; }
};

struct Inverse 
{
	static const Types::integer_type	kind() 	{ return FFTW_BACKWARD; } 
};

// complex discrete Fourier transform
template<class DT = Forward,
		 Types::size_type R = 1>
struct DFT : 
	public DT
{ 
	typedef DT		direction_type;
	
	typedef Types::integer_type		integer_type;
	typedef Types::size_type		size_type;
	typedef Types::flag_type		flag_type;
	
	static const Types::size_type	Rank = R;
};		

template<class T,class P> struct Dispatcher;

// DFT
template<>
struct Dispatcher<double,DFT<Forward,1> >
{
	typedef DFT<Forward,1>	transform_type;
	
	typedef transform_type::flag_type		flag_type;
	typedef transform_type::integer_type	integer_type;
	typedef transform_type::size_type		size_type;
	
	typedef fftw_plan							plan_type;
	typedef double								value_type;
	typedef std::complex<double>				input_type;
	typedef std::complex<double>				output_type;
	typedef input_type*							input_pointer_type;
	typedef output_type*						output_pointer_type;
	typedef fftw_complex*						fftw_pointer_type;
	
	static plan_type dispatch(const size_type size[],input_pointer_type in,output_pointer_type 
out,flag_type flags)
	{
		return fftw_plan_dft_1d(size[0],
								reinterpret_cast<fftw_pointer_type>(in),
								reinterpret_cast<fftw_pointer_type>(out),
								transform_type::kind(),
								flags);
	}
};

template<>
struct Dispatcher<double,DFT<Inverse,1> >
{
	typedef DFT<Inverse,1>	transform_type;
	
	typedef transform_type::flag_type		flag_type;
	typedef transform_type::integer_type	integer_type;
	typedef transform_type::size_type		size_type;

	typedef fftw_plan							plan_type;
	typedef double								value_type;
	typedef std::complex<double>				input_type;
	typedef std::complex<double>				output_type;
	typedef input_type*							input_pointer_type;
	typedef output_type*						output_pointer_type;
	typedef fftw_complex*						fftw_pointer_type;
	
	static plan_type dispatch(const size_type size[],input_pointer_type in,output_pointer_type 
out,flag_type flags)
	{
		return fftw_plan_dft_1d(size[0],
								reinterpret_cast<fftw_pointer_type>(in),
								reinterpret_cast<fftw_pointer_type>(out),
								transform_type::kind(),
								flags);
	}
};

template<>
struct Dispatcher<double,DFT<Forward,2> >
{
	typedef DFT<Forward,2>	transform_type;
	
	typedef transform_type::flag_type		flag_type;
	typedef transform_type::integer_type	integer_type;
	typedef transform_type::size_type		size_type;
	
	typedef fftw_plan							plan_type;
	typedef double								value_type;
	typedef std::complex<double>				input_type;
	typedef std::complex<double>				output_type;
	typedef input_type*							input_pointer_type;
	typedef output_type*						output_pointer_type;
	typedef fftw_complex*						fftw_pointer_type;
	
	static plan_type dispatch(const size_type size[],input_pointer_type in,output_pointer_type 
out,flag_type flags)
	{
		return fftw_plan_dft_2d(size[0],size[1],
								reinterpret_cast<fftw_pointer_type>(in),
								reinterpret_cast<fftw_pointer_type>(out),
								transform_type::kind(),
								flags);
	}
};

template<>
struct Dispatcher<double,DFT<Inverse,2> >
{
	typedef DFT<Inverse,2>	transform_type;
	
	typedef transform_type::flag_type		flag_type;
	typedef transform_type::integer_type	integer_type;
	typedef transform_type::size_type		size_type;

	typedef fftw_plan							plan_type;
	typedef double								value_type;
	typedef std::complex<double>				input_type;
	typedef std::complex<double>				output_type;
	typedef input_type*							input_pointer_type;
	typedef output_type*						output_pointer_type;
	typedef fftw_complex*						fftw_pointer_type;
	
	static plan_type dispatch(const size_type size[],input_pointer_type in,output_pointer_type 
out,flag_type flags)
	{
		return fftw_plan_dft_2d(size[0],size[1],
								reinterpret_cast<fftw_pointer_type>(in),
								reinterpret_cast<fftw_pointer_type>(out),
								transform_type::kind(),
								flags);
	}
};

// FFTW support
template<class T> struct FFTW;

template<>
struct FFTW<double>
{
	typedef double		value_type;
	typedef fftw_plan	plan_type;
	
	static void cleanup()								
											{ return fftw_cleanup(); }
	static void destroy_plan(plan_type plan)				
										{ return fftw_destroy_plan(plan); }
	static void execute(plan_type plan)					
											{ return fftw_execute(plan); }
	static void free(void *ptr)							
											{ return fftw_free(ptr); }
	static void* malloc(int size)						
											{ return fftw_malloc(size); }
};

// allocator
template<class D>
class Allocator
{
	public:
		typedef D	dispatcher_type;
		
		typedef typename dispatcher_type::size_type		
		size_type;
		typedef typename dispatcher_type::value_type		
	value_type;
		typedef typename dispatcher_type::input_type		
	input_type;
		typedef typename dispatcher_type::output_type	
		output_type;
		typedef typename dispatcher_type::input_pointer_type 
	input_pointer_type;
		typedef typename dispatcher_type::output_pointer_type
	output_pointer_type;
		
		Allocator();
		~Allocator();

		void reserve(const size_type sz)	{ size_ = sz; }
		
		input_pointer_type in()				{ allocate(); return in_; }
		output_pointer_type out()			{ allocate(); return out_; }
		
	private:
		bool allocated() const				{ return allocated_; }
		
		void allocate();
		void deallocate();
		
		size_type			size_;
		bool				allocated_;
		input_pointer_type	in_;
		output_pointer_type	out_;
};

template<class D>
Allocator<D>::Allocator() :
	size_(0),
	allocated_(false),
	in_(NULL),
	out_(NULL)
{ }

template<class D>
Allocator<D>::~Allocator()
{
	deallocate();
}

template<class D>
void
Allocator<D>::allocate()
{
	if (!allocated())
	{
		in_ = reinterpret_cast<input_pointer_type>(FFTW<value_type>::
malloc(size_*sizeof(input_type)));
		out_ = reinterpret_cast<output_pointer_type>(FFTW<value_type>::
malloc(size_*sizeof(output_type)));
		
		allocated_ = true;
	}
}

template<class D>
void
Allocator<D>::deallocate()
{
	if (allocated())
	{
		FFTW<value_type>::free(reinterpret_cast<void*>(in_));
		FFTW<value_type>::free(reinterpret_cast<void*>(out_));
		
		size_ = 0;
		allocated_ = false;
		in_ = out_ = NULL;
	}	
}

// FFTW plan
template<class T,
		 class P,
		 class A = Allocator<Dispatcher<T,P> > >
class FFT
{
	public:
		typedef FFT<T,P,A>	this_type;
		typedef T			value_type;
		typedef P			transform_type;
		typedef A			allocator_type;
		
		static const typename transform_type::size_type	Rank = transform_type::Rank;
		
		typedef typename transform_type::direction_type	
		direction_type;
		
		typedef Dispatcher<value_type,transform_type>	
		dispatcher_type;
		typedef typename dispatcher_type::integer_type	
		integer_type;
		typedef typename dispatcher_type::size_type		
		size_type;
		typedef typename dispatcher_type::flag_type		
		flag_type;
		
		typedef typename dispatcher_type::plan_type		
		plan_type;
		typedef typename dispatcher_type::input_type		
	input_type;
		typedef typename dispatcher_type::output_type	
		output_type;
		typedef typename dispatcher_type::input_pointer_type
	input_pointer_type;
		typedef typename dispatcher_type::output_pointer_type
	output_pointer_type;
		
		FFT(size_type sz);
		FFT(size_type sz0,size_type sz1);
		FFT(size_type sz0,size_type sz1,size_type sz2);
		FFT(const size_type sz[]);
		~FFT();
		
		input_type& input(const size_type idx)			
		{ return alloc_.in()[idx]; }
		const input_type& input(const size_type idx) const	
	{ return alloc_.in()[idx]; }
		output_type& output(const size_type idx)			
	{ return alloc_.out()[idx]; }
		const output_type& output(const size_type idx) const
	{ return alloc_.out()[idx]; }
		
		void estimate()								
			{ set_strategy(FFTW_ESTIMATE); }
		void measure()								
			{ set_strategy(FFTW_MEASURE); }
		void patient()								
			{ set_strategy(FFTW_PATIENT); }
		void exhaustive()							
			{ set_strategy(FFTW_EXHAUSTIVE); }
	
		void execute()								
			{ FFTW<value_type>::execute(get_plan()); }	
	
		
	 	size_type size() const						
			{ return alloc_.size(); }	
		
	private:
		BOOST_STATIC_ASSERT((boost::is_same<input_type,typename allocator_type::
input_type>::value));
		BOOST_STATIC_ASSERT((boost::is_same<output_type,typename allocator_type::
output_type>::value));
		
		FFT(const this_type&);
		this_type& operator=(const this_type&);
		
		// internal state
		enum 
		{ 
			PLAN_VALID = 0
		};
			
		void clear()								
			{ if (valid()) FFTW<value_type>::cleanup(); options_.reset(); state_.reset(); }
		void set_strategy(const flag_type);
		
		bool valid() const							
			{ return state_[PLAN_VALID]; }
		
		plan_type get_plan();
		
		plan_type			plan_;
		size_type			size_[Rank];
		allocator_type		alloc_;
		flag_type			flags_;
		
		std::bitset<8>		options_,
							state_;
};

template<class T,
		 class P,
		 class A>
FFT<T,P,A>::FFT(size_type sz)
{ 
	BOOST_STATIC_ASSERT(Rank == 1);
	
	size_[0] = sz;
	
	alloc_.reserve(sz);

	clear();
}

template<class T,
		 class P,
		 class A>
FFT<T,P,A>::FFT(size_type sz0,size_type sz1)
{
	BOOST_STATIC_ASSERT(Rank == 2);
	
	size_[0] = sz0;
	size_[1] = sz1;
	
	alloc_.reserve(sz0*sz1);
	
	clear();
}

template<class T,
		 class P,
		 class A>
FFT<T,P,A>::FFT(size_type sz0,size_type sz1,size_type sz2)
{
	BOOST_STATIC_ASSERT(Rank == 3);
	
	size_[0] = sz0;
	size_[1] = sz1;
	size_[2] = sz2;
	
	alloc_.reserve(sz0*sz1*sz2);

	clear();
}

template<class T,
		 class P,
		 class A>
FFT<T,P,A>::FFT(const size_type sz[])
{ 
	size_type	tmp = 1;
	
	for (size_type i=0;i< Rank;++i) 
	{
		size_[i] = sz[i];
		tmp *= sz[i];
	}
	
	alloc_.reserve(tmp);
	
	clear();
}

template<class T,
		 class P,
		 class A>
FFT<T,P,A>::~FFT()
{
	if (valid()) 
	{
		FFTW<value_type>::destroy_plan(get_plan());
	}
		
	flags_ = 0;
	
	options_.reset();
	state_.reset();
}

template<class T,
		 class P,
		 class A>
inline void
FFT<T,P,A>::set_strategy(const flag_type f)
{
	if (!(flags_ & f))
	{
		flags_ &= ~(FFTW_ESTIMATE | FFTW_MEASURE | FFTW_PATIENT | FFTW_EXHAUSTIVE);
		flags_ |= f;
		state_.reset();
	}
}

template<class T,
		 class P,
		 class A>
typename FFT<T,P,A>::plan_type 
FFT<T,P,A>::get_plan()		
{ 
	if (!valid())
	{
		plan_ = dispatcher_type::dispatch(size_,alloc_.in(),alloc_.out(),flags_);
		
		if (plan_ == NULL) throw std::exception();
	} 
	
	return plan_;
}

} //  namespace fft

int main()
{
	using namespace fft;
	
    fftw_complex *in, *out;
    fftw_plan p;
    
    int	N = 64;
    
    in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * N));
    out = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * N));
 
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    
    fftw_execute(p); /* repeat as needed */
    
    fftw_destroy_plan(p);

    fftw_free(in); 
    fftw_free(out);
    
    FFT<double,DFT<Forward> >	d128fp(64);
    FFT<double,DFT<Inverse> >	d128ip(64);
    
    for (int i=0;i<64;++i)
    {
    	d128fp.input(i) = std::cos(i/8.0+(3.141592/4));
    }
    
    d128fp.execute();
    d128ip.execute();
    
    for (int i=0;i<64;++i)
    {
    	std::cout << d128fp.input(i) << "\t" 
    			  << d128fp.output(i) << std::endl;
	}

    FFT<double,DFT<Forward,2> >	d128fp2(8,8);
    FFT<double,DFT<Inverse,2> >	d128ip2(8,8);
    
	return 0;
}