62
As so often, the choice depends on (1) the problem you are trying to solve, (2) the skills you have, and (3) the people you work with (unless it's a solo project). I'll leave (3) aside for the moment because it depends on everyone's individual situation.
Problem dependence: Fortran excels at array processing. If your problem can be described in terms of ...
50
Language designers face many choices. Ken Kennedy emphasized two: (1) better abstractions and (2) higher- or lower-level (less or more machine-like) code. While functional languages like Haskell and Scheme focus on the former, traditional scientific-computing languages like Fortran and C/C++ focused on the latter. Saying that one language is faster than ...
37
I think that both C++ and Fortran are good enough and work well.
However I think that Fortran is better for numeric scientific computing, for algorithms that can be expressed using arrays and don't need other sophisticated data structures, so in fields like finite differences/elements, PDE solvers, electronic structure calculations. Fortran is a domain ...
37
The difference in your timings seems to be due to the manual unrolling of the unit-stride Fortran daxpy. The following timings are on a 2.67 GHz Xeon X5650, using the command
./test 1000000 10000
Intel 11.1 compilers
Fortran with manual unrolling: 8.7 sec
Fortran w/o manual unrolling: 5.8 sec
C w/o manual unrolling: 5.8 sec
GNU 4.1.2 compilers
Fortran ...
33
Pick your poison. I recommend using Homebrew. I have tried all of these methods except for "Fink" and "Other Methods". Originally, I preferred MacPorts when I wrote this answer. In the two years since, Homebrew has grown a lot as a project and has proved more maintainable than MacPorts, which can require a lot of PATH hacking.
Installing a version that ...
31
I'm also throwing my two cents in kind of late, but I've only just seen this thread and I feel that, for posterity, there are a few points that desperately need to be made.
Note in the following that I will talk about C and not C++. Why? Well, otherwise it's apples and oranges to compare a full-fledged dynamically typed object-oriented language with ...
29
Ease of learning
Python and Fortran are both relatively easy-to-learn languages. It's probably easier to find good Python learning materials than good Fortran learning materials because Python is used more widely, and Fortran is currently considered a "specialty" language for numerical computing.
I believe the transition from Python to Fortran would be ...
27
It's a bit of a popular misnomer that there is a "version" of Fortran to know. With rare exception, the latest Fortran standards (and compilers) retain excellent backwards compatibility with older standards. This is with good reason: not many people would use Fortran today if it weren't for the large amounts of legacy code still in use. That is to say, a ...
23
The design of Fortran allows the compiler to perform stronger optimizations in some cases, optimizations that are not generally available to C.
One famous example is the handling of aliasing. In Fortran, you can access a specific memory area only though the specific symbol associated with that memory area. This knowledge allows the compiler to employ smart ...
23
I am not super familiar with f2py internals, but I am very familiar with wrapping Fortran. F2py just automates some or all of the things below.
You first need to export to C using the iso_c_binding module, as described for example here:
http://fortran90.org/src/best-practices.html#interfacing-with-c
Disclaimer: I am the main author of the fortran90.org ...
20
Another way to do this is to first explicitly specify the precision you desire in the variable using the SELECTED_REAL_KIND intrinsic and then use this to define and initialize the variables. Something like:
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
REAL(dp) :: x
x = 1.0_dp
A nice advantage to doing it this way is that you can store the definition ...
19
You can use the Python builtin ctypes module as described on fortran90.org. It is pretty straight forward and doesn't require any external dependencies. Also, the ndpointer arg type helper is very handy.
19
Let me try and break down your requirements:
Maintainability
Reading/writing text data
Strong interfaces/capability for LU factorizations
Sparse linear solvers
Performance and scalability to large data
From this list, I would consider the following languages:
C, C++, Fortran, Python, MATLAB, Java
Julia is a promising new language, but the community is ...
18
The column major layout is the scheme used by Fortran and that's why it's used in LAPACK and other libraries.
In general it is much more efficient in terms of memory bandwidth usage and cache performance to access the elements of an array in the order in which they're laid out in memory. Depending on how your matrices are stored, you'll want to pick ...
18
GFortran is free and open source and usually works pretty well.
16
I'm coming late to this party, so it's hard for me to follow the back-and-forth from all above. The question is big, and I think if you are interested it could be broken up into smaller pieces. One thing I got interested in was simply the performance of your daxpy variants, and whether Fortran is slower than C on this very simple code.
Running both on my ...
16
From my 15 years of thinking about scientific software: If your code runs 25% faster because you write it in Fortran, but it takes you 4 times as long to write it (no STL, difficulty implementing complex data structures, etc), then Fortran only wins if you spend a significant fraction of your day twiddling thumbs and waiting for your computations to finish. ...
16
Matrix exponentials of skew-Hermitian matrices are cheap to compute:
Suppose $A$ is your skew-Hermitian matrix, then $iA$ is Hermitian, and via zheevd and friends you can get the decomposition
$$iA = U \Lambda U^H,$$
where $U$ is the unitary eigenvector matrix and $\Lambda$ is real and diagonal. Then, trivially,
$$A = U (-i \Lambda) U^H.$$
Once you ...
14
To make more robust comparisons (on linux), you can :
1) On Intel CPUs the turbo overclocks your CPU. This is controlled by the temperature of the CPU, so it can behave differently from one run to the other. On Linux, you can block the frequency of the CPU as follows. For example, for 2.4GHz:
echo 1 > /sys/module/processor/parameters/ignore_ppc
for ...
13
My approach has been to use C++ for everything but computational kernels, which are usually best written in assembly; this buys you all of the performance of the traditional HPC approach but allows you to simplify the interface, e.g., by overloading computational kernels like SGEMM/DGEMM/CGEMM/ZGEMM into a single routine, say Gemm. Clearly the abstraction ...
13
You look at how deal.II (http://www.dealii.org/) does it -- there, dimension independence lies at the very heart of the library, and is modeled as a template argument to most data types. See, for example, the dimension-agnostic Laplace solver in the step-4 tutorial program:
http://www.dealii.org/developer/doxygen/deal.II/step_4.html
See also
https://...
12
You declared the variables as double precision, but you initialized them with single precision values.
You could have written:
X=1.0d0
Y=1.0d-1
Barron's answer below is another way of making a literal double precision, with the advantage that it allows you to change the precision of your variables at a later time.
12
I don't think Fortran is that close to the metal (see other answer) but it tends to optimize very easily. Loops are simple, and the language readily supports vectorization extensions (okay when I used it in my first job we were targeting a wide range of vector big iron).
There is also the large factor of inertia. A lot of numeric code is in Fortran, so ...
12
The question highlights that most "plain" programming languages (C, Fortran, at least) do not allow you to do this cleanly. An added constraint is that you want notational convenience and good performance.
Therefore, instead of writing a dimension-specific code, consider writing a code that generates a dimension-specific code. This generator is dimension-...
12
First of all, thanks for posting this question/challenge! As a disclaimer, I'm a native C programmer with some Fortran experience, and feel most at home in C, so as such, I will focus only on improving the C version. I invite all Fortran hacks to have their go too!
Just to remind newcomers about what this is about: The basic premise in this thread was that ...
11
In vacuum without considering any existing software, there's no reason to prefer column major over row major from the code point of view. However, most mathematical literature is written in a way that groups vectors into a matrix by storing them as columns instead of rows. For example when you write the full eigenvalue equation $AX=X\Lambda$, the $X$ matrix ...
11
There are a few other ways to do this, with advantages and drawbacks:
MPI_WTIME: This is a high resolution wall-clock. It is probably the most `trusted' option; it just works. The downside is that if your program doesn't already use MPI, you'll have to wrap MPI around it (which isn't hard).
Use a fortran intrinsic (as you have): This is probably the ...
10
There are modern CFD codes which you can look into. For example ...
Fluidity: General purpose multiphase CFD (FE) code; Even does fully unstructured AMR
WRF: Next-generation numerical model weather prediction system from NCAR
Code Saturne: General purpose CFD (FV) code; Some features listed on Wikipedia
There are similar modern Fortran codes for ...
10
With the Intel compiler of any modern vintage, -O3 -vec-report3. Optimization level three guarantees that it's trying to vectorize, and the vector report will tell you what it's doing.
The GNU page on vectorization says that it's on by default at optimization level 3, but I can't find the equivalent of vec-report.
Only top voted, non community-wiki answers of a minimum length are eligible
