Hello,
I have been trying to use a collection of Fortran subroutines to return a 2D
array of calculated values to my R code, calling a Fortran wrapper
subroutine from R. I've done this successfully before with C & C++ code.
The Fortran wrapper subroutine which is to be called by R takes a set of
input arguments & then should return an array of 2 columns & the specified
number of rows.
I've tested the wrapping subroutine from another Fortan main program & it
does work as expected, so the Fortran works, but my problem has been with
the correct syntax for retrieving the output double array from within R.
The wrapping Fortran subroutine is defined;
subroutine xypos_parallax_r(k,year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,y)
k is the number of rows & y an assumed size double array, y(2,*), the other
arguments are other input variables for the Fortran.
I've tried putting the Fortran function call inside an R function, with this
syntax;
#Source positions in lens plane under annual parallax
xypos_parallax <- function(year,ra,dec,ti,model_par) {
if (!is.loaded('xypos_parallax_r')){
dyn.load('parallax.so')}
returned_data = .Fortran('xypos_parallax_r',
as.integer(length(ti)),
as.integer(year),
as.double(ra),
as.double(dec),
as.double(ti),
as.double(model_par$t0),
as.double(model_par$tE),
as.double(model_par$alpha),
as.double(model_par$u0),
as.double(model_par$piee),
as.double(model_par$pien),
as.double(array("double",dim=c(2,length(ti)))) )[[12]]
}
The input arguments to the Fortran call include the number of rows of the
output array, with the final argument in the .Fortran call intended to be
the output array itself. ti is a vector of same length as the number of rows
needed in the output array.
Testing this R function calling the Fortran code returns;
Error in xypos_parallax(year, ra, dec, ti, model_par) :
NA/NaN/Inf in foreign function call (arg 12)
In addition: Warning message:
In xypos_parallax(year, ra, dec, ti, model_par) :
NAs introduced by coercion
The Fortran code does work as intended, so the problem must be with how I've
written the R function making the Fortran call, but I can't see where I've
gone wrong.
Would anyone have any idea, or experience with returning multi-dimensional
arrays from external code into R?
Thanks in advance
--
View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862.html
Sent from the R help mailing list archive at Nabble.com.
Trouble returning 2D array into R from Fortran
9 messages · paulfjbrowne, Berend Hasselman
See inline.
On 20-10-2012, at 17:18, paulfjbrowne <paulfj.browne at gmail.com> wrote:
Hello,
I have been trying to use a collection of Fortran subroutines to return a 2D
array of calculated values to my R code, calling a Fortran wrapper
subroutine from R. I've done this successfully before with C & C++ code.
The Fortran wrapper subroutine which is to be called by R takes a set of
input arguments & then should return an array of 2 columns & the specified
number of rows.
I've tested the wrapping subroutine from another Fortan main program & it
does work as expected, so the Fortran works, but my problem has been with
the correct syntax for retrieving the output double array from within R.
The wrapping Fortran subroutine is defined;
subroutine xypos_parallax_r(k,year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,y)
k is the number of rows & y an assumed size double array, y(2,*), the other
arguments are other input variables for the Fortran.
I've tried putting the Fortran function call inside an R function, with this
syntax;
#Source positions in lens plane under annual parallax
xypos_parallax <- function(year,ra,dec,ti,model_par) {
if (!is.loaded('xypos_parallax_r')){
dyn.load('parallax.so')}
returned_data = .Fortran('xypos_parallax_r',
as.integer(length(ti)),
as.integer(year),
as.double(ra),
as.double(dec),
as.double(ti),
as.double(model_par$t0),
as.double(model_par$tE),
as.double(model_par$alpha),
as.double(model_par$u0),
as.double(model_par$piee),
as.double(model_par$pien),
as.double(array("double",dim=c(2,length(ti)))) )[[12]]
}
Why as.double(array("double",dim=c(2,length(ti))?
You are passing a character array to your fortran subroutine.
Use matrix(0,nrow=2,ncol=length(ti)) which will be a matrix of double precision's.
Berend
The input arguments to the Fortran call include the number of rows of the output array, with the final argument in the .Fortran call intended to be the output array itself. ti is a vector of same length as the number of rows needed in the output array. Testing this R function calling the Fortran code returns; Error in xypos_parallax(year, ra, dec, ti, model_par) : NA/NaN/Inf in foreign function call (arg 12) In addition: Warning message: In xypos_parallax(year, ra, dec, ti, model_par) : NAs introduced by coercion The Fortran code does work as intended, so the problem must be with how I've written the R function making the Fortran call, but I can't see where I've gone wrong. Would anyone have any idea, or experience with returning multi-dimensional arrays from external code into R? Thanks in advance -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
I will look into using inline, but since the Fortran code is several thousand lines long & is comprised of multiple subroutines, compiling it into a shared object & dynamically loading it into R is probably the easier solution. I have also noticed a strange numerical problem when calling the routine from within R as opposed to in its native Fortran. For example, for the same set of input parameters, R will output; (0.0315031507927081, 0.220339391742628) Whereas the Fortran code internally outputs; (0.0350479965640488472, 0.22220883087569138) All variables are defined in the Fortran code as double, and are passed in as double in the .Fortran call in R. I am a bit puzzled at the discrepancy between the two calls' output. -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4646874.html Sent from the R help mailing list archive at Nabble.com.
Look at my comments in between your post.
On 20-10-2012, at 19:18, paulfjbrowne <paulfj.browne at gmail.com> wrote:
I will look into using inline, but since the Fortran code is several thousand lines long & is comprised of multiple subroutines, compiling it into a shared object & dynamically loading it into R is probably the easier solution.
You misread my "see inline". I did not mean look at the inline package. I meant that my comments were inline (I try never to top post).
I have also noticed a strange numerical problem when calling the routine from within R as opposed to in its native Fortran. For example, for the same set of input parameters, R will output; (0.0315031507927081, 0.220339391742628) Whereas the Fortran code internally outputs; (0.0350479965640488472, 0.22220883087569138) All variables are defined in the Fortran code as double, and are passed in as double in the .Fortran call in R. I am a bit puzzled at the discrepancy between the two calls' output.
Without further information and code, it is impossible to give any kind of sensible advice. In the Fortran code you could try to check the input arguments. Are you sure that from R you are passing exactly the same values as "in its native Fortran". Berend
-- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4646874.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
1 day later
I haven't determined the root cause of the numerical accuracy bug yet (I haven't determined a way to compare the values used when the code is run from R as opposed to from a Fortran test), but I did have another issue with retrieving the matrix from the Fortran code that I hoped you might be able to help with. I'm very much a beginner user of R, so any help would be appreciated. The attached R file contains the R function that calls the Fortran subroutines & code to make a test run. The attached Fortran file is the Fortran wrapper subroutine called by R & a second subroutine that it then calls to calculate the two elements of each row of the output matrix. When running this code, the matrix returned by the R function is incorrect in that the returned elements are being stored sequentially down the first row & then the second row of the 2D output matrix, rather than sequentially across the two columns as desired. Might you be able to point me to where I've written the code wrong so that the output is being stored incorrectly? Thanks in advance for any help parallax_Rhelp.f90 <http://r.789695.n4.nabble.com/file/n4647005/parallax_Rhelp.f90> xypos_parallax.R <http://r.789695.n4.nabble.com/file/n4647005/xypos_parallax.R> -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647005.html Sent from the R help mailing list archive at Nabble.com.
Your attachments did not make it through to the official R-help mailing list. Nabble R is not R-help. I downloaded the files from the links to Nabble. See further down for my comments.
On 22-10-2012, at 14:24, paulfjbrowne wrote:
I haven't determined the root cause of the numerical accuracy bug yet (I haven't determined a way to compare the values used when the code is run from R as opposed to from a Fortran test), but I did have another issue with retrieving the matrix from the Fortran code that I hoped you might be able to help with. I'm very much a beginner user of R, so any help would be appreciated. The attached R file contains the R function that calls the Fortran subroutines & code to make a test run. The attached Fortran file is the Fortran wrapper subroutine called by R & a second subroutine that it then calls to calculate the two elements of each row of the output matrix.
Butthe code for the geta subroutine is missing. I commented out the if branch using the geta subroutine. So all calculations are for dtau=0 and du0=0.
When running this code, the matrix returned by the R function is incorrect in that the returned elements are being stored sequentially down the first row & then the second row of the 2D output matrix, rather than sequentially across the two columns as desired.
I do not understand what you mean by this. Your fortran routine declares a Nx2 matrix i.e. a matrix with 2 columns. The code do i=1,k call xypos_parallax(year,ra,dec,ti(i),t0,tE,alpha,u0,piee,pien,y1,y2) y(i,1) = y1 y(i,2) = y2 end do fills the columns in row i sequentially from row 1 to row k. Fortran stores arrays column-wise. In the R function calling your Fortran you declare the output matrix as a matrix with 2 columns. R also stores matrix in column order. So you are getting exactly what you have programmed. Running your code as follows for a short input vector, this is the result (with geta disabled):
ti<-seq(from=5767.0299999999997,to=5797.0,length=10) ra<-264.5590833 dec<--27.13613889 year<-2011 y<-xypos_parallax(year,ra,dec,ti,model_par) y
[,1] [,2] [1,] -0.118969682 -0.203707575 [2,] -0.103819120 -0.153139480 [3,] -0.088668558 -0.102571385 [4,] -0.073517996 -0.052003290 [5,] -0.058367434 -0.001435194 [6,] -0.043216872 0.049132901 [7,] -0.028066310 0.099700996 [8,] -0.012915748 0.150269091 [9,] 0.002234814 0.200837187 [10,] 0.017385376 0.251405282 Would you by any chance mean transposed output like this
t(y)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.1189697 -0.1038191 -0.08866856 -0.07351800 -0.058367434 -0.04321687
[2,] -0.2037076 -0.1531395 -0.10257138 -0.05200329 -0.001435194 0.04913290
[,7] [,8] [,9] [,10]
[1,] -0.02806631 -0.01291575 0.002234814 0.01738538
[2,] 0.09970100 0.15026909 0.200837187 0.25140528
Might you be able to point me to where I've written the code wrong so that the output is being stored incorrectly?
As I said above the output is not being stored incorrectly. Berend
Thanks in advance for any help parallax_Rhelp.f90 <http://r.789695.n4.nabble.com/file/n4647005/parallax_Rhelp.f90> xypos_parallax.R <http://r.789695.n4.nabble.com/file/n4647005/xypos_parallax.R> -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647005.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Yes, I had stripped out the Fortran from the geta subroutine downwards because is wasn't very important in returning a test output matrix. I think that I should make explicit what problem I believe I am seeing in what I have coded, and what form of the output matrix I would prefer. When I run the R code
#Source positions in lens plane under annual parallax
xypos_parallax <- function(year,ra,dec,ti,model_par) {
if (!is.loaded('xypos_parallax_r')){
dyn.load('parallax.so')}
returned_data = .Fortran('xypos_parallax_r',
as.integer(length(ti)),
as.integer(year),
as.double(ra),
as.double(dec),
as.double(ti),
as.double(model_par$t0),
as.double(model_par$tE),
as.double(model_par$alpha),
as.double(model_par$u0),
as.double(model_par$piee),
as.double(model_par$pien),
matrix(0,nrow=length(ti),ncol=2) )[[12]]
}
model_par<-list(t0=5781.49344,
tE=63.0814,
alpha=1.2797,
u0=0.0555,
rho=0.0063,
d=0.7119,
q=0.0020496,
piee=-0.900,
pien=-0.800)
ti<-seq(from=5770.0,to=5797.0,length=10)
ra<-264.5590833
dec<--27.13613889
year<-2011
y<-xypos_parallax(year,ra,dec,ti,model_par)
calling the Fortran subroutines;
!************************************************************************************************** !XYPOS_PARALLAX_R: R wrapper subroutine !************************************************************************************************** subroutine xypos_parallax_r(k,year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,y) implicit none integer::year,i,k double precision,dimension(k)::ti double precision::ra,dec double precision::t0,tE,alpha,u0 double precision::piee,pien double precision::y1,y2 double precision,dimension(2,k)::y do i=1,k call xypos_parallax(year,ra,dec,ti(i),t0,tE,alpha,u0,piee,pien,y1,y2) y(1,i) = y1 y(2,i) = y2 end do return end !************************************************************************************************** !XYPOS_PARALLAX: Source positions under annual parallax !************************************************************************************************** subroutine xypos_parallax(year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,xx,yy) implicit none double precision::ra,dec double precision::ti,t0,tE,alpha,u0 double precision::piee,pien double precision::xx,yy integer::year double precision::qn, qe, thjd, t0dum double precision::dtau,du0 thjd = ti t0dum = t0 IF((piee .ne. 0) .or. (pien .ne. 0))then call geta(qn,qe,thjd,ra,dec,t0dum,year) !the following equations are still to be checked. dtau = piee*qn + pien*qe du0 = -piee*qe + pien*qn du0 = -du0 ELSE dtau=0 du0=0 END IF xx = ((ti-t0)/te+dtau)*cos(alpha)-(u0+du0)*sin(alpha) yy = ((ti-t0)/te+dtau)*sin(alpha)+(u0+du0)*cos(alpha) return end
The matrix y which R returns is;
-
*
0.09712576
*
-0.036362910
*
-0.17291343
*
0.067752761 -0.08716623 -0.020590942 -0.12098203 0.109764867 -0.07617879 -0.003468398 -0.07086033 0.149555601 -0.06409356 0.015060440 -0.02263996 0.187057573 -0.05084316 0.035047997 0.02359192 0.222208831
But from Fortan, via
do i=1,npoint write(*,*) y(1,i),y(2,i) end do
The following matrix is reported;
*
-9.71257556964655877E-002 -0.17291342658178052
*
-8.71662283852910891E-002 -0.12098202580263227 -7.61787919231983052E-002 -7.08603339220746364E-002 -6.40935606206668174E-002 -2.26399573881042940E-002 -5.08431642024918945E-002 2.35919215412063593E-002 -3.63629100550238102E-002 6.77527605328501203E-002 -2.05909424429263285E-002 0.10976486692565257 -3.46839842574317436E-003 0.14955560107302845 1.50604398227174099E-002 0.18705757282679777 3.50479966174077409E-002 0.22220883066716274
The second is the form of the matrix that should be stored in the returned matrix y in the R function call, if possible. Do you see how the elements are being stored in the wrong order in the matrix that R is storing, versus the one that is read out in the Fortran? Is it possible to get the form of the second matrix? Thanks in advance. -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647106.html Sent from the R help mailing list archive at Nabble.com.
On 23-10-2012, at 12:33, paulfjbrowne wrote:
Yes, I had stripped out the Fortran from the geta subroutine downwards because is wasn't very important in returning a test output matrix.
But then any results are not reproducible. So there is no purpose in trying to reproduce your results. See inline for further comments.
I think that I should make explicit what problem I believe I am seeing in what I have coded, and what form of the output matrix I would prefer. When I run the R code
#Source positions in lens plane under annual parallax
xypos_parallax <- function(year,ra,dec,ti,model_par) {
if (!is.loaded('xypos_parallax_r')){
dyn.load('parallax.so')}
returned_data = .Fortran('xypos_parallax_r',
as.integer(length(ti)),
as.integer(year),
as.double(ra),
as.double(dec),
as.double(ti),
as.double(model_par$t0),
as.double(model_par$tE),
as.double(model_par$alpha),
as.double(model_par$u0),
as.double(model_par$piee),
as.double(model_par$pien),
matrix(0,nrow=length(ti),ncol=2) )[[12]]
}
Output matrix defined with dimensions length(ti) rows and 2 columns!!! This conflicts with the use in the Fortran code. See further on.
model_par<-list(t0=5781.49344,
tE=63.0814,
alpha=1.2797,
u0=0.0555,
rho=0.0063,
d=0.7119,
q=0.0020496,
piee=-0.900,
pien=-0.800)
ti<-seq(from=5770.0,to=5797.0,length=10)
ra<-264.5590833
dec<--27.13613889
year<-2011
y<-xypos_parallax(year,ra,dec,ti,model_par)
calling the Fortran subroutines;
!************************************************************************************************** !XYPOS_PARALLAX_R: R wrapper subroutine !************************************************************************************************** subroutine xypos_parallax_r(k,year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,y) implicit none integer::year,i,k double precision,dimension(k)::ti double precision::ra,dec double precision::t0,tE,alpha,u0 double precision::piee,pien double precision::y1,y2 double precision,dimension(2,k)::y
You are now changing the rules of the game. Array y is now 2 rows and k columns. And k == length(ti). In your previous post, which I used, y was an array of k rows and 2 columns!
do i=1,k call xypos_parallax(year,ra,dec,ti(i),t0,tE,alpha,u0,piee,pien,y1,y2) y(1,i) = y1 y(2,i) = y2 end do
In your previous post this loop was do i=1,k call xypos_parallax(year,ra,dec,ti(i),t0,tE,alpha,u0,piee,pien,y1,y2) y(i,1) = y1 y(i,2) = y2 So you seem to have changed things. But how can you expect to get interpretable results if you let R know that your matrix is Nx2 when in the Fortran code you declare the array as 2xN? Since the total number of elements is 2N in both cases you are not getting into memory access trouble.
return end !************************************************************************************************** !XYPOS_PARALLAX: Source positions under annual parallax !************************************************************************************************** subroutine xypos_parallax(year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,xx,yy) implicit none double precision::ra,dec double precision::ti,t0,tE,alpha,u0 double precision::piee,pien double precision::xx,yy integer::year double precision::qn, qe, thjd, t0dum double precision::dtau,du0 thjd = ti t0dum = t0 IF((piee .ne. 0) .or. (pien .ne. 0))then call geta(qn,qe,thjd,ra,dec,t0dum,year) !the following equations are still to be checked. dtau = piee*qn + pien*qe du0 = -piee*qe + pien*qn du0 = -du0 ELSE dtau=0 du0=0 END IF xx = ((ti-t0)/te+dtau)*cos(alpha)-(u0+du0)*sin(alpha) yy = ((ti-t0)/te+dtau)*sin(alpha)+(u0+du0)*cos(alpha) return end
The matrix y which R returns is;
-
*
0.09712576
*
-0.036362910
*
-0.17291343
*
0.067752761 -0.08716623 -0.020590942 -0.12098203 0.109764867 -0.07617879 -0.003468398 -0.07086033 0.149555601 -0.06409356 0.015060440 -0.02263996 0.187057573 -0.05084316 0.035047997 0.02359192 0.222208831
But from Fortan, via
do i=1,npoint write(*,*) y(1,i),y(2,i) end do
So in Fortran you want an array of N rows and 2 columns printed as column 1 and 2 of row 1 column 1 and 2 of row 2 .... In other words you want the output matrix as defined in your Fortran program to be printed in Transposed form. You should get what you want if in the R calling function you declare the matrix correctly. I.e. length(ti) rows and 2 columns. R is not storing anything in the wrong order; it is taking what the Fortran code is delivering. R receives a matrix and you are telling R the dimensions. Thereafter R is doing what you tell it to do. Use t(y) which should display what you want. Berend
The following matrix is reported;
*
-9.71257556964655877E-002 -0.17291342658178052
*
-8.71662283852910891E-002 -0.12098202580263227 -7.61787919231983052E-002 -7.08603339220746364E-002 -6.40935606206668174E-002 -2.26399573881042940E-002 -5.08431642024918945E-002 2.35919215412063593E-002 -3.63629100550238102E-002 6.77527605328501203E-002 -2.05909424429263285E-002 0.10976486692565257 -3.46839842574317436E-003 0.14955560107302845 1.50604398227174099E-002 0.18705757282679777 3.50479966174077409E-002 0.22220883066716274
The second is the form of the matrix that should be stored in the returned matrix y in the R function call, if possible. Do you see how the elements are being stored in the wrong order in the matrix that R is storing, versus the one that is read out in the Fortran? Is it possible to get the form of the second matrix? Thanks in advance. -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647106.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Thank you for your help. In any case, I eventually realized that I had reversed the Fortran/R column major order in the Fortran, which is why elements were returned incorrectly. As you point out, I did not run into memory access issues because of the self-similar array sizes. -- View this message in context: http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647167.html Sent from the R help mailing list archive at Nabble.com.