[Rcpp-devel] Multiplication of ComplexVector?

[Romain Francois]

Le 17/08/10 23:09, baptiste auguie a écrit :
> Hi,
>
> On 17 August 2010 10:04, Romain Francois<romain at r-enthusiasts.com>  wrote:
>> 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
>
> Thanks for the nice examples.
>
>>
>>> 2- when is it not safe? ;)
>>
>> When people don't know what they are doing.
>
> Oops, that'd be me :)
>
>>
>> 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.
>>
>
> I'll need to get bitten by the bug before I can really absorb the
> meaning of this. I basically have no clue as to who owns the data /
> pointers, underneath layers of Rcpp and Armadillo.

Rcpp runs the show, so as long as your Rcpp object lives longer than the 
armadillo object, you are fine.

armadillo has no clue it gets the data from Rcpp and will not interfere, 
it will not attempt to delete the data or whatever, just access it.

The way RcppArmadillo is usually used : "get SEXP, convert to Rcpp data 
(e.g. NumericVector), create a armadillo object, do some stuff with the 
armadillo object, convert the result back to Rcpp somehow and return 
something to R", this should be ok.

Romain

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