Computational Physics - Department of Physics

6.3 Programming Details 165

for(k=0 ; k < n ; k++){

a[i][j]+=b[i][k]*c[k][j]

}

}

}

and in Fortran we have

DOj=1, n

DOi=1, n

DOk = 1, n

a(i,j)=a(i,j)+b(i,k)*c(k,j)

ENDDO

ENDDO

ENDDO

However, Fortran has an intrisic function called MATMUL, and the above three loops can
be coded in a single statementa=MATMUL(b,c). Fortran contains several array manipulation
statements, such as dot product of vectors, the transpose ofa matrix etc etc. The outer
product of two vectors is however not included in Fortran. The coding of Eq. (6.6) takes then
the following form in C++

for(i=0 ; i < n ; i++){

for(j=0 ; j < n ; j++){

a[i][j]+=x[i]*y[j]

}

}

and in Fortran we have

DOj=1, n

DOi=1, n

a(i,j)=a(i,j)+x(j)*y(i)

ENDDO

ENDDO

A matrix-matrix multiplication of a generaln×nmatrix with

a(i,j) =a(i,j)+b(i,k)∗c(k,j),

in its inner loops requires a multiplication and an addition. We define now a flop (floating
point operation) as one of the following floating point arithmetic operations, viz addition,
subtraction, multiplication and division. The above two floating point operations (flops) are
donen^3 times meaning that a general matrix multiplication requires 2 n^3 flops if we have
a square matrix. If we assume that our computer performs 109 flops per second, then to
perform a matrix multiplication of a 1000 × 1000 case should take two seconds. This can be
reduced if we multiply two matrices which are upper triangular such as

A=







a 11 a 12 a 13 a 14

0 a 22 a 23 a 24

0 0 a 33 a 34

0 0 0 a 44





.

The multiplication of two upper triangular matricesBCyields another upper triangular matrix
A, resulting in the following C++ code

for(i=0 ; i < n ; i++){

for(j=i ; j < n ; j++){

for(k=i ; k < j ; k++){