[Rcpp-devel] Multiplication of ComplexVector?
Romain Francois
romain at r-enthusiasts.com
Tue Aug 17 09:16:17 CEST 2010
Le 17/08/10 03:24, Dirk Eddelbuettel a écrit :
>
> On 16 August 2010 at 19:43, Dirk Eddelbuettel wrote:
> | Hi Christian,
> |
> | Thanks for your interest in Rcpp, and for posting here.
> |
> | On 16 August 2010 at 17:19, Christian Gunning wrote:
> | | Dear list,
> | |
> | | I'm trying to use the ComplexVector class, and I'm having difficulty
> | | with multiplication. I sugar multiplication, element-by-element
> | | multiplication (no * operator for either), and using temporary
> | | Rcomplex variables (don't really understand what I'm doing here).
> |
> | Well, NumericVector et al aren't really made for all possible math ops done
> | on vectors, real or complex. They are first and foremost storage types. But
> | for the math we do have Armadillo and RcppArmadillo. So in the short run, I
> | would suggest the following two-step procedure:
> |
> | i) have a good look at Armadillo (http://arma.sf.net) and use the cx_mat
> | or cx_vec types for mulitplication; write a small test program to see
> | that everything works
> |
> | ii) use RcppArmadillo to from R via Rcpp / RcppArmadillo to this
> | functionality
> |
> | | I was able to successfully import get("*") as an Rcpp::Function and do
> | | the multiplication through R. Is this the best option for now?
> |
> | I am sure we can do better than this but I am not sure all pieces are in
> | place just yet.
>
> Looks like I wasn't all that far off. This builds:
>
> library(inline)
>
> code<- '
> ComplexVector y1(i), y2(ii);
> //arma::cx_vec a1(y1.begin(), y1.size(), false);
> //arma::cx_vec a2(y2.begin(), y2.size(), false);
> arma::cx_vec a1(y1.size());
> arma::cx_vec a2(y2.size());
> arma::cx_vec a3 = a1 * a2;
> arma::cx_vec a4 = a1 + a2;
> List z = List::create(a1, a2, a3, a4);
> return z;
> '
>
> fun<- cxxfunction(signature(i="compex", ii="comples"), code, plugin="RcppArmadillo")
>
> The only trouble is that nobody has written the corresponding 'glue' code to
> make
>
> arma::cx_vec a1(y1.begin(), y1.size(), false);
cx_mat holds an array of std::complex<double>, not an array of Rcomplex.
so the compiler is unhappy.
However, Rcomplex and std::complex<double> have the same layout in
memory, so you are one reinterpret_cast away :
require( inline )
code<- '
ComplexVector y1(i), y2(ii);
arma::cx_vec a1( reinterpret_cast< std::complex<double>*
>(y1.begin()), y1.size(), false);
arma::cx_vec a2( reinterpret_cast< std::complex<double>*
>(y2.begin()), y2.size(), false);
arma::cx_vec a3 = a1 % a2; // element-wise multiplication in armadillo
arma::cx_vec a4 = a1 + a2;
List z = List::create(a1, a2, a3, a4);
return z;
'
fun<- cxxfunction(signature(i="complex", ii="complex"), code,
plugin="RcppArmadillo")
fun( 1:10*(1+1i), (1:10)*2+2i )
> happen: create an Armadillo complex vector from an Rcpp::ComplexVector. We
> can init by scalar size, what you'd need to insert for now is a simply (and
> very pedestrian) copy-loop.
>
> Hth, Dirk
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/bzoWrs : Rcpp svn revision 2000
|- http://bit.ly/b8VNE2 : Rcpp at LondonR, oct 5th
`- http://bit.ly/aAyra4 : highlight 0.2-2
More information about the Rcpp-devel
mailing list