Column major versus Row major¶
Computer memory is a flat, 1D line of addresses. When a 2D grid/matrix is saved into memory, it must be flattened into a single sequence of numbers.
The difference between Fortran and Python/C/C++ comes down to whether that sequence goes column-by-column or row-by-row .
Visual Comparison¶
Consider this 2D matrix:
| Language | Order Type | Memory Storage Layout |
|---|---|---|
| Fortran | Column-Major | [1, 3, 2, 4] (Column 1 first, then Column 2) |
| C / C++ / Python | Row-Major | [1, 2, 3, 4] (Row 1 first, then Row 2) |
Key Differences¶
1. Which Index Changes Fastest in Memory?¶
- Fortran (Column-Major): The first index changes fastest. Moving from
A(1, 1)toA(2, 1)jumps to the very next slot in memory. - C / C++ / Python (Row-Major): The last index changes fastest. Moving from
A[0][0]toA[0][1]jumps to the very next slot in memory.
2. Writing Efficient Loops¶
To process arrays as fast as possible, code should step through memory sequentially (stride-1 access). This means the innermost loop must control the index that changes fastest in memory:
Fortran (Loop columns on the outside, rows on the inside):
Fortran
! Correct for Fortran: 'i' (first index) is in the inner loop
do j = 1, cols
do i = 1, rows
A(i, j) = A(i, j) + 1.0
end do
end do
C / C++ (Loop rows on the outside, columns on the inside):
C++
// Correct for C/C++: 'j' (second index) is in the inner loop
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
A[i][j] = A[i][j] + 1.0;
}
}
Interoperability Rule of Thumb¶
When passing a matrix from C/Python to Fortran (or interfacing via C-bindings/f2py), a matrix with shape (M, N) in C/Python appears as shape (N, M) in Fortran , or as the transposed matrix if shapes are kept identical.
Example of array storarge in column major order¶
program main
implicit none
real(8) :: frac(4, 3), fracc(3,4)
integer :: i
! Column 1 (x's): 0.0, 0.0, 0.1, 0.1
! Column 2 (y's): 0.2, 0.22, 0.2, 0.29
! Column 3 (z's): 0.3, 0.31, 0.3, 0.322
frac = reshape([0.0d0, 0.0d0, 0.1d0, 0.1d0, &
0.2d0, 0.22d0, 0.2d0, 0.29d0, &
0.3d0, 0.31d0, 0.3d0, 0.322d0], [4, 3])
fracc = reshape(frac, [3,4])
! Note: Correct syntax is `print *,` (no comma between `print` and `*`)
print *, "--- Unformatted Output (print *, frac) ----- "
print* " <column stored values of frac matrix are output>"
print *, frac
print *
print *, "--- Formatted Row-by-Row Output ---"
do i = 1, 4
write(*, '(3F10.4)') frac(i, :)
end do
print *, "--- Formatted Row-by-Row Output for fracc---"
do i = 1, 3
write(*, '(4F10.4)') fracc(i, :)
end do
end program main
output: note the garbage numbers after 16 decimal place
--- Unformatted Output (print *, frac) -----
`<column stored values of frac matrix are output>`
0.0000000000000000 0.0000000000000000 0.10000000000000001 0.10000000000000001 0.20000000000000001 0.22000000000000000 0.20000000000000001 0.28999999999999998 0.29999999999999999 0.31000000000000000 0.29999999999999999 0.32200000000000001
--- Formatted Row-by-Row Output ---
0.0000 0.2000 0.3000
0.0000 0.2200 0.3100
0.1000 0.2000 0.3000
0.1000 0.2900 0.3220
--- Formatted Row-by-Row Output for fracc---
0.0000 0.1000 0.2000 0.3100
0.0000 0.2000 0.2900 0.3000
0.1000 0.2200 0.3000 0.3220
Matrix-matrix multiplication in Fortran and Python¶
1. Matrix Multiplication Mechanics¶
The matrix multiplication @ computes the matrix product of frac and lattice.
- Shapes:
frachas shape \((4, 3)\) andlatticehas shape \((3, 3)\). The resultcarthas shape \((4, 3)\). - Mathematical Formula: Each element \((i, j)\) in
cartis computed as the dot product of row \(i\) offracand column \(j\) oflattice:
Row-Vector Interpretation¶
In physical terms, each row \(i\) of frac contains fractional coordinates \([x_i, y_i, z_i]\), and the rows of lattice are the basis vectors \(\mathbf{a}\), \(\mathbf{b}\), and \(\mathbf{c}\):
The \(i\)-th row of cart is constructed as a linear combination of the lattice vectors:
Step-by-Step Numerical Example (Row 0)¶
For the first atom (frac[0] = [0, 0.2, 0.3]):
2. Difference Between Row-Based and Column-Based Multiplications¶
In textbook linear algebra and physics, coordinates are usually written as column vectors :
where \(\mathbf{L}_{\text{col}}\) has lattice vectors stored as columns :
Comparison Table¶
| Property | Row-Vector Representation (frac @ lattice) | Column-Vector Representation (Lcol⋅r) |
|---|---|---|
| Atom Shape | \((N, 3)\)—**\(N\)**rows, 3 spatial coordinates | \((3, N)\)— 3 spatial coordinates,**\(N\)**columns |
| Lattice Vectors | Stored as**rows**in lattice |
Stored as**columns**in**\(L_{\text{col}}\)** |
| Formula | \(\text{cart} = \text{frac} \cdot L_{\text{row}}\) | \(\text{cart}^T = L_{\text{col}} \cdot \text{frac}^T\) |
| Transpose Relation | \(L_{\text{row}} = L_{\text{col}}^T\) | \(L_{\text{col}} = L_{\text{row}}^T\) |
Why NumPy uses Row Vectors¶
NumPy arrays use C-contiguous memory layout (row-major). Storing \(N\) atomic coordinates as shape \((N, 3)\) keeps each atom's \([x, y, z]\) contiguous in memory, making row-based operations (cart = frac @ lattice) faster and more intuitive in data scientific pipelines (such as ASE or PyMatGen).
3. How This is Done in Fortran¶
Fortran uses column-major memory layout and 1-based indexing .
Method A: Using Fortran's Intrinsic MATMUL¶
If frac is shape (4, 3) and lattice is shape (3, 3), the syntax is identical to Python using MATMUL:
Fortran
program main
implicit none
real(8) :: frac(4, 3), lattice(3, 3), cart(4, 3)
frac = reshape([0.0d0, 0.0d0, 0.1d0, 0.1d0, &
0.2d0, 0.22d0, 0.2d0, 0.29d0, &
0.3d0, 0.31d0, 0.3d0, 0.322d0], [4, 3])
lattice = reshape([2.0d0, 10.0d0, 10.1d0, &
0.0d0, 1.22d0, 0.2d0, &
1.3d0, 0.31d0, 0.3d0], [3, 3])
! Matrix multiplication
cart = matmul(frac, lattice)
end program main
Method B: Explicit Nested Loops¶
Fortran
integer :: i, j, k
do j = 1, 3
do i = 1, 4
cart(i, j) = 0.0d0
do k = 1, 3
cart(i, j) = cart(i, j) + frac(i, k) * lattice(k, j)
end do
end do
end do
Method C: High-Performance BLAS (DGEMM)¶
For large production codes, calling BLAS DGEMM is standard:
Fortran