title: Using arma::solve with RcppParallel: how to avoid stack imbalance?
body: <p>I want to use arma::solve within a rcppParallel worker.
  The problem is that the code results in a crash of R, probably due to stack imbalance.</p>
  
  <p>To avoid stack imbalance issues one should use RMatrix objects instead of arma::mat objects in the worker as far as I understand the decription of the rcppParallel package.
  However arma::solve only accepts arma::mat objects. </p>
  
  <p>How is this done correctly? What is the correct way to deal with this issue?</p>
  
  <p>Thank you very much for any suggestions/help/corrections!</p>
  
  <p>Below you find my approach which unfortunately easily results in crash of R (especially when executed serveral times): </p>
  
  <p>// [[Rcpp::depends(RcppParallel)]]</p>
  
  <pre><code>// [[Rcpp::depends(RcppArmadillo)]]
  #include &lt;RcppArmadillo.h&gt;
  #include &lt;RcppParallel.h&gt;
  #include &lt;iostream&gt;
  
  using namespace Rcpp;
  using namespace RcppParallel;
  
  
  struct MakeITER : public Worker {
    // input matrix to read from
    const arma::mat MAT;
    const arma::colvec VEC;
    const int p;
     // output matrix to write to
    arma::mat res1;
  
    MakeITER( arma::mat res, const arma::mat MAT0, 
    const arma::colvec VEC0, int p0)
    :  res1(res), MAT(MAT0), VEC(VEC0), p(p0){
      res1.zeros();
    }
  
    void operator()(std::size_t begin, std::size_t end) {
      int j, m;   
  
      arma::colvec coef(p);
  
      for(m = begin; m &lt; end; m++){              
        try{
              coef = arma::solve(MAT, VEC);
              for(j=0; j&lt;p; j++) {
                      res1(j,m) = coef[j];
              } ;
        }catch(...){             
              for(j=0; j&lt;p; j++) {
                res1(j,m) =  arma::datum::nan; 
              };
          }   
      }    
    }
  };
  
  // [[Rcpp::export]]
  arma::mat fuc( arma::mat MAT0,  arma::colvec VEC0, int nmkt, int p0) {
    // result:
    arma::mat res(p0,nmkt);  
  
    MakeITER makeit(res, MAT0,  VEC0, p0);
    parallelFor(0, nmkt, makeit);
  
    return(  res );
  
  };
  </code></pre>
  
  <p>In R:</p>
  
  <pre><code># setup:
  p1   = 5  # dimension of matrix to solve
  nrep1= 70 # number of repetitions
  set.seed(123)
  MAT0 =  matrix(rnorm(p1^2)/10,p1,p1)
  VEC0 =  matrix(rnorm(p1)/10,p1,1)
  
  # Runs okay:
  for(i in 1:10){
    r1=fuc(MAT0, VEC0, nrep1, p1)  ;str(r1)
  }
  
  
  # Results in crash of R (especially when executed serveral times):
  for(i in 1:10){
    r1=fuc(MAT0*0, VEC0, nrep1, p1)  ;str(r1)
  }
  </code></pre>
  
  <p>In Rstudio (R version 3.1.2 (2014-10-31, Platform: x86_64-w64-mingw32/x64 (64-bit))
  I typically receive a message like this:</p>
  
  <pre><code>Warning: stack imbalance in '.Call', 35 then 33
  Warning: stack imbalance in '=', 29 then 27
  Warning: stack imbalance in '{', 26 then 24
  </code></pre>
  
  <p>and then:</p>
  
  <pre><code>Error: Unable to establish connection with R session
  </code></pre>
  
  <p>and there is a fatal error.</p>