[Rcpp-devel] Integrating RcppParallel and RcppEigen

[JJ Allaire]

The initialization of your const fields needs to happen within the
initializer list (rather than inline where the fields are declared).
If you substitute this code for the code in your example everything
will compile:

  // QR Decomposition
  const CPivQR PQR;
  const Permutation Pmat;
  const int r;

  // output
  RMatrix<double> betamat;
  Eigen::VectorXd betahat;

  // initialize with input and output
  CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X,
Rcpp::NumericMatrix betamat)
    : y(y), X(X), PQR(X), Pmat(PQR.colsPermutation()), r(PQR.rank()),
betamat(betamat) {}



On Tue, Apr 14, 2015 at 4:41 AM, Joshua Wiley <jwiley.psych at gmail.com> wrote:
> Hi,
>
> I'm interested in combining faster linear model estimation via QR
> decompositions from RcppEigen with JJs new(er) RcppParallel package to do
> cross validation and bootstrapping on linear models.
>
> As a first step, I've been trying to merge the JSS paper on RcppEigen (
> http://www.jstatsoft.org/v52/i05/paper ) with about the documentation for
> RcppParallel ( http://rcppcore.github.io/RcppParallel/#examples )
>
> ignoring bootstrapping and cross validation for the time being and just
> getting a linear model to run in the parallel framework.  I've toyed with
> each separately successfuly, and think I am somewhat close to getting them
> together, but I am running into an error:
>
> partest.cpp:28:20: error: 'X' is not a type
> partest.cpp:29:26: error: 'PQR' is not a type
> partest.cpp:29:29: error: expected ',' or '...' before '.' token
> partest.cpp:30:15: error: 'PQR' is not a type
> partest.cpp:30:18: error: expected ',' or '...' before '.' token
> partest.cpp: In member function 'virtual void CVLM::operator()(std::size_t,
> std::size_t)':
> partest.cpp:44:28: warning: comparison between signed and unsigned integer
> expressions [-Wsign-compare]
> make: *** [partest.o] Error 1
>
> My code is below (run via sourceCpp() ), and I can see exactly where the
> first error occurs, but I don't understand why this is OK when run outside
> of RcppParallel, but so very wrong in conjunction with RcppParallel.
>
> Any pointers on either what's wrong or where I should be reading would be
> deeply appreciated.
>
> Josh (code below)
>
>
> // [[Rcpp::depends(RcppParallel)]]
> // [[Rcpp::depends(RcppEigen)]]
> #include <Rcpp.h>
> #include <RcppEigen.h>
> #include <RcppParallel.h>
>
> using namespace std;
> using namespace Rcpp;
> using namespace RcppParallel;
>
> using Eigen::Lower;
> using Eigen::Map;
> using Eigen::MatrixXd;
> using Eigen::Upper;
> using Eigen::VectorXd;
> typedef Map<MatrixXd> MapMatd;
> typedef Map<VectorXd> MapVecd;
> typedef Eigen::ColPivHouseholderQR<MatrixXd> CPivQR;
> typedef CPivQR::PermutationType Permutation;
>
> struct CVLM : public Worker
> {
>   // design matrix and outcome
>   const Eigen::VectorXd y;
>   const Eigen::MatrixXd X;
>
>   // QR Decomposition
>   const CPivQR PQR(X); // ERROR HAPPENS HERE
>   const Permutation Pmat(PQR.colsPermutation());
>   const int r(PQR.rank());
>
>
>   // output
>   RMatrix<double> betamat;
>   Eigen::VectorXd betahat;
>
>   // initialize with input and output
>   CVLM(const Eigen::VectorXd y, const Eigen::MatrixXd X, Rcpp::NumericMatrix
> betamat)
>     : y(y), X(X), betamat(betamat) {}
>
>   void operator()(std::size_t begin, std::size_t end) {
>
>     //betahat = PQR.solve(y);
>     for(int j = begin; j < end; j++) {
>       betamat(j, j) = X(j, j) + y[j];
>     }
>   }
> };
>
> // [[Rcpp::export]]
> NumericMatrix parallelFit(Eigen::VectorXd dv, Eigen::MatrixXd x) {
>
>   // allocate the output matrix
>   NumericMatrix betamat(x.rows(), x.cols());
>
>   // pass input and output
>   CVLM cvLM(dv, x, betamat);
>
>   // parallelFor to do it
>   parallelFor(0, x.rows(), cvLM);
>
>   return betamat;
> }
>
>
>
>
>
>
>
>
>
>
>
>
> --
> Joshua F. Wiley
> http://joshuawiley.com/
> ---
> Postdoctoral Research Fellow
> Mary MacKillop Institute for Health Research
> Australian Catholic University
> ---
> Senior Analyst, Elkhart Group Ltd.
> http://elkhartgroup.com
> Office: 260.673.5518
>
> _______________________________________________
> 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-devel