rrecursioninlinercppstrassen

calling a cxxfunction itself recursion


I have a problem with my cxxfunction in R. And I would like to call this itself, unfortunately, the compiler gives me this error message:

'Matmult2' is not declared in this scope

The Problem is also to call Matmult2 in a cxxfunction.

My original problem is related to the strassen algorithm, I would like to call recursively with RCPP / Inline. Many thanks for your help!

`

Matmult2<- cxxfunction(signature( a="NumericMatrix",b="NumericMatrix"), plugin="RcppArmadillo",includes = rcpp_inc,
  body='

              Rcpp::NumericMatrix verg(a);
              int n = verg.nrow();

              arma::mat A = Rcpp::as< arma::mat >(a);
              arma::mat B = Rcpp::as< arma::mat >(b);

              arma::mat a11 = A(span(0,(n/2)-1), span(0,(n/2)-1));
              arma::mat a12 = A(span(0,(n/2)-1), span((n/2),n-1));
              arma::mat a21 = A(span((n/2),n-1), span(0,(n/2)-1));
              arma::mat a22 = A(span((n/2),n-1), span((n/2),n-1));

              arma::mat b11 = B(span(0,(n/2)-1), span(0,(n/2)-1));
              arma::mat b12 = B(span(0,(n/2)-1), span((n/2),n-1));
              arma::mat b21 = B(span((n/2),n-1), span(0,(n/2)-1));
              arma::mat b22 = B(span((n/2),n-1), span((n/2),n-1));


              if ( n  < 780) return Rcpp::wrap( A * B );


              arma::mat P1(Matmult2(a11,(b12-b22)));
              arma::mat P2(Matmult2((a11+a12),(b22)));
              arma::mat P3(Matmult2((a21+a22),(b11)));
              arma::mat P4(Matmult2((a22),(b21-b11)));
              arma::mat P5(Matmult2((a11+a22),(b11+b22)));
              arma::mat P6(Matmult2((a12-a22),(b21+b22)));
              arma::mat P7(Matmult2((a11-a21),(b11+b12)));

              arma::mat c11((P5+P4-P2+P6));
              arma::mat c12((P1+P2));
              arma::mat c21((P3+P4));
              arma::mat c22((P5+P1-P3-P7));

              mat C1 = join_rows(c11, c12);
              mat C2 = join_rows(c21, c22);
              mat C = join_cols(C1, C2);

              return Rcpp::wrap( C );

')

`


Solution

  • You can't do direct recursion because the glue code we insert creates a wrapper function that is called---but you do not want the wrapper to be called recursively, you want your function to be called. So you would need too: a wrapper called from R, and a C(++)-only recursive function that is called from your wrapper, and which calls itself.

    It also works if you switch to Rcpp Attributes. Here we use the standard bearer of recursive functions, the Fibonacci sequence:

    R> library(Rcpp)
    R> cppFunction("double fib(double n) { if (n < 2) return n; return(fib(n-1) + fib(n-2)); }")
    R> sapply(0:10, fib)
     [1]  0  1  1  2  3  5  8 13 21 34 55
    R> 
    

    (And I use double instead of int as the latter overflows earlier.)