6.3 Programming Details 163
double∗∗A =⇒double∗A[ 0 ... 3 ]
A[ 0 ][ 0 ] A[ 0 ][ 1 ] A[ 0 ][ 2 ] A[ 0 ][ 3 ]
A[ 1 ][ 0 ] A[ 1 ][ 1 ] A[ 1 ][ 2 ] A[ 1 ][ 3 ]
A[ 2 ][ 0 ] A[ 2 ][ 1 ] A[ 2 ][ 2 ] A[ 2 ][ 3 ]
A[ 3 ][ 0 ] A[ 3 ][ 1 ] A[ 3 ][ 2 ] A[ 3 ][ 3 ]
A[ 0 ]
A[ 1 ]
A[ 2 ]
A[ 3 ]
Fig. 6.2Conceptual representation of the allocation of a matrix in C++.
First of all, when we think of ann×nmatrix in Fortran and C++, we typically would
have a mental picture of a two-dimensional block of stored numbers. The computer stores
them however as sequential strings of numbers. The latter could be stored as row-major
order or column-major order. What do we mean by that? Recalling that for our matrix el-
ementsai j,irefers to rows and jto columns, we could store a matrix in the sequence
a 11 a 12 ...a 1 na 21 a 22 ...a 2 n...annif it is row-major order (we go along a given rowiand pick up all
column elementsj) or it could be stored in column-major ordera 11 a 21 ...an 1 a 12 a 22 ...an 2 ...ann.
Fortran stores matrices in the latter way, i.e., by column-major, while C++ stores them
by row-major. It is crucial to keep this in mind when we are dealing with matrices, because
if we were to organize the matrix elements in the wrong way, important properties like the
transpose of a real matrix or the inverse can be wrong, and obviously yield wrong physics.
Fortran subscripts begin typically with 1 , although it is no problem in starting with zero, while
C++ starts with 0 for the first element. This means thatA( 1 , 1 )in Fortran is equivalent to
A[ 0 ][ 0 ]in C++. Moreover, since the sequential storage in memory means that nearby matrix
elements are close to each other in the memory locations (andthereby easier to fetch) ,
operations involving e.g., additions of matrices may take more time if we do not respect the
given ordering.
To see this, consider the following coding of matrix addition in C++ and Fortran. We have
n×nmatricesA,BandCand we wish to evaluateA=B+Caccording to Eq. (6.2). In C++
this would be coded like
for(i=0 ; i < n ; i++){
for(j=0 ; j < n ; j++){
a[i][j]=b[i][j]+c[i][j]
}
}