[Rcpp-devel] Multiplication of ComplexVector?

[Romain Francois]

Le 17/08/10 10:00, baptiste auguie a écrit :
> On 17 August 2010 09:20, Romain Francois<romain at r-enthusiasts.com>  wrote:
>> Le 17/08/10 07:43, baptiste auguie a écrit :
>>>
>>> Hi,
>>>
>>> On 17 August 2010 03:24, Dirk Eddelbuettel<edd at debian.org>    wrote:
>>>>
>>>> On 16 August 2010 at 19:43, Dirk Eddelbuettel wrote:
>>>
>>> [...]
>>>
>>>> The only trouble is that nobody has written the corresponding 'glue' code
>>>> to
>>>> make
>>>>
>>>>       arma::cx_vec a1(y1.begin(), y1.size(), false);
>>>>
>>>> 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.
>>>
>>> I'm confused, isn't
>>>
>>>   arma::cx_colvec a1 = Rcpp::as<    arma::cx_vec>( y1 );
>>>
>>> aimed at doing this kind of conversion? (that's what I use, following
>>> Romain's tip)
>>>
>>> baptiste
>>
>> I don't think you are confused. This works well, but it does require some
>> extra copying, which Dirk was trying to spare.
>>
>> Using the "advanced" constructor in armadillo means that both the
>> ComplexVector and the cx_vec use the same memory. although it says loud and
>> clear in armadillo documentation :
>>
>> """
>> This is faster, but can be dangerous unless you know what you're doing!
>> """
>
>
> Thanks. Two questions then,
>
> 1- reinterpret_cast<  std::complex<double>*>  seems like a good option,
> does it work with real / complex matrices as well? The only things
> that came up from googling it were related to RcppGSL, not
> RcppArmadillo.

Sure.

code<- '
     ComplexMatrix y1(i), y2(ii);
     arma::cx_mat a1( reinterpret_cast< std::complex<double>* 
 >(y1.begin()), y1.nrow(), y1.ncol(), false);
     arma::cx_mat a2( reinterpret_cast< std::complex<double>* 
 >(y2.begin()), y2.nrow(), y2.ncol(), false);
     return wrap( a1 + a2);
     '
fun<- cxxfunction(signature(i="complex", ii = "complex"), code, 
plugin="RcppArmadillo")

x1 <- matrix( 1:10*(1+1i), 5, 2 )
x2 <- matrix( 1*1:10*(1+1i), 5, 2 )
fun( x1, x2 )

You don't need reinterpret_cast for converting NumericMatrix to a 
Mat<double> aka mat because they both use the same type and layout for 
storage : double with column major order layout

> 2- when is it not safe? ;)

When people don't know what they are doing.

For example : R has no knowledge that armadillo shares the memory with 
it, so if for some reason R thinks it can garbage collect the SEXP that 
came in ("i" and "ii"), it won't ask around for permission and so the 
armadillo matrix/vector will contain a pointer to wonderland, and then 
killing the jabbawockee bug is not easy.

So if all you do is share the memory inside the same function/scope you 
should be fine.

> Thanks!
>
> baptiste

And sorry it did not occur to me before.

> PS: i wouldn't want to hijack the original query; I'll open a new
> thread if this seems to raise further questions.

-- 
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