Fortran is really a nice language when it comes to scientific computations. So nice that I have some difficulties when using a different one like Python.
If you do science with Python, you always end up with using numpy, a wonderful extension to work with arrays and matrices efficiently. But... there is a but... it is not statically typed and can create some funny results, like when you multiply a matrix by a vector with a dot product.
>>> import numpy as np
>>> A = np.asarray([[1.0, 2.0], [3.0, 4.0]])
>>> B = np.asarray([5.0, 6.0])
>>> np.dot(A, B)
array([ 17., 39.])
>>> np.dot(B, A)
array([ 23., 34.])
What is happening here is that when you run np.dot(A, B)
, as B
is
a simple array, it is replicated two times to match the size of A
,
so you effectively run: 5x1 + 6x2
and again 5x3 + 6x4
.
When you run np.dot(B, A)
, you run in columns: 5x1 + 6x3
and
5x2 + 6x4
.
Numpy is smart, but deep into your computational code, an order swap or an error in dimensions can result in really hard to track errors. Fortran will simply send you back to the drawing board telling you that the dimensions do not match and that you should work a bit more.