Recently I am trying to add a Fortran function into an existing R package that contains C++ code and is build under Rcpp package. I have successfully added the Fortran function and build the R package. But when I try to run an example for the R package for different times, each time the function returns a different value, which a bit confusing since the Fortran function I added has no undeterministic parameter or contains self iteration. Also sometimes when I try to run the function, R crashes.
This is a part of the Fortran function that are in the file "src/k2s.f95"
REAL*8 FUNCTION k2(n1,n2,n3,vec1,length1,a1,vec2,length2,a2)
integer :: n1,n2,n3,length1,length2
integer, dimension(length1) :: vec1
double precision :: a1,a2
double precision, dimension(length2):: vec2
double precision :: P(n1+2)
...
k2 = -1.0
IF ((n1<1) .OR. (n2<1) .OR. (length2 .NE. n1+n2-1) .OR. (n3<1) .OR. (n3 > 3)) RETURN
k2 = -2.0
IF(MINVAL(vec1).LE.0) RETURN
IF(SUM(vec1).NE.(n1+n2)) RETURN
IF(MINVAL(vec2).LE.0) RETURN
...
END FUNCTION K2
And I have used this answer Integrate Fortran, C++ with R to construct the package. Therefore I have written the following C++ file under src dictionary. This is the content of C++ file "k2s2.cpp"
#include "Rcpp.h"
extern "C" {
double k2_(int *n1, int *n2, int *n3, int vec1[], int *length1, double *a1, double vec2[], int *length2, double *a2);
}
// [[Rcpp::export]]
double K2_fortran(int n1, int n2, int n3, Rcpp::IntegerVector vec1, double a1, Rcpp::NumericVector vec2, double a2)
{
int length1 = vec1.size();
int length2 = vec2.size();
double q = 0;
pval = k2_(&n1,&n2,&n3,vec1.begin(),&length1,&a1,vec2.begin(),&length2,&a2);
return q;
}
I think it should be fine. Then I build the package (with Rcpp automatically generates the file "src/RcppExports.cpp" and "R/RcppExports.R" to link the cpp function K2_fortran to R) and load all the functions, then I get the R function K2_fortran. And I run the following example.
K2_fortran(120, 150, 1, c(80,70,40,80), 0.1, rep(1,269), 1e-6)
And the result is -2. Then I rerun the code for the second time, I got the expected result 0.175. For the third time, I got 1(some other returned value when there is an error). Then I run it again, a fatal error occurs in R and closes the R.
If I reopen it and run this example 4 times again. The result will be again -2, 0.175, 1 and a fatal error. I am not sure what has happened.
