Trouble returning 2D array into R from Fortran
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.