Dear R experts
I work in computational fluid dynamics in 2D: I have a 200-by-200
array of fluid properties such as density and velocity and these
evolve in time (the precise equations depend on the problem).
Up to now, I've been using Fortran and the code is very very messy.
It works, but a professional programmer friend of mine saw the source
code once, and had to be strapped down for his own (and my!) safety.
It's time to rewrite it in a more object-oriented manner.
I have a couple of options: port the thing to c++, or deal directly in
R. The R route has a number of advantages. For example, I deal with
fluxes between adjacent chessboard-square-style elements. The Fortran
idiom is:
do 10 i=1,200
do 20 j=1,200
flux(i,j) = ( r(i,j)*u(i+1,j) + r(i+1,j)*u(i,j) ) /2
20 continue
10 continue
where r is the density and u the x-component of velocity. I might
need to do this or similar-looking things such as u(i+1,j)*v(i,j+1)
perhaps a hundred times in my Fortran code. But in R we could have
flux = ( r*right(u) + right(r)*u ) /2
[where right <- function(x){cbind(x[,-1],NA)} ]
To my mind, the functional form is much better: it's vectorized and
clear and terse. My question is, is it fast? (or more precisely, how
much slower would this nice approach be than my clunky old Fortran).
I guess I could tolerate a factor of two or three, and wait for a
shiny next-generation PC.
any comments anyone?
--
Robin Hankin, Lecturer,
School of Geography and Environmental Science
Tamaki Campus
Private Bag 92019 Auckland
New Zealand
r.hankin at auckland.ac.nz
tel 0064-9-373-7599 x6820; FAX 0064-9-373-7042
as of: Thu Oct 3 09:30:01 NZST 2002
This (linux) system up continuously for: 399 days, 16 hours, 12 minutes
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Robin Hankin <r.hankin at auckland.ac.nz> writes:> > do 10 i=1,200 > do 20 j=1,200 > flux(i,j) = ( r(i,j)*u(i+1,j) + r(i+1,j)*u(i,j) ) /2 > 20 continue > 10 continue > > where r is the density and u the x-component of velocity. I might > need to do this or similar-looking things such as u(i+1,j)*v(i,j+1) > perhaps a hundred times in my Fortran code. But in R we could have > > flux = ( r*right(u) + right(r)*u ) /2 > > [where right <- function(x){cbind(x[,-1],NA)} ] > > To my mind, the functional form is much better: it's vectorized and > clear and terse. My question is, is it fast? (or more precisely, how > much slower would this nice approach be than my clunky old Fortran). > I guess I could tolerate a factor of two or three, and wait for a > shiny next-generation PC. > > > any comments anyone?I would be highly surprised if it turned out to be even remotely fast in R... This kind of thing generally needs a compiled language like C(++) or Fortran. Preferably with good optimisers. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
[This was sent to Robin earlier and is posted here in case any others are interested in modern Fortran.] At 09:38 AM 10/03/2002 +1200, Robin Hankin wrote:>I work in computational fluid dynamics in 2D: I have a 200-by-200 >array of fluid properties such as density and velocity and these >evolve in time (the precise equations depend on the problem). >Up to now, I've been using Fortran and the code is very very messy.>[...]>flux = ( r*right(u) + right(r)*u ) /2 > >[where right <- function(x){cbind(x[,-1],NA)} ]Just to remind you, you can do this kind of thing in modern Fortran. Intrinsic vector and matrix operations have been available in the Fortran language for about ten years (since Fortran 90). In my experience, Fortran is usually much faster than R. It will probably be easier to port your code to modern Fortran than to C and its derivatives, and also less error prone. If the code is messy, why not clean it up? It is possible to write good code, and messy code, in any language. Some comparisons of Fortran and C++, and information on modern Fortran as a scientific programming language, can be found at http://www.cts.com.au/compare.html http://citeseer.nj.nec.com/cache/papers/cs/7409/http:zSzzSzwww.cs.colorado.eduzSz~zornzSzcs5535zSzFall-1995zSzprojectszSzwharton.pdf/wharton95should.pdf http://www.kcl.ac.uk/kis/support/cit/fortran/f90home.html http://www.ibiblio.org/pub/languages/fortran/ch1-2.html Whatever path you choose, good luck! -- Michael Prager <Mike.Prager at noaa.gov> NOAA Center for Coastal Fisheries and Habitat Research Beaufort, North Carolina 28516 USA http://shrimp.ccfhrb.noaa.gov/~mprager/ Standard Disclaimers: * Opinions expressed here are personal and are not otherwise represented. * Any use of tradenames does not constitute a NOAA or NMFS endorsement. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._