[R] Trouble returning 2D array into R from Fortran

[Berend Hasselman]

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.