32
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 ...
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 ...
24
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 ...
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.
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 ...
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 ...
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
Fortran's allocatable variables are automatically deallocated when the variable goes out of scope (see http://www.fortran90.org/src/best-practices.html#allocatable-arrays). This means that it is not possible to create a memory leak by failing to deallocate an allocatable array. This is one of the big benefits of using allocatable arrays rather than pointers.
...
10
Python is a very slow, high level language. For fast number crunching you'll have to write the main compute kernels in low level languages like C/C++ which means that now you have to learn not one but at least two languages. You'll also have to deal with additional headache associated with debugging/installation/maintenance etc. Most people use Python as a ...
9
I don't have experience with Octrees, but whenever there is some nice C++ library that I want to use in Fortran, I simply write a simple C driver --- typically a few C functions that do exactly what I need. Then I call them from Fortran using the iso_c_binding module. This has the great advantage that you reuse a well tested library with a community around ...
9
Java has been around for almost 20 years now as a major programming language, but it hasn't caught on in scientific computing so far. I think that's a good indicator for what's going to happen in the future.
My take is that the issue isn't speed. Most people are probably willing to give up 20% of performance (or even a factor of 2) if they would be vastly ...
9
You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq.
JuliaDiffEq and DifferentialEquations.jl has a large feature set dedicated to efficiently integrating computationally-difficult differential equations. It has a simple ...
8
I would stay away from Fortan, or if you must, use a reasonably new version (2003 rather than 77). A lot of physics software (Monte Carlo simulations in particular) is written in Fortran, simply because the projects were originally started in the 80s.
That being said, python and Fortran are two very different language, and what they should be used for is ...
8
There is no built-in Fortran functionality to do linear interpolation. You could either use a library or write your own routine. I haven't tried compiling or testing and my fortran may be a bit rusty, but something like the following should work.
subroutine interp1( xData, yData, xVal, yVal )
! Inputs: xData = a vector of the x-values of the data to be ...
8
I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS array, otherwise the interface will create an internal copy with the correct storage order.
Python 3.6.2 (default, Jul 22 2017, 21:19:22)
[GCC 7.1.1 20170516] on linux
Type "help", "copyright", "...
8
I think it's generally true that there are no advantages of Fortran 77 over either newer versions of Fortran or in fact any number of other programming languages that are widely used in scientific computing.
The reason it's still used is because there are millions of lines of code around that are written in Fortran 77. Now, recall that it takes a good ...
7
It is important to preallocate correctly. This is almost certainly the reason why your assembly was slow. If you are starting with a dense matrix representation, just scan through it once counting the number of nonzeros per row, then call MatSeqAIJSetPreallocation(). See this FAQ. The option MAT_IGNORE_ZERO_ENTRIES is really intended to be used when there is ...
7
It appears you are only taking half a time step on the fourth stage (where you should take a full one). Since you aren't taking the proper time step, this is turning into a complicated but relatively useless explicit scheme.
7
While not as free or available as GFortran that Bill mentions, the Intel Fortran Compiler also works F77 (or atleast works for the legacy code I work with). I am just putting this here to give an alternative, I would still recomment GFortran for most use cases simply due to the cost of ifort.
7
There are three issues that are likely to cause such problems in pseudospectral methods:
Gibbs oscillations
Aliasing
Time step too large
In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, ...
7
There is almost no good reason to write your own dense matrix manipulation routines until you have compared to fast libraries (MKL, FLAME, MAGMA, etc). Writing these libraries is challenging work that depends highly on the target architecture and requires deeper concepts than naive nested do loops. There's a whole literature on the subject that you should ...
7
One issue causing the jagged spectra at high wavenumbers is under sampling there. For example consider the 2D analogue of your binning procedure:
You don't want to sample from the red zones as they will become increasingly under-sampled as you move past a radius of size $|k_{x}|=|k_{y}|$.
6
I don't actually know how fortran compilers work, but based on the language features, I can take a guess.
Dynamic arrays in fortran come with meta-data to work with intrinsic functions like shape, size, lbound, ubound, and allocated or associated (allocatable vs pointers). For large arrays, the size of the meta-data is negligible, but for tiny arrays, like ...
6
What have you tried so far? This completely naive implementation manages to compute 7 (maybe 7.5) digits in 2.5 seconds on my laptop:
#include <iostream>
#include <complex>
#include <cmath>
#include <iomanip>
int main ()
{
const double alpha = 1;
std::cout.precision(16);
std::complex<double> sum = 0;
for (unsigned int ...
6
All you have to do is the following:
!alloc_test.f90
subroutine f(x, z, n)
implicit none
! Argument Declarations !
integer :: n
real*8, intent(in) :: x(n)
real*8, intent(out) :: z(n)
! Variable Declarations !
real*8, allocatable :: y(:)
! Variable Initializations !
allocate(y(n))
! Statements !
y(:) = 1.0
z = x + y
deallocate(y)
...
6
You seem to be very set on using Fortran. Octrees, when implemented efficiently, are rather complex data structures and, as such, better suited to programming languages that have more support for this, such as C/C++. There are a number of very high quality implementations in C/C++ that you could use.
6
You seem to have the wrong declaration for stat. It must be declared as an array of size MPI_STATUS_SIZE.
integer stat(MPI_STATUS_SIZE)
6
I think that the Metcalfe book is a good annotated guide to the modern features of Fortran, very useful if you already know an earlier version of the language but not a good introduction to programming in the language.
For my money the Chiver's and Sleightholme book and the Brainerd book are both good enough, though I don't think either is as good as the ...
6
This is probably not the answer you are looking for, but I wanted to state it anyway: Your choice of programming languages introduces two difficulties you will encounter as your program grows.
First, you will find that there is a cost to trying to couple different languages. They may interoperate, but there is always some quirk to it and as your software ...
Only top voted, non community-wiki answers of a minimum length are eligible
