|
|
Ublas : |
Subject: Re: [ublas] Want to allocate my own C arrays
From: Heiko Bauke (heiko.bauke_at_[hidden])
Date: 2010-04-09 10:22:46
Hi,
On Thu, 8 Apr 2010 15:29:04 -0400
Paul Dugas <paul_at_[hidden]> wrote:
> I'm poking around to see if I can find a way to marry uBLAS vector and
> matrix data with the DFT implementation provided by the FFTW library.
> FFTW operates on 1, 2 and 3 dimensional C-style arrays of real and
> complex values. There's also support for multi- operations that can
> perform DFTs on multiple arrays at the same time; i.e. arrays and
> arrays and arrays of matrices. There is a C++ wrapper for FFTW but
> it's still using the C-style arrays for the data. The algorithms I
> need to implement would ideally utilize uBLAS for much of the math and
> call wrapped FFTW routines to perform DFTs.
>
> I started looking at the c_vector<> and c_matrix<> templates. The way
> they expose the underlying C-array with the data() method seemed
> promising. Then I noticed it's not new'ing the array. Turns out FFTW
> operates better if the values are aligned just so in memory. Is there
> something along the lines of c_vector<> that takes an allocator
> template parameter?
it is not difficult to transform data of a ublas vector or matrix. Just
determine the address of the first element and call the FFTW function:
typedef std::complex<double> complex;
typedef boost::numeric::ublas::vector<complex> complex_vector;
int N=1024;
complex_vector v_in(N), v_out(N);
fftw_plan p;
p=fftw_plan_dft_1d(N,
reinterpret_cast<fftw_complex *>(&v_in(0)),
reinterpret_cast<fftw_complex *>(&v_out(0)),
FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
In order to align data in memory use a custom allocator. The following
code snipped utilizes an allocator that I have written for my own FFTW
wrapper library:
#include <cstdlib>
#include <cstddef>
#include <iostream>
#include <limits>
#include <boost/numeric/ublas/vector.hpp>
#include <fftw3.h>
namespace fftw {
template<typename T>
struct traits;
template<>
struct traits<float> {
typedef float float_type;
typedef float real_complex_type;
};
template<>
struct traits<double> {
typedef double float_type;
typedef double real_complex_type;
};
template<>
struct traits<long double> {
typedef long double float_type;
typedef long double real_complex_type;
};
template<>
struct traits<std::complex<float> > {
typedef float float_type;
typedef fftwf_complex real_complex_type;
};
template<>
struct traits<std::complex<double> > {
typedef double float_type;
typedef fftw_complex real_complex_type;
};
template<>
struct traits<std::complex<long double> > {
typedef long double float_type;
typedef fftwl_complex real_complex_type;
};
template<>
struct traits<fftwf_complex> {
typedef float float_type;
typedef fftwf_complex real_complex_type;
};
template<>
struct traits<fftw_complex> {
typedef double float_type;
typedef fftw_complex real_complex_type;
};
template<>
struct traits<fftwl_complex> {
typedef long double float_type;
typedef fftwl_complex real_complex_type;
};
namespace detail{
inline void * malloc(::size_t n, float) {
return fftwf_malloc(n);
}
inline void * malloc(::size_t n, double) {
return fftw_malloc(n);
}
inline void * malloc(::size_t n, long double) {
return fftwl_malloc(n);
}
inline void free(void *p, float) {
return fftwf_free(p);
}
inline void free(void *p, double) {
return fftw_free(p);
}
inline void free(void *p, long double) {
return fftwl_free(p);
}
}
template<typename T>
class allocator {
public:
// type definitions
typedef T value_type;
typedef T * pointer;
typedef const T * const_pointer;
typedef T & reference;
typedef const T & const_reference;
typedef std::size_t size_type;
typedef std::ptrdiff_t difference_type;
// rebind allocator to type U
template<typename U>
struct rebind {
typedef allocator<U> other;
};
// return address of values
pointer address(reference value) const {
return &value;
}
const_pointer address(const_reference value) const {
return &value;
}
/* constructors and destructor
* - nothing to do because the allocator has no state
*/
allocator() throw() {
}
allocator(const allocator &) throw() {
}
template<typename U>
allocator(const allocator<U> &) throw() {
}
~allocator() throw() {
}
// return maximum number of elements that can be allocated
size_type max_size() const throw() {
return std::numeric_limits<std::size_t>::max()/sizeof(T);
}
// allocate but don't initialize num elements of type T
pointer allocate(size_type num, const void *hint=0) {
pointer ret=reinterpret_cast<pointer>(detail::malloc(num*sizeof(T),
typename traits<T>::float_type()));
return ret;
}
// initialize elements of allocated storage p with value value
void construct(pointer p, const T &value) {
// initialize memory with placement new
new(reinterpret_cast<void *>(p)) T(value);
}
// destroy elements of initialized storage p
void destroy(pointer p) {
// destroy objects by calling their destructor
p->~T();
}
// deallocate storage p of deleted elements
void deallocate(pointer p, size_type num) {
detail::free(reinterpret_cast<void *>(p),
typename traits<value_type>::float_type());
}
};
// return that all specializations of this allocator are interchangeable
template<typename T1, class T2>
bool operator==(const allocator<T1> &, const allocator<T2>&) throw() {
return true;
}
template<typename T1, class T2>
bool operator!=(const allocator<T1>&, const allocator<T2>&) throw() {
return false;
}
}
int main() {
typedef std::complex<double> complex;
typedef boost::numeric::ublas::unbounded_array<complex, fftw::allocator<complex> > complex_storage_type;
typedef boost::numeric::ublas::vector<complex, complex_storage_type> complex_vector;
int N=1024;
complex_vector v_in(N), v_out(N);
fftw_plan p;
p=fftw_plan_dft_1d(N,
reinterpret_cast<fftw_complex *>(&v_in(0)),
reinterpret_cast<fftw_complex *>(&v_out(0)),
FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
return EXIT_SUCCESS;
}
Heiko
-- -- Was man vergißt, hat man im Grunde nicht erlebt. -- (Ernst R. Hauschka, dt.Aphoristiker, 1926-) -- Number Crunch Blog @ http://numbercrunch.de -- Heiko Bauke @ http://www.mpi-hd.mpg.de/personalhomes/bauke