[R] Trouble returning 2D array into R from Fortran
Berend Hasselman
bhh at xs4all.nl
Tue Oct 23 17:29:26 CEST 2012
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.
More information about the R-help
mailing list