[Rcpp-devel] Multiplication of ComplexVector?

[baptiste auguie]

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.


Thanks,

baptiste