Computational Physics - Department of Physics

(Axel Boer) #1

6.3 Programming Details 167


P 1 =A 11 B 11 ,U 1 =P 1 +P 2 ,
P 2 =A 12 B 21 ,U 2 =P 1 +P 4 ,
P 3 =S 1 T 1 , U 3 =U 2 +P 5 ,
P 4 =S 2 T 2 , U 4 =U 3 +P 7 ,
P 5 =S 3 T 3 , U 5 =U 3 +P 3 ,
P 6 =S 4 B 22 ,U 6 =U 2 +P 3 ,
P 7 =A 22 T 4 ,U 7 =U 6 +P 6 ,

(6.10)

to find the resultingCmatrix as [
C 11 C 12
C 21 C 22


]

=

[

U 1 U 7

U 4 U 5

]

. (6.11)

In spite of the seemingly additional work, we have reduced the number of multiplications
from eight to seven. Since the multiplications are the computational bottleneck compared to
addition and subtraction, the number of flops are reduced.
In the case of squaren×nmatrices withnequal to a power of two,n= 2 m, the divided
blocks will haven 2 = 2 m−^1. Lettingf(m)be the number of flops needed for the full matrix and
applying Strassen recursively we find the total number of flops to be


f(m) = 7 f(m− 1 ) = 72 f(m− 2 ) =···= 7 mf( 0 ), (6.12)

wheref( 0 )is the one floating-point operation needed for multiplication of two numbers (two
20 × 20 matrices). For large matrices this can prove efficient, yielding a much better scaling,


O( 7 m) =O

(

2 log^27
m)
=O

(

2 mlog^27

)

=O

(

nlog^27

)

≈O

(

n^2.^807

)

, (6.13)

effectively saving 7 / 8 = 12 .5%each time it is applied.


6.3.3.2 Fortran Allocate Statement and Mathematical Operations on Arrays


An array is declared in the declaration section of a program,module, or procedure using the
dimension attribute. Examples include


REAL,DIMENSION(10) :: x,y
REAL,DIMENSION(1:10) :: x,y
INTEGER,DIMENSION(-10:10) :: prob
INTEGER,DIMENSION(10,10) :: spin

The default value of the lower bound of an array is 1. For this reason the first two state-
ments are equivalent to the first. The lower bound of an array can be negative. The last two
statements are examples of two-dimensional arrays.
Rather than assigning each array element explicitly, we canuse an array constructor to
give an array a set of values. An array constructor is a one-dimensional list of values, sepa-
rated by commas, and delimited by "(/" and "/)". An example is


a(1:3) = (/ 2.0, -3.0, -4.0 /)

is equivalent to the separate assignments


a(1) = 2.0
a(2) = -3.0
a(3) = -4.0
Free download pdf