[R] Trouble returning 2D array into R from Fortran
Berend Hasselman
bhh at xs4all.nl
Mon Oct 22 18:45:41 CEST 2012
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.
More information about the R-help
mailing list