White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI
Подождите немного. Документ загружается.
240 CHAPTER 6. HIGH PERFORMANCE COMPUTING
Figure 6.1.4: Bit Adder
6.1.2 Applied Area
Vector pipelines were introduced so as to make greater use of the register level
hardware. We will focus on the operation of floating point addition, which
requires four distinct steps for each addition. The segments of the device that
execute these steps are only busy for one fourth of the time to perform a floating
point add. The objective is to d esign computer hardware so that all of the
segments will be busy most of the time. In the case of the four segment floating
point adder this could give a speedup possibly close to four.
A vector pipeline is a register level device, which is usually in either the
control unit or the arithmetic logic unit. It has a collection of distinct hardware
modules or segments that execute the steps of an operation and each segment is
vector floating point adder in the arithmetic logic unit. The first pair of floating
point numbers is denoted by D1, and this pair enters the p ipeline in the upper
left in the figure. Segment CE on D1 is done during the first clock cycle. During
the second clock cycle D1 moves to segment AE, and the second pair of floating
point numbers D2 enters segment CE. Continue this process so that after three
clock cycles the pip eline is full and a floating point add is produced every clo ck
cycle. So, for large number of floating point adds with four segments the ideal
speedup is four.
6.1.3 Model
A discrete model for the speedup of a particular pipeline is as follows. Such
mo dels are often used in the design phase of computers. Also they are used to
determine how to best utilize vector pipelines on a selection of applications.
Vector Pipeline Timing Model.
Let N = time for one clock cycle,
Q = number of items of data,
O = number of segments,
V
y
= startup time for the vector pipeline operation and
V
v
= startup time for a serial operation.
Tserial = serial time = V
v
+ (ON)Q .
© 2004 by Chapman & Hall/CRC
required to be busy once the deviceisfull. Figure 6.1.5 depicts afour segment
6.1. VECTOR COMPUTERS AND MATRIX PRODUCTS 241
Figure 6.1.5: Vector Pipeline for Floating Point Add
Tvector = vector time = [V
y
+ (O 1)N] + NQ.
Speedup = Tserial/Tvector.
The vector startup time is usually much larger than the serial startup time.
So, for small amounts of data (small N), the serial time may be smaller than
the vector time! The vector pipeline does not start to become eective u ntil
it is full, and this takes
V
y
+ (O 1)N clock cycles. Note that the speedup
approaches the number of segments
O as Q increases.
Another important consideration in the use of vector pipelines is the orderly
and timely input of data. An example is matrix-vector multiplication and the
Fortran programming language. Fortran stores arrays by listing the columns
from left to right. So, if one wants to input data from the rows of a matrix,
then this will be in stride equal to the length of the columns. On the other
hand, if one inputs data from columns, then this will be in stride equal to one.
For this reason, when vector pipelines are used to do matrix-vector products,
the ji version performs much better than the i j version. This can have a very
significant impact on the e
ective use of vector computers.
6.1.4 Implementation
In order to illustrate the benefits of vector pipelines, consider the basic matrix-
vector product. The ij method uses products of rows times the column vector,
and the ji method uses linear combinations of the columns. The advantage of
the two methods is that it often allows one to make use of particular properties
of a computer such as communication speeds and local versus main memory
size. We shall use the following notation:
{ = [{
m
] is a column vector where m = 1> ===> q and
D = [d
l>m
] is a matrix where l = 1> ===> p are the row numbers
and
m = 1> ===> q are the column numbers.
© 2004 by Chapman & Hall/CRC
242 CHAPTER 6. HIGH PERFORMANCE COMPUTING
Matrix-vector Product (ij version) Ax = d.
for
l = 1> p
g
l
= 0
for m = 1> q
g
l
= g
l
+ d
l>m
{
m
endloop
endloop.
An alternate way to do matrix-vector products is based on the following
reordering of the arithmetic operations. Consider the case where
q = p = 3
and
d
1>1
{
1
+ d
1>2
{
2
+ d
1>3
{
3
= g
1
d
2>1
{
1
+ d
2>2
{
2
+ d
2>3
{
3
= g
2
d
3>1
{
1
+ d
3>2
{
2
+ d
3>3
{
3
= g
3
=
This can be written in vector form as
5
7
d
1>1
d
2>1
d
3>1
6
8
{
1
+
5
7
d
1>2
d
2>2
d
3>2
6
8
{
2
+
5
7
d
1>3
d
2>3
d
3>3
6
8
{
3
=
5
7
g
1
g
2
g
3
6
8
=
In other words the product is a linear combination of the columns of the matrix.
This amounts to reversing the order of the nested loop in the above algorithm.
The matrix-matrix products are similar, but they will have three nested loops.
Here there are 6 di
erent orderings of the loops, and so, the analysis is a little
more complicated.
Matrix-vector Product (ji version) Ax = d.
g = 0
for
m = 1> q
for l = 1> p
g
l
= g
l
+ d
l>m
{
m
endloop
endloop.
Carolina Supercomputing Center. All calculations were for a matrix-vector
product where the matrix was 500 × 500 for 2(500)
2
= 5 · 10
5
floating point
operations. In the ji Fortran code the inner loop was vectorized or not vectorized
by using the Cray directives !dir$ vector and !dir$ novector just before the loop
30 in the following code segment:
© 2004 by Chapman & Hall/CRC
The calculations in Table 6.1.2 were done on the Cray T916 at the North
6.1. VECTOR COMPUTERS AND MATRIX PRODUCTS 243
Table 6.1.2: Matrix-vector Computation Times
Method Time x 10
4
sec.
ji (vec. on, -O1) 008.86
ji (vec. o, -O1) 253.79
ij (vec. on, -O1) 029.89
ij (vec. o, -O1) 183.73
matmul (-O2) 033.45
ji (vec. on, -O2) 003.33
do 20 j = 1,n
!dir$ vector
do 30 i = 1,n
prod(i) = prod(i) + a(i,j)*x(j)
30 continue
20 continue.
The first two calculations indicate a speedup of over 28 for the ji method. The
next two calculations illustrate that the ij method is slower than the ji method.
This is because Fortran stores numbers of a two dimensional array by columns.
Since the ij method gets rows of the array, the input into the vector pipe will
be in stride equal to n. The fifth calculation used the f90 intrinsic matmul
for matrix-vector products. The last computation used full optimization, -O2.
The loops 20 and 30 were recognized to be a matrix-vector product and an
optimized BLAS2 subroutine was u sed. BLAS2 is a collection of basic lin-
ear al gebra subroutines with order
q
2
an overall speedup equal to about 76. The floating point operations per second
for this last computation was (5 10
5
)/ (3.33 10
4
) or about 1,500 megaflops.
6.1.5 Assessment
Vector pip es can be used to do a number of operations or combination of oper-
ations. The potential speedup dep ends on the number of segments in the pip e.
If the computations are more complicated, then the speedups will decrease or
the code may not be vectorized. There must be an orderly input of data into
the pipe. Often there is a special memory called a register, which is used to
input data into the pip e. Here it is important to make optimal use of the size
and speed of these registers. This can be a time consuming e
ort, but many of
the BLAS subroutines have been optimized for certain computers.
Not all computations have independent parts. For example, consider the
following calculations
d(2) = f + d(1) and d(3) = f + d(2). The order of the
two computations will give di
erent results! This is called a data dependency.
Compilers will try to detect these, but they may or may not find them. Basic
iterative methods such as Euler or Newton methods have data dependencies.
© 2004 by Chapman & Hall/CRC
operations, see http://www.netlib.org or
http://www.netlib.org/blas/sgemv.f. This gave an additional speedupover2for
244 CHAPTER 6. HIGH PERFORMANCE COMPUTING
Some calculations can be reordered so that they have independent parts. For
example, consider the traditional sum of four numbers ((
d(1) + d(2)) + d(3))) +
d(4)= This can be reordered into partial sums (d(1) + d(2)) + (d(3) + d(4)) so
that two processors can be used.
6.1.6 Exercises
1. Write a Fortran or C or MATLA B code for the ij and ji matrix vector
products.
2. Write a Fortran or C or M
ATLA B code so that the BLAS subroutine
sgemv() is used.
3. On your computing environment do calculations for a matrix-vector prod-
the computing times.
4. Consider the explicit finite di
erence methods in Sections 1.2-1.5.
(a). Are the calculations in the outer time loops independent, vectoriz-
able and why?
(b). Are the calculations in the inner space loops independent, vectoriz-
able and why?
6.2 Vector Computations for Heat Diusion
6.2.1 Introduction
Consider heat conduction in a thin plate, which is thermally insulated on its
surface. The mo del of the temperature will have the form x
n+1
= Dx
n
+ e
for the time dependent case and x = Dx + e for the steady state case. In
general, the matrix D can be extremely large, but it will also have a special
structure with many more zeros than nonzero components. Here we will use
vector pipelines to execute this computation, and we will also extend the model
to heat di
usion in three directions.
6.2.2 Applied Area
Previously we considered the model of heat diusion in a long thin wire and in
a thin cooling fin. The temperature was a function of one or two space variables
and time. A more realistic model of temperature requires it to be a function of
three space variables and time. Consider a cooling fin that has di
usion in all
three space directions as discussed in Section 4.4. The initial temperature of
the fin will be given and one hot surface will be specified as well as the other
five cooler surfaces. The objective is to predict the temperature in the interior
of the fin in order to determine the e
ectiveness of the cooling fin.
© 2004 by Chapman & Hall/CRC
uct similar to those in Table 6.1.2 that were done on the CrayT916. Compare
6.2. VECTOR COMPUTATIONS FOR HEAT DIFFUSION 245
6.2.3 Model
The model can be formulated as either a continuous model or as a discrete
model. For appropriate choices of time and space steps the solutions should
be close. In order to generate a 3D time dependent model for heat transfer
diusion, the Fourier heat law must be applied to the x, y and z directions.
The continuous and discrete 3D models are very similar to the 2D versions.
In the continuous 3D model the temperature
x will depend on four variables,
x({> |> }> w). In (6.2.1) (Nx
}
)
}
models the diusion in the z direction where
the heat is entering and leaving the top and bottom of the volume
{|}=
Continuous 3D Model for x = x({> |> }> w)=
fx
w
(Nx
{
)
{
(Nx
|
)
|
(Nx
}
)
}
= i (6.2.1)
x({> |> }> 0) = given and (6.2.2)
x({> |> }> w) = given on the boundary. (6.2.3)
Explicit Finite Di
erence 3D Model for x
n
l>m>o
x(l{> m|> o}> nw)=
x
n+1
l>m>o
= (w@f)i
n
l>m>o
+ (1 )x
n
l>m>o
+w@{
2
(x
n
l+1>m>o
+ x
n
l1>m>o
)
+w@|
2
(x
n
l>m
+1>o
+ x
n
l>m
1>o
)
+
w@}
2
(x
n
l>m>o+1
+ x
n
l>m>o1
) (6.2.4)
= (N@f)w(2@{
2
+ 2@|
2
+ 2@}
2
)>
l> m> o = 1> ==> q 1 and n = 0> ==> pd{n 1,
x
0
l>m>o
= given> l> m> o = 1> ==> q 1 and (6.2.5)
x
n
l>m>o
= given> n = 1> ===> pd{n, l> m> o on the boundary grid. (6.2.6)
Stability Condition.
1
((N@f)w(2@{
2
+ 2@|
2
+ 2@}
2
)) A 0.
6.2.4 Method
The computation of the above explicit model does not require the solution of a
linear algebraic system at each time step, but it does require the time step to be
suitably small so that the stability condition holds. Since there are three space
directions, there are three indices associated with these directions. Therefore,
there will be four nested loops with the time loop being the outer loop. The
inner three loops can be in any order because the calculations for
x
n+1
l>m>o
are
independent with respect to the space indices. This means one could use vector
or multiprocessing computers to do the inner loops.
© 2004 by Chapman & Hall/CRC
246 CHAPTER 6. HIGH PERFORMANCE COMPUTING
Table 6.2.1: Heat Diusion Vector Times
Loop Length Serial Time Vector Time Speedup
30 (Alliant) 04.07 2.60 1.57
62 (Alliant) 06.21 2.88 2.16
126 (Alliant) 11.00 4.04 2.72
254 (Alliant) 20.70 6.21 3.33
256(Cray) 1.36 .0661 20.57
512(Cray) 2.73 .1184 23.06
6.2.5 Implementation
The following code segment for 2D diusion was run on the Cray T916 com-
puter. In code the Cray directive !dir$ vector before the beginning of loop 30
instructs the compiler to use the vector pipeline on loop 30. Note the index for
loop 30 is i, and this is a row nu mber in the array u. This ensures that the
vector pipe can sequentially access the data in u.
do 20 k = 1,maxit
!dir$ vector
do 30 i = 2,n
u(i,k+1) = dt*f + alpha*(u(i-1,k) + u(i+1,k))
+ (1 - 2*alpha)*u(i,k)
30 continue
20 continue.
Table 6.2.1 contains vector calculations for an older computer called the
Alliant FX/40 and for the Cray T916. It indicates increased speedup as the
length of loop 30 increases. This is because the startup time relative to the
execution time is decreasing. The Alliant FX/40 has a vector pipe with four
segments and has a limit of four for the speedup. The Cray T916 has more
segments in its vector pipeline with speedups of about 20. Here the speedup
for the Cray T916 is about 20 because the computations inside lo op 30 are a
little more involved than those in the matrix-vector product example.
The M
AT LAB code heat3d.m is for heat diusion in a 3D cooling fin, which
has initial temperature equal to 70, and with temperature at the boundary { = 0
equal to 370 for the first 50 time steps and then set equal to 70 after 50 time
steps. The other temperatures on the boundary are always equal to 70. The
code in heat3d.m generates a 4D array whose entries are the temperatures for
3D space and time. The input data i s given in lines 1-28, the finite di
erence
method is executed in the four nested loops in lines 37-47, and some of the
output is graphed in the 3D plot, using the M
ATLA B command slice in line
50. The M
AT LAB commands slice and pause allow the user to view the heat
moving from the hot mass towards the cooler sides of the fin. This is much
co e
!cient, 1-alpha, in line 41 must be positive so that the stability condition
© 2004 by Chapman & Hall/CRC
more interesting than the single grayscale plot in Figure 6.2.1 at time 60. The
6.2. VECTOR COMPUTATIONS FOR HEAT DIFFUSION 247
holds.
MATLAB Code heat3d.m
1. % Heat 3D Diusion.
2. % Uses the explicit method.
3. % Given boundary conditions on all sides.
4. clear;
5. L = 2.0;
6. W = 1.0;
7. T = 1.0;
8. Tend = 100.;
9. maxk = 200;
10. dt = Tend/maxk;
11. nx = 10.;
12. ny = 10;
13. nz = 10;
14. u(1:nx+1,1:ny+1,1:nz+1,1:maxk+1) = 70.; % Initial temperature.
15. dx = L/nx;
16. dy = W/ny;
17. dz = T/nz;
18. rdx2 = 1./(dx*dx);
19. rdy2 = 1./(dy*dy);
20. rdz2 = 1./(dz*dz);
21. cond = .001;
22. spheat = 1.0;
23. rho = 1.;
24. a = cond/(spheat*rho);
25. alpha = dt*a*2*(rdx2+rdy2+rdz2);
26. x = dx*(0:nx);
27. y = dy*(0:ny);
28. z = dz*(0:nz);
29. for k=1:maxk+1 % Hot side of fin.
30. time(k) = (k-1)*dt;
31. for l=1:nz+1
32. for i=1:nx+1
33. u(i,1,l,k) =300.*(time(k)
?50)+ 70.;
34. end
35. end
36. end
37. for k=1:maxk % Explicit method.
38. for l=2:nz
39. for j = 2:ny
40. for i = 2:nx
41. u(i,j,l,k+1) =(1-alpha)*u(i,j,l,k) ...
42. +dt*a*(rdx2*(u(i-1,j,l,k)+u(i+1,j,l,k))...
© 2004 by Chapman & Hall/CRC
248 CHAPTER 6. HIGH PERFORMANCE COMPUTING
Figure 6.2.1: Temperature in Fin at t = 60
43. +rdy2*(u(i,j-1,l,k)+u(i,j+1,l,k))...
44. +rdz2*(u(i,j,l-1,k)+u(i,j,l+1,k)));
45. end
46. end
47. end
48. v=u(:,:,:,k);
49. time(k)
50. slice(x,y,z,v,.75,[.4 .9],.1)
51. colorbar
52. pause
53. end
6.2.6 Assessment
The explicit time discretization is an example of a method that is ideally vector-
izable. The computations in the inner space loops are independent so that the
inner loops can be executed using a vector or multiprocessing computer. How-
ever, the stability condition on the step sizes can still be a s erious constraint. An
alternative is to use an implicit time discretization, but as indicated in Section
4.5 this generates a sequence of linear systems, which require some additional
computations at each time step.
© 2004 by Chapman & Hall/CRC
6.3. MULTIPROCESSORS AND MASS TRANSFER 249
6.2.7 Exercises
1. Duplicate the calculations in heat3d.m. Experiment with the slice para-
meters in the MATLAB command slice.
2. In heat3d.m experiment with dierent time mesh sizes, pd{n = 150> 300
and 450. Be sure to consider the stability constraint.
3. In heat3d.m exp eriment with di
erent space mesh sizes, q{ or q| or
q} = 10> 20 and 40. Be sure to consider the stability constraint.
4. In heat3d.m experiment with dierent thermal conductivities N = frqg =
=01> =02 and .04. Be sure to make any adjustments to the time step so that the
stability condition holds.
5. Suppose heat is being generated at a rate of 3 units of heat per unit
volume per unit time.
(a). Modify heat3d.m to implement this source of heat.
(b). Experiment with di
erent values for this heat source i = 0> 1> 2 and
3
=
6.3 Multiprocessors and Mass Transfer
6.3.1 Intro duction
Since computations for 3D heat diusion require four nested loops, the compu-
tational demands increase. In such cases the use of vector or multiprocessing
computers could be very e
ective. Another similar application is the concentra-
tion of a pollutant as it is dispersed within a deep lake. Here the concentration
is a function of time and three space variables. This problem, like heat di
usion
in 3D, will also require more computing p ower. In this section we will describe
and use a multiprocessing computer.
A multiprocessing computer is a computer with more than one "tightly"
coupled CPU. Here "tightly" means that there is relatively fast communication
among the CPUs. There are several classification schemes that are commonly
used to describe various multiprocessors: memory, communication connections
and data streams.
Two examples of the memory classification are shared and distributed. The
shared memory multiprocessors communicate via the global shared memory,
sor. Shared memory multiprocessors often have in-code directives that indicate
the code segments to be executed concurrently. The distributed memory mul-
tiprocessors communicate by explicit message passing, which must be part of
memory computers. In these depictions each node could have several CPUs,
for example some for computation and one for communication. Another illus-
tration is each no de could be a shared memory computer, and the IBM/SP is
a particular example of this.
© 2004 by Chapman & Hall/CRC
and Figure 6.1.2 is adepiction of afour processor shared memory multiproces-
the computer code. Figures 6.3.1 and 6.3.2 illustrate three types of distributed