[Rcpp-devel] RcppMAtrix<complex> ?

Davor Cubranic cubranic at stat.ubc.ca
Fri Apr 23 19:55:04 CEST 2010 

Baptiste,

Do have a look at RcppArmadillo, it makes matrix code very easy to write 
(and read!). For example, this is a function callable from R to add two 
matrices:

#include <RcppArmadillo.h>

using namespace Rcpp;

RcppExport SEXP add_mat(SEXP y_in,
			SEXP x_in) {
  NumericMatrix y_r(y_in);
  arma::mat y(y_r.begin(), y_r.nrow, y_r.ncol(), false);
  
  NumericMatrix x_r(x_in);
  arma::mat x(x_r.begin(), x_r.nrow(), x_r.ncol(), false);

  arma::mat result = x+y;
  return wrap(result);
}

In this case, the matrices are real, not complex. but I believe simply 
replacing 'arma::mat' with 'arma::cx_mat' and 'NumericMatrix' with 
'ComplexMatrix' will work for the complex case.

Davor


On April 22, 2010 01:22:58 pm baptiste auguie wrote:
> OK, thanks for the information. I guess I was lead to believe that
> such operations were of common use because I worked previously with
> colleagues who had defined their own C++ complex class (and I'm
> guessing it was precisely for this purpose of operator overloading).
> Sadly the code was not open source. I'll look into armadillo,
> hopefully it provides an alternative. I might also need to rethink
> what portion of the code I should really be porting out of R; this
> function was clearly the bottleneck of my code but it looks like it
> will be painful to write it in a lower-level language.
> 
> Thanks,
> 
> baptiste
> 
> On 22 April 2010 18:42, Romain Francois <romain at r-enthusiasts.com> 
wrote:
> > Le 22/04/10 18:19, baptiste auguie a écrit :
> >> Thanks for the example, it is useful indeed. However, I am having
> >> difficulties with operations involving whole matrices, rather than
> >> matrix elements.
> >
> > We don't currently have those.
> >
> >> Again I must warn that I don't know C++ ; but the
> >> addition of two matrices does not seem to work out-of-the-box, as
> >> well as more complicated functions. The dispatch of these
> >> functions might not exist for the complex matrix class, or maybe
> >> it is not implemented in Rcpp?
> >>
> >> The operations I would need to perform with complex matrices are,
> >>
> >> +, -, *, transpose,
> >> as well as operations on 1-column and one-row matrices (==
> >> vectors?) such as exp().
> >
> > Those are things you'd typically do in R, not in C/C++
> >
> >> Working component by component is not a very attractive option
> >
> > That is what you usually do in a C/C++ world
> >
> >> so I'm
> >> hoping there is an easy way to define operations between matrices,
> >> matrices and vectors, and matrices and scalars.
> >
> > One thing you can do perhaps is look into armadillo, which we wrap
> > nicely with the RcppArmadillo package. and thanks to Doug, the
> > wrapping is even nicer now since armadillo is packed up inside
> > RcppArmadillo (but this version is not released yet).
> >
> > We will not implement these operators soon because it is very easy
> > to not do it right. armadillo does it nicely.
> >
> >> Thanks,
> >>
> >> baptiste
> >>
> >> PS: The problem with my previous email was in the gmail vs.
> >> googlemail domain, as Dirk pointed out (I had had that same
> >> problem before, but I forgot!).
> >>
> >> On 22 April 2010 13:17, Romain Francois<romain at r-enthusiasts.com> 
 wrote:
> >>> Thank you for reposting here.
> >>>
> >>> It is not trivial to see what is happening in your example, so
> >>> I'll just give you some tools.
> >>>
> >>> The old api (which the classicRcppMatrixExample example uses)
> >>> does not have
> >>> support for complex vectors or matrices.
> >>>
> >>> The new api does have support for complex vectors and complex
> >>> matrices. The
> >>> unit test file runit.ComplexVector.R does indeed contain some
> >>> very basic examples of using ComplexVector, but not complex
> >>> matrices. However, you can
> >>> use Rcpp::ComplexMatrix.
> >>>
> >>> Here is an example that calculates the sum of the real part of
> >>> the elements
> >>> of a complex matrix diagonal and the sum of the imaginary part:
> >>>
> >>> require( Rcpp )
> >>> require( inline)
> >>>
> >>> fx<- cfunction( signature( x = "matrix" ), '
> >>>        /* grab the R object as a complex matrix */
> >>>        ComplexMatrix m(x) ;
> >>>        double re_sum = 0.0 ;
> >>>        double im_sum = 0.0 ;
> >>>        for( int i=0; i<m.ncol(); i++){
> >>>                re_sum += m(i,i).r ;
> >>>                im_sum += m(i,i).i ;
> >>>        }
> >>>        return List::create(
> >>>                _["sum real part"] = re_sum,
> >>>                _["sum imag part"] = im_sum
> >>>                ) ;
> >>>
> >>> ', Rcpp = TRUE, includes = "using namespace Rcpp;" )
> >>>
> >>> x<- diag( (1-2i)*1:5 )
> >>> fx( x )
> >>>
> >>> Let us know if this gives you enough to get started.
> >>>
> >>> Romain
> >>>
> >>> Le 22/04/10 12:59, baptiste auguie a écrit :
> >>>> Dear all,
> >>>>
> >>>> I'm hoping to port some R code to C++ to make it faster. The
> >>>> code makes heavy use of matrices of complex numbers, and playing
> >>>> with the RcppExamples package this morning I got the impression
> >>>> that it's not currently implemented in Rcpp. I basically took
> >>>> the example from classicRcppMatrixExample and tried to change
> >>>> the types from double to complex. I must confess that I don't
> >>>> know C++, so maybe I missed something obvious.
> >>>>
> >>>> Attached is my dummy example, as well as the R code I'm trying
> >>>> to port to give you an idea. Any suggestions would be greatly
> >>>> appreciated!
> >>>>
> >>>> Best regards,
> >>>>
> >>>> baptiste
> >>>>
> >>>>
> >>>> PS: This message initially failed to reach the list; in the
> >>>> meantime I got the suggestion from Romain and Dirk to have a
> >>>> look at the class Rcpp::ComplexMatrix and the example in
> >>>> runit.ComplexVector.R
> >>>
> >>> --
> >>> Romain Francois
> >>> Professional R Enthusiast
> >>> +33(0) 6 28 91 30 30
> >>> http://romainfrancois.blog.free.fr
> >>>
> >>> |- http://bit.ly/9aKDM9 : embed images in Rd documents
> >>> |- http://tr.im/OIXN : raster images and RImageJ
> >>> |- http://tr.im/OcQe : Rcpp 0.7.7
> >
> > --
> > Romain Francois
> > Professional R Enthusiast
> > +33(0) 6 28 91 30 30
> > http://romainfrancois.blog.free.fr
> >
> > |- http://bit.ly/9aKDM9 : embed images in Rd documents
> > |- http://tr.im/OIXN : raster images and RImageJ
> > |- http://tr.im/OcQe : Rcpp 0.7.7
> 
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-dev
> el
>

More information about the Rcpp-devel mailing list