=Paper=
{{Paper
|id=Vol-1482/084
|storemode=property
|title=Parallel algorithm for solving large-scale dynamic general equilibrium models
|pdfUrl=https://ceur-ws.org/Vol-1482/084.pdf
|volume=Vol-1482
}}
==Parallel algorithm for solving large-scale dynamic general equilibrium models==
https://ceur-ws.org/Vol-1482/084.pdf
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
Parallel algorithm for solving large-scale dynamic
general equilibrium models
N. B. Melnikov1, A. P. Gruzdev1, M. G. Dalton2, B. C. O’Neill3
Lomonosov Moscow State University1
National Oceanic and Atmospheric Administration, USA2
National Center for Atmospheric Research, USA3
We present a parallel algorithm for computing an equilibrium path in a large-scale eco-
nomic growth model. We exploit the special block structure of the nonlinear systems
of equations common in such models. Our algorithm is based on an iterative method
of Gauss-Seidel type with prices of different time periods calculated simultaneously
rather than recursively. We have implemented the parallel algorithm in OpenMP
and MPI programming environments. The numerical results show that speedup im-
proves almost linearly as number of nodes increases. Different methods for solving an
individual block: Newton-type methods, Krylov subspace methods and trust-region
methods, give similar results for the speedup.
1. Introduction
Dynamic general equilibrium (GE) models are used to quantify the effects on economy due
to demographic, technological and climate change (see, e.g., [1, 2] and refs. therein). Dynamic
GE is described by large-scale systems of nonlinear equations. One of the the most popular early
methods for solving such nonlinear systems was the Fair-Taylor method [3], which is a variation
of the Gauss-Seidel method. Recent contributions have concentrated on Newton-type methods
that are now more widely used (see, e.g., [4]).
In this paper, we give the Fair-Taylor method a fresh look. We use the special block
structure of the equations system common in the dynamic GE models to describe and implement
a parallel algorithm. The idea of taking block structure into account was used previously for
multi-country models (see, e.g., [5]). We implement parallel calculation of individual time blocks
using OpenMP and MPI environments for systems with shared and distributed memory. For
the solution of the time blocks, we compare the performance of different Newton-type methods:
LU -factorization, Krylov methods and trust-region methods. The algorithm is demonstrated in
the Population–Environmental–Technology (PET) model (see [1, 2] and refs. therein).
The paper is organized as follows. In Section 2 we present the block Gauss-Seidel method
and briefly describe the Newton-type methods that we use for solving the equations system of
individual blocks. In section 3, by example of the PET model, we derive the nonlinear system
of equations that determines the GE. In section 4 we implement the sequential and two parallel
versions (OpenMP and MPI) of the algorithm. Finally, in section 5 we present and discuss the
numerical results.
2. Block Gauss-Seidel method
First, we briefly recall the Gauss-Seidel method for a linear equation system Ax = b, where
A is a nondegenerate n × n matrix with nonzero diagonal elements (see, e.g., [6]). Assume that
the kth iterate xk = (xk1 , . . . , xkn )T and i − 1 components xk+1 k+1
1 , . . . , xi−1 of the (k + 1)th iterate
xk+1 have been determined. Then to obtain xk+1 i the ith equation of the system
i−1
X n
X
aij xk+1
j + aii xi + aij xkj = bi ,
j=1 j=i+1
84
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
is solved with respect to xi and the solution is taken as xk+1
i . We write method in a more
compact form by splitting the matrix as A = D + L + U , where D is diagonal and L and U are
strictly lower and upper triangular parts. Then
xk+1 = M xk + c,
where M = −(D + L)−1 U is the iteration matrix and c = (D + L)−1 b. The usual termination
criterium is ||rk ||/||b|| < ε, where rk = b − Axk is the residual. Since
rk = (D + L)(xk+1 − xk ),
we can use the properly scaled error kxk+1 − xk k to terminate iterations.
The same idea is generalized to a nonlinear equations system f (x) = 0, where f = (f1 , . . . , fn ).
To obtain xk+1
i the ith equation
fi (xk+1 k+1 k k
1 , . . . , xi−1 , xi , xi+1 , . . . , xn ) = 0
is solved with respect to xi and the solution is taken as xk+1
i . An advantage of the method is
that one does not need to calculate the Jacobian of the system but only values of the functions
fi (x).
Now we assume that fi and xi are vectors:
fi = (fi1 , . . . , fij , . . . , fil ), xi = (xi1 , . . . , xij , . . . , xil ),
where the blocks fij and xij are also vectors having the same dimension. Let the jth block of
functions fij depend only on the jth block of variables xij . Then to obtain xk+1
i the equations
fij (xk+1 k+1 k k
1 , . . . , xi−1 , x, xi+1 . . . , xn ) = 0, j = 1, l, (1)
are solved independently with respect to x and the solution is taken as xk+1
ij .
To solve the equations system for a separate block (1), which we denote F (x) = 0 for
now, it is reasonable to use a method with superlinear convergence. Here, we describe three
methods that are further used for comparison: Newton’s method, Krylov subspace method and
trust-region method. In each of the three methods, at the current iterate xk , the linearized
system
F ′ (xk )sk = −F (xk ) (2)
is solved with respect to the step sk in order to obtain the next iterate xk+1 = xk + sk .
In Newton’s method, a LU -factorization of the Jacobian F ′ (xk ) is used (see, e.g., [6]). The
inputs of the algorithm are the initial iterate x0 and a termination tolerance tol.
input : x0 , tol
output : the solution xk
1 k ← 0;
2 while ||F (xk )|| > tol do
3 solve (2) with respect to sk using the LU-factorization;
4 xk+1 ← xk + sk ;
5 k ← k + 1;
6 end
Figure 1. Newton’s method.
In the Newton Iterative Solver (NITSOL) [7, 8], instead of the direct method for solving
the linearized system, an inexact Newton’s method with backtracking is used. The termination
criterium is
||F ′ (xk )sk + F (xk )|| ≤ ηk ||F (xk )||, (3)
85
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
where ηk ∈ [0, 1) is the tolerance. To obtain the step sk , the generalized minimal residual (GM-
RES) Krylov subspace method is used (see, e.g., [9]). Once an initial sk has been determined,
it is tested and, if necessary, reduced in length through backtracking until an acceptable step is
obtained. The parameter t ∈ (0, 1) is used to judge sufficient reduction θ ∈ (θmin , θmax ). The
Jacobian is calculated using the difference approximation, just as in Newton’s method.
input : x0 , tol;
ηmax ∈ [0, 1), t ∈ (0, 1), 0 < θmin < θmax < 1
output : the solution xk
1 k ← 0;
2 while ||F (xk )|| > tol do
3 solve (2) with respect to sk using a Krylov subspace method (termination
criterium (3));
4 while ||F (xk + sk )|| > [1 − t(1 − ηk )]||F (xk )|| do
5 ηk ← 1 − θ(1 − ηk );
6 sk ← θsk , θ ∈ [θmin , θmax ];
7 end
8 xk+1 ← xk + sk ;
9 k ← k + 1;
10 end
Figure 2. Inexact Newton with backtracking.
The third method solves the system of nonlinear equations using the trust-region method
to obtain the current iterate (see, e.g., [10]). In the routine NEQBF, the following problem
||F (xk ) + Bk sk ||2 −
→ min, (4a)
sk
||sk || ≤ δk . (4b)
is solved to get the direction sk , where δk is the size of the trust region and Bk is the approximate
Jacobian evaluated at the current point xk . Then, the function at the point xk+1 = xk + sk
is evaluated and used to decide whether the new point xk should be accepted. Problem (4) is
solved approximately by the double dogleg method. The Jacobian is approximated by Broyden’s
formula
([F (xk+1 ) − F (xk )] − Bk sk )sTk
Bk+1 = Bk + .
sTk sk
input : x0 , tol;
δ0 , θ ∈ (0, 1)
output : the solution xk
1 k ← 0;
2 while ||F (xk )|| < tol do
3 δk ← δ0 ;
4 while f (xk+1 ) ≥ f (xk ) do
5 solve (4) with respect to sk using the double dogleg method;
6 δk ← θδk ;
7 xk+1 ← xk + sk ;
8 end
9 k ← k + 1;
10 end
Figure 3. Trust-region method.
86
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
3. Dynamic general equilibrium model
In this section, we show that a dynamic GE model leads to a block equations system.
To be specific, we consider the one region PET model (see, e.g., [1, 2] for details). PET is a
neoclassical growth models with three types of agents: consumers, producers and government.
Consumers take the prices as given and solve the utility maximization problem over the infinite
time-horizon under their budget constraints (rational expectations). Producers minimize costs
given the levels of production and prices at each moment in time. The government plays a
neutral role redistributing the capital through taxes and transfers. Prices are determined by
the markets clearing conditions. The GE is defined as a solution to the nonlinear equations
system that consists of the first-order optimality conditions for consumers and producers and
supply-equals-demand conditions for markets.
3.1. Optimal economic growth
The consumer side of the model consists of Nc heterogeneous groups of households. Demand
of the ith consumer group is determined by the intertemporal utility function
∞
1X t
Ui (ci ) = β nit ui (cit ),
ψ
t=0
where t = 0, 1, ... is time, cit is the level of consumption (here and henceforth, small letters indi-
cate the per capita values), β ∈ (0, 1) is the discount rate, ψ ∈ (−∞, 1) is the time substitution
rate and nit is the size of population. The instantaneous utility function ui (cit ) is given by the
constant-elasticity of substitution (CES) function
ψ
Nc ρ
X
ρ
ui (cit ) = (µijt cijt ) ,
j=1
with constant electivity σ = 1/(1−ρ). Here, the index j = 1, Nc labels goods, ρ ∈ (−∞, 1)\ {0} is
the substitution parameter among goods and µijt is the preference coefficient. Capital dynamics
is described as
(1 + νit ) ki,t+1 = (1 − δ) kit + xit , ki0 > 0, (5)
where xit is investment, δ ∈ (0, 1) is the capital depreciation coefficient, 1 + νit = ni,t+1 /nit is
the population growth coefficient (νit is the growth rate). Budget constraint of the ith consumer
group at time t is
Nc
X
pjt cijt + qt xit = (1 − θit ) ωt lit + (1 − φit ) rt kit + git , (6)
j=1
where pjt is the price of jth consumer good, qt and rt are prices of investments and capital, ωt
is the wage rate, git is the government transfers, lit is the labor supply, θit and φit are the tax
rates on capital and labor incomes.
Variables kit , cijt and xit are determined by maximizing
Ui (ci ) −
→ max , (7)
cit ,kit ,xit
under constraints (5) and (6). It is convenient to split the solution of problem (5)–(7) into two
steps. The first one is to maximize the utility ui (cit ) at time t given the expenditures mit :
1
Nc ρ Nc
X X
ρ
(µij cij ) − → max, pj cij = mi ,
cij
j=1 j=1
87
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
with respect to consumer demand of different types of goods cij , i.e. a static problem (we omit
the index t where it does not lead to confusion). Solving the dual problem
1
Nc Nc ρ
X X
→ min,
pj cij − (µij cij )ρ = ci ,
cij
j=1 j=1
we obtain
ρ−1
Nc Nc ρ ρ
X X pj ρ−1
min pj cij = p̄i ci , p̄i = .
µij
j=1 j=1
The second step is to solve problem (5)–(7) for savings as a function of time (dynamics).
Rewriting (5)–(7) in terms of p̄it and cit :
∞
1X t
β nit cit ψ −
→ max ,
ψ cit ,kit
t=0
pit cit + qt xit = (1 − θit ) ωt lit + (1 − φit ) rt kit + git ,
(1 + νit ) ki,t+1 = (1 − δ) kit + xit ,
we obtain the first-order optimality condition (Euler equation)
qt ψ−1 qt+1 (1 − δ) + (1 − φi,t+1 )rt+1
cit =β ci,t+1 ψ−1 . (8)
pit pi,t+1
The transversality conditions
lim λit kit = 0, (9)
t−
→∞
where λit is the Lagrange multiplier, ensures the existence of the optimal trajectory (see,
e.g., [11]).
3.2. Duality theory for production functions
The production side of the PET model consists of competitive firms divided into sectors.
Each sector produces either a final good: NC consumer goods and “investment good”, or an
intermediate good: NE energy types and the rest, which is called materials (NX = NC +1+NE +1
is the total number of production sectors).
Production of the good X is determined by the the inputs of capital K, labor L, energy
composite Ē and materials M according to a CES function
1
X = γX αK (GK K)ρX + αL (GL L)ρX + αĒ (GĒ Ē)ρX + αM (GM M )ρX ρX
, (10)
where GI is the productivity factor for I = K, L, Ē, M and γX normalizes αI to sum to unity.
The parameters αI and GI can be sector and time dependent. The producer of the good X
minimizes the costs
PK K + PL L + PĒ Ē + (1 + τM )PM M −
→ min
K,L,Ē,M
under the given level of production (10) (τM is the ad-valorem tax on the use of materials). The
solution has the form
min PK K + PL L + PĒ Ē + (1 + τM )PM M = PX X,
88
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
where the price PX is 1
1
ρX 1
ρX
1 1−ρX PK ρX −1
1−ρX PL ρX −1
PX = αK + αL +
γX GK GL
1
ρX 1
ρX !
1−ρX PĒ ρX −1
1−ρX (1 + τM )PM ρX −1
+ αĒ + αM .
GĒ GM
The input-output coefficients AI = I/X for I = K, L, Ē are given by
1
1 PI ρX −1
AI = ρ
,
αI (γX GI ) PX
X
and for I = M by
1
1 (1 + τM )PM ρX −1
AM = ρ
.
αM (γX GM ) X PX
Similarly, production of the energy composite Ē is determined according to
! 1
X ρE
Ē = γE αEi (GEi Ei )ρE , (11)
i
where Ei , i = 1, NE are different energy types. Solving the cost-minimization problem
X
(1 + τEi )PEi Ei −
→ min
Ei
i
under the given the level of production (11), we derive the prices
ρ
ρE ! EρE
−1
1
1 X 1−ρE (1 + τEi )PEi ρE −1
PĒ = α Ei
γE GE i
i
and input-output shares AEi = Ei /Ē for the energy use:
1
1 (1 + τEi )PEi ρE −1
AEi = ρ
,
αEi (γE GEi ) E PĒ
where τEi is the ad-valorem tax on the energy use.
Government revenues include capital and labor income taxes on households and specific
taxes on output plus ad-valorem taxes on energy use and materials for each industry,
Nd NX NE
!
X X X
GREVt = ni (φi rki + θi ωli ) + τX j Xj + τESj PESj ESj + τMj Mj =
i=1 j=1 s=1
Nd NX NE
!
X X X
= (φi PK Ki + θi PL Li ) + τX j + τEsj PEsj AEsj AĒj + τMj AMj Xj .
i=1 j=1 s=1
Government expenditures are the sum of purchases and transfers, GEXPt = GPt + GTt . Gov-
ernment purchases are assumed to be constant in the current prices: GPt /p̄t = (1 + ξ)t GP0 /p̄0 ,
1
By the envelope theorem, the price is identified with the marginal cost of production (see, e.g., [12]).
89
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
where p̄t = (1/Nd ) N
P d
i p̄it . It is assumed that GP0 is given by a Cobb-Douglas production
function, so that the input-output coefficients are
αK αL αM αE i
AGP
K = , AGP
L = , AGP
M = , AGP
Ei = , i = 1, NE .
PK PL PM PEi
Similarly, government baseline transfers are neutral. There is additional lump-sum adjustment
LSAt to balance the government budget:
d N
p̄t X
GTt = (1 + ξ)t nit gi0 + LSAt .
p̄0
i=1
3.3. Equilibrium conditions
Aggregate supply for capital and labor at each point in time are determined by summing
over levels supplied by each consumer group:
Nd
X Nd
X
K AS = ni ki , LAS = ni li .
i=1 i=1
Aggregate demand for capital and labor are
NX NX
AjK Xj + AGP AjL Xj + AGP
X X
K AD = K GP, LAD = L GP.
j=1 j=1
Similarly, for consumer goods and investment.
T
Total demand for intermediate goods z = X1AD , . . . , XN
AD
E +1 is the sum of demands
to produce other intermediate goods, and those to produce final goods. Total demand for
intermediate inputs by final goods producers and government is equal to
T
NX NX
Aj1 XjAD + AGP AjNE +1 XjAD + AGP
X X
y= 1 GP, . . . , NE +1 GP
= (Y1 , . . . , YNE +1 )T .
NE +2 NE +2
Since the production functions are homogeneous functions, equilibrium implies
E +1
NX NX
Aji XjAD + Aji XjAD + AGP
X
XiAD = i GP, i = 1, NE + 1.
j=1 j=NE +2
or in the compact form z = Az + y, where A = (Aji ) is the (NE + 1) × (NE + 1) matrix of
input-output coefficients. Hence outputs of the intermediate goods sectors z are determined by
the final goods production y as
z = (I − A)−1 y,
where I is the (NE + 1) × (NE + 1) unity matrix.
A competitive equilibrium in each time moment is defined by markets clearing for factors
of production, capital and labor, final goods and a balanced government budget:
K AD = K AS , LAD = LAS , X AD = X AS , GREV = GEXP.
Note that the consumer prices pj and q are the gross (after-output-tax) prices PX of consumer
goods and investment good in industry. The rental rate of capital r and wage rate w are equal
to the prices of capital PK and labor PL in industry.
90
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
3.4. Computation of equilibrium
The PET model, as well as other dynamic GE models used in application, evaluates changes
during a transition period, and does not address possible effects on the long-run equilibrium.
Therefore, the infinite time horizon is replaced by a sufficiently large finite time interval [0, T ].
As t ≥ T the system is assumed to follow the balanced growth path, i.e. the quantities xit , kit ,
lit , etc. grow exponentially with the growth rate ξ. Proceeding to the limit t → − ∞ in (8), we
obtain
(1 − φ)r
(1 + ξ)1−ψ = β 1 − δ + . (12)
q
The boundary condition (12) at t = T is used instead of transversality condition at infinity (9).
This way, we get a two-point (discrete) boundary-value problem for each consumer group.
The nonlinear equations system describing the GE can be written as
f (K, P ) = 0,
where K is the capital and P are prices (at all time moment t = 0, T ). For brevity, we do
not write consumption, investment and transfers since they can be obtained from K and P .
The Gauss-Seidel method is applied to this system as follows. Assume that the sth iteration of
capital K s has been determined. Then solving the system
f (K s , P ) = 0
with respect to P , we obtain the (s + 1)th iteration of prices, P s+1 = P ; the previous iteration
P s is used as the initial point of the method. Different time blocks of this part of the algorithm,
called the inner loops, can be calculated in parallel.
Once P s+1 has been determined, the next iteration of capital is obtained by solving the
system
f (K, P s+1 ) = 0
with respect to K and setting K s+1 = K. The part of the algorithm that calculates K s+1 from
K s is called the outer loop (see [1, 3] for details of the capital iteration).
4. Algorithm implementation
First, we describe the sequential version (Algorithm 4). The input is the initial guess of the
capital K 0 as a function of time and the outputs are the capital K and prices P as functions of
time. Here, tol is the tolerance and maxit is the maximum number of iterations in the outer loop;
T is the length of the time interval. Two types of arrays are used to store capital, consumption,
prices, etc. as functions of time. The arrays of the first type, called stor, contain the results of
the it-th iterate of the outer loop over all time moments t. The arrays of the second type, called
dyn, contain the results of the it-th iterate of the outer loop for at the current t. Variable dif f
is the error at the current iterate of the outer loop.
In the OpenMP version, steps of the loop over time are performed in the parallel region
(Algorithm 5). All dyn and stor arrays are shared by all threads. The parameters arrays
are transferred through the copyin clause. The time steps are distributed into blocks using
schedule(auto).
In the MPI version, we use only one communicator with number of nodes nsize. Each node
has its own copy of the dyn and stor arrays. Moreover, each node works with its own part
of the array, and the rest entries are set to zero. To synchronize the dyn and stor arrays, we
use the collective communication routine MPI AllReduce with the operation MPI SUM that
sums up different copies of an array and sends the result back to all nodes (lines 20 and 23 of
Algorithm 6).
91
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
input : K 0
output : K, P
1 marker: if diff > tol and it < numIt then
2 for t ← 0 to T do
3 Calculate prices P (t) (inner loop);
4 Update dyn arrays;
5 end
6 Update stor arrays;
7 it ← (it + 1);
8 diff ← update (K, P );
9 end
10 goto marker
Figure 4. Sequential version.
input : K 0
output : K, P
1 marker: if diff > tol and it < numIt then
2 omp parallel default(private)
3 omp shared(dyn arrays, stor arrays)
4 omp copyin(parameters)
5 omp for
6 for t ← 0 to T do
7 Calculate prices P (t) (inner loop);
8 Update dyn arrays;
9 end
10 omp end parallel
11 Update stor arrays;
12 it ← (it + 1);
13 diff ← update (K, P );
14 end
15 goto marker
Figure 5. OpenMP version.
To calculate diff in the outer loop, we normalize the error using the variable norm. In the
sequential version this variable is calculated in the first step of the loop over time by one of the
functions. The variable normCounter is used to calculate the number of calls of this function.
For better comparison of with the sequential version, we keep the same way of calculating norm
in the MPI version. To this end, at each iteration of the outer loop, all nodes wait till the node
with rank zero is completed. The zero node sends norm and normCounter using MPI SEND
to all other nodes and they resume their work (lines 5 – 7 of Algorithm 6). This approach does
not affect the productivity much since most of the time is spent on collective operations (lines
20 and 23).
5. Results and discussion
The input data of the PET model consists of projections over time and base year data.
Projections are used for the population, i.e. the number of consumers in different groups, and
for productivity growth coefficients over the interval [0, T ]. The base year data divides into
production data from input-output tables, i.e. production shares, and consumer survey data,
i.e. capital, expenditures, savings, etc. of different consumer groups.
The calculations were carried out over T = 200 years with the time step of one year. The
tolerance is tol = 10−5 for the outer loop and tol′ = 10−6 for the inner loops. The total number
of production sectors is NX = 10. Then the number of variables is of the order T × NX ∼ 103 .
92
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
input : K 0
output : K, P
1 MPI INIT();
2 nsize ← MPI COMM SIZE();
3 rank ← MPI COMM RANK();
4 marker: if diff > tol and it < numIt then
5 if rank! = 0 then
6 MPI RECV for norm and normCounter from node rank;
7 end
8 for t ← (rank) to T by nsize do
9 Calculate prices P (t) (inner loop);
10 Update dyn arrays;
11 end
12 Update stor arrays;
13 it ← (it + 1);
14 diff ← update (K, P );
15 if rank = 0 then
16 for rankId ← 1 to (nsize − 1) do
17 MPI SEND norm and normCounter to node rankId;
18 end
19 end
20 MPI AllReduce(MPI SUM) for all dyn arrays;
21 end
22 goto marker
23 MPI AllReduce(MPI SUM) for all stor arrays;
24 MPI FINALIZE();
Figure 6. MPI version.
With a more detailed multisector production and several regions the number of variables can be
of order of hundred thousands.
For the OpenMP runs, we used Intel Visual Fortran Compiler XE 12.0. Calculations were
carried out on Intel Core i5 (4 cores, 2.70 GHz). For comparison, solution of the equations
system in the inner loop was obtained by Newton’s method, NITSOL and NEQBF. For each
method we made runs with optimization (-o3) and without it. The speedup graphs are presented
at Fig. 7. As we see, parallelization with four cores gives a more than three times speedup.
Calculations with the MPI version were carried out using Fortran compiler IBM XL Fortran
Advanced Edition. Calculations were carried out on Blue Gene/P V 11.1. Solution of the
equations system in the inner loop was obtained by Newton’s method and NITSOL. For each
method we made runs with optimization (-o3) and without it. The speedup graphs are presented
at Fig. 8. Here, parallelization with 200 nodes and one rank per node gives a speedup of up to
45 times.
The numerical results imply that the parallel version of the algorithm gives a substantial
improvement both with OpenMP and MPI. The efficiency of the parallel method is especially
pronounced with the MPI version, where the number of ranks can be as large as the number
of time steps. The increase of the speedup is almost linear and does not depend much on the
iterative method used in the inner loop and compiler optimization.
The results allow us to analyze the scalability. We see that the parallelization method
shows good strong scaling, i.e. the speedup increases with the number of processors for the fixed
problem size. This justifies the use of parallelization in the Fair-Taylor method.
Work is underway to determine the reasons why some of the methods used in the inner loop
perform better than others. Future work could address the effectiveness of the MPI version with
several ranks per node and hybrid realizations, such as OpenMP+MPI, with several processes
per node (see, e.g., [13, 14]).
93
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
4
Newton (no opt.)
Newton (-o3)
3.5 NITSOL (no opt.)
NITSOL (-o3)
NEQBF (no opt.)
NEQBF (-o3)
3
Speedup
2.5
2
1.5
1
1 2 3 4
Threads
Figure 7. The speedup for the OpenMP parallel version of the algorithm.
50
Newton (no opt.)
45 Newton (-o3)
NITSOL (no opt.)
40 NITSOL (-o3)
35
30
Speedup
25
20
15
10
5
0
1 40 80 120 160 200
Nodes
Figure 8. The speedup for the MPI parallel version of the algorithm.
94
Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org
References
1. Dalton M., O’Neill B., Prskawetz A., Jiang L., Pitkin J. Population aging and future carbon
emissions in the United States // Energy economics. 2008. Vol. 30, N. 2. P. 642–675.
2. Melnikov N.B., O’Neill B.C., Dalton M.G. Accounting for the household heterogeneity in
dynamic general equilibrium models // Energy economics. 2012. Vol. 34, N. 5. P. 1475–1483.
3. Fair R., Taylor J. Solution and maximum likelihood estimation of dynamic nonlinear rational
expectations models // Econometrica. 1983. Vol. 51, N. 4. P. 1169–1185.
4. Brayton F. Two practical algorithms for solving rational expectations models // Finance and
Economics Discussion Series. 2011–44, U.S. Federal Reserve Board.
5. Faust J., Tryon R. A distributed block approach to solving near-block-diagonal systems with
an application to a large macroeconometric model // Computational Economics. 1995. Vol. 8,
N. 4. P. 303–316.
6. Ortega J.M., Rheinboldt W.C. Iterative Solution of Nonlinear Equations in Several Variables.
SIAM, 2000. 572 p.
7. Eisenstat S., Walker H. Globally convergent inexact Newton methods // SIAM Journal on
Optimization. 1994. Vol. 4, N. 2. P. 393—422.
8. Pernice M., Walker H. Nitsol: A Newton Iterative solver for Nonlinear systems // SIAM
Journal on Scientific Computing. 1998. Vol. 19, N. 1. P. 302–318.
9. Kelley C.T. Iterative Methods for Linear and Nonlinear Equations. SIAM, 1995. 166 p.
10. Conn A.R., Gould N.I.M., Toint P.L. Trust-region methods. SIAM, 2000. 959 p.
11. Stokey N.L., Lucas R.E., Prescott E.C. Recursive Methods in Economic Dynamics. Harvard
University Press, 1989. 608 p.
12. Varian H.R. Microeconomic Analysis. Norton, 1992. 563 p.
13. Hasert M., Klimach H., Roller S. CAF versus MPI - Applicability of Coarray Fortran to a
Flow Solver // Springer. 2011. Vol. 6960, P. 228-236.
14. Hoefler T., Dinan J., Buntinas D., Balaji P., Barrett B., Brightwell R., Gropp W., Kale V.,
Thakur R. MPI+MPI: A New Hybrid Approach to Parallel Programming with MPI Plus
Shared Memory // Computing. 2013. Vol. 95, P. 1121-1136.
95