We present a parallel algorithm for computing an equilibrium path in a large-scale economic 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 improves 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 2. Block Gauss-Seidel method i−1 X aij xjk+1 + aiixi + j=1
n X aij xjk = bi, j=i+1 is solved with respect to xi and the solution is taken as xik+1. 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, rk = (D + L)(xk+1 − xk), where M = −(D + L)−1U is the iteration matrix and c = (D + L)−1b. The usual termination criterium is ||rk||/||b|| < ε, where rk = b − Axk is the residual. Since we can use the properly scaled error kxk+1 − xkk to terminate iterations.
The same idea is generalized to a nonlinear equations system f (x) = 0, where f = (f1, . . . , fn). To obtain xik+1 the ith equation
fi(x1k+1, . . . , xik−+11, xi, xik+1, . . . , xkn) = 0 is solved with respect to xi and the solution is taken as xik+1. 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 xik+1 the equations fij(x1k+1, . . . , xik−+11, x, xik+1 . . . , xkn) = 0, j = 1, l, are solved independently with respect to x and the solution is taken as xikj+1.
To solve the equations system for a separate block (
F ′(xk)sk = −F (xk) 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 (
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
where ηk ∈ [0, 1) is the tolerance. To obtain the step sk, the generalized minimal residual
(GMRES) 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 ∈ (
ηmax ∈ [0, 1), t ∈ (
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) + Bksk||2 −→ min,
sk
||sk|| ≤ δk.
(4a)
(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 (
Bk+1 = Bk + ([F (xk+1) − F (xk)] − Bksk)skT .
skT sk input : x0, tol;
δ0, θ ∈ (
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
where t = 0, 1, ... is time, cit is the level of consumption (here and henceforth, small letters
indicate the per capita values), β ∈ (
Ui(ci) = 1 X∞ βtnitui(cit), ψ t=0 Nc ui(cit) = X(μijtcijt)ρ 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,
where xit is investment, δ ∈ (
Nc X pjtcijt + qtxit = (1 − θit) ωtlit + (1 − φit) rtkit + git, 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.
under constraints (
max , cit,kit,xit −→ max, cij Nc X(μij cij )ρ j=1
1 ρ
Nc
X pj cij = mi,
j=1
(
Nc X pj cij −→ min, j=1 cij Nc X(μij cij )ρ j=1
1 ρ Nc min X pj cij = p¯ici, j=1
Nc p¯i = X j=1 pj μij = ci,
ρ−1 ρ−ρ1 ρ
The second step is to solve problem (
pit pi,t+1
lim λitkit = 0, 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 E¯ and materials M according to a CES function where GI is the productivity factor for I = K, L, E¯, 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 + PLL + PE¯ E¯ + (1 + τM )PM M −→
min
K,L,E¯,M
under the given level of production (
= PX X,
(
PX =
GE¯ PE¯ ρXρX−1 + α M1−1ρX (1 + τM )PM
GM ρXρX−1 ! ,
(
PE¯ = 1 γE
1 X α E1−iρE i
ρE ρE−1 (1 + τEi)PEi ρE−1 ! ρE
GEi and input-output shares AEi = Ei/E¯ for the energy use:
AEi =
1 αEi(γEGEi)ρE
PE¯
1 (1 + τEi)PEi ρE−1 , 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 GREVt = X ni(φirki + θiωli) + X i=1 j=1
NE τXj Xj + X τESj PESj ESj + τMj Mj s=1 ! =
Nd NX = X(φiPK Ki + θiPLLi) + X i=1 j=1
NE τXj + X τEsj PEsj AEsj AE¯j + τMj AMj s=1 !
Xj.
Government expenditures are the sum of purchases and transfers, GEXPt = GPt + GTt. Government purchases are assumed to be constant in the current prices: GPt/p¯t = (1 + ξ)tGP0/p¯0, 1By the envelope theorem, the price is identified with the marginal cost of production (see, e.g., [12]). where p¯t = (1/Nd) PiNd p¯it. It is assumed that GP0 is given by a Cobb-Douglas production function, so that the input-output coefficients are
AGKP = αK
PK ,
ALGP = αL
, PL
AGMP = αM ,
PM
AEGiP = αEi ,
PEi i = 1, NE.
Similarly, government baseline transfers are neutral. There is additional lump-sum adjustment
Nd GTt = p¯t (1 + ξ)t X nitgi0 + 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 KAS = X niki, i=1
Nd LAS = X nili.
i=1
NX KAD = X AjK Xj + AGKP GP, j=1
NX LAD = X AjLXj + ALGP GP. j=1
Total demand for intermediate goods z = X1AD, . . . , XNAED+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
NX NX y = X Aj1XjAD + A1GP GP, . . . , X
NE+2 NE+2
j ANE+1XjAD + AGNPE+1GP
= (Y1, . . . , YNE+1)T .
T T Since the production functions are homogeneous functions, equilibrium implies
NE+1 XiAD = X AijXjAD + j=1
NX
X j=NE+2
AijXjAD + AiGP GP, i = 1, NE + 1. or in the compact form z = Az + y, where A = (Aij) 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)−1y, 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: KAD = KAS,
LAD = LAS,
XAD = XAS,
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. 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 (
q (1 + ξ)1−ψ = β 1 − δ + .
(
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 Ks has been determined. Then solving the system f (K, P ) = 0, f (Ks, 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 Ks+1 = K. The part of the algorithm that calculates Ks+1 from
4. Algorithm implementation
First, we describe the sequential version (Algorithm 4). The input is the initial guess of the capital K0 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). input : K0 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 input : K0 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
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. input : K0 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 MPI AllReduce(MPI SUM) for all dyn arrays;
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]). 2
3
1.5 1 1 p u d ee 25 p S 50 45 40 35 30 20 15 10 5 0
NITSOL (-o3) 1 40 80 120 160 200