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.

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 x k = (x k 1 , . . . , x k n ) T and i − 1 components x k+1 1 , . . . , x k+1 i−1 of the (k + 1)th iterate x k+1 have been determined. Then to obtain x k+1 i the ith equation of the system

i−1 j=1 a ij x k+1 j + a ii x i + n j=i+1 a ij x k j = b i ,

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org is solved with respect to x i and the solution is taken as x k+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

x k+1 = M x k + c, where M = −(D + L) −1 U is the iteration matrix and c = (D + L) −1 b. The usual termination criterium is ||r k ||/||b|| < ε, where r k = b − Ax k is the residual. Since r k = (D + L)(x k+1 − x k ),

we can use the properly scaled error x k+1 − x k to terminate iterations.

The same idea is generalized to a nonlinear equations system f (x) = 0, where f = (f 1 , . . . , f n ). To obtain x k+1 i the ith equation

f i (x k+1 1 , . . . , x k+1 i−1 , x i , x k i+1 , . . . , x k n ) = 0

is solved with respect to x i and the solution is taken as x k+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 f i (x). Now we assume that f i and x i are vectors:

f i = (f i1 , . . . , f ij , . . . , f il ), x i = (x i1 , . . . , x ij , . . . , x il ),

where the blocks f ij and x ij are also vectors having the same dimension. Let the jth block of functions f ij depend only on the jth block of variables x ij . Then to obtain x k+1 i the equations

f ij (x k+1 1 , . . . , x k+1 i−1 , x, x k i+1 . . . , x k n ) = 0, j = 1, l,(1)

are solved independently with respect to x and the solution is taken as x k+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 x k , the linearized system

F ′ (x k )s k = −F (x k )(2)

is solved with respect to the step s k in order to obtain the next iterate x k+1 = x k + s k . In Newton's method, a LU -factorization of the Jacobian F ′ (x k ) is used (see, e.g., [6]). The inputs of the algorithm are the initial iterate x 0 and a termination tolerance tol. 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 ′ (x k )s k + F (x k )|| ≤ η k ||F (x k )||,(3)

where η k ∈ [0, 1) is the tolerance. To obtain the step s k , the generalized minimal residual (GM-RES) Krylov subspace method is used (see, e.g., [9]). Once an initial s k 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 x k

1 k ← 0; 2 while ||F (x k )|| > tol do 3 solve(2

) with respect to s k using a Krylov subspace method (termination criterium (3)); 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

4 while ||F (x k + s k )|| > [1 − t(1 − η k )]||F (x k )|| do 5 η k ← 1 − θ(1 − η k ); 6 s k ← θs k , θ ∈ [θmin, θmax]; 7 end 8 x k+1 ← x k + s k ; 9 k ← k + 1;||F (x k ) + B k s k || 2 − → min s k ,(4a)||s k || ≤ δ k .(4b)

is solved to get the direction s k , where δ k is the size of the trust region and B k is the approximate Jacobian evaluated at the current point x k . Then, the function at the point x k+1 = x k + s k is evaluated and used to decide whether the new point x k should be accepted. Problem ( 4) is solved approximately by the double dogleg method. The Jacobian is approximated by Broyden's formula

B k+1 = B k + ([F (x k+1 ) − F (x k )] − B k s k )s T k s T k s k . input : x0, tol; δ0, θ ∈ (0, 1) output : the solution x k 1 k ← 0; 2 while ||F (x k )|| < tol do 3 δ k ← δ0; 4 while f (x k+1 ) ≥ f (x k ) do 5 solve(

4) with respect to s k using the double dogleg method; Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org

6 δ k ← θδ k ; 7 x k+1 ← x k + s k ; 8 end 9 k ← k + 1;

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.

Optimal economic growth

The consumer side of the model consists of N c heterogeneous groups of households. Demand of the ith consumer group is determined by the intertemporal utility function

U i (c i ) = 1 ψ ∞ t=0 β t n it u i (c it ),

where t = 0, 1, ... is time, c it is the level of consumption (here and henceforth, small letters indicate the per capita values), β ∈ (0, 1) is the discount rate, ψ ∈ (−∞, 1) is the time substitution rate and n it is the size of population. The instantaneous utility function u i (c it ) is given by the constant-elasticity of substitution (CES) function

u i (c it ) =   Nc j=1 (µ ijt c ijt ) ρ   ψ ρ

, with constant electivity σ = 1/(1−ρ). Here, the index j = 1, N c labels goods, ρ ∈ (−∞, 1)\ {0} is the substitution parameter among goods and µ ijt is the preference coefficient. Capital dynamics is described as

(1 + ν it ) k i,t+1 = (1 − δ) k it + x it , k i0 > 0,(5)

where

x it is investment, δ ∈ (0, 1) is the capital depreciation coefficient, 1 + ν it = n i,p jt c ijt + q t x it = (1 − θ it ) ω t l it + (1 − φ it ) r t k it + g it ,(6)

where p jt is the price of jth consumer good, q t and r t are prices of investments and capital, ω t is the wage rate, g it is the government transfers, l it is the labor supply, θ it and φ it are the tax rates on capital and labor incomes. Variables k it , c ijt and x it are determined by maximizing

U i (c i ) − → max c it ,k it ,x it ,(7)

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 u i (c it ) at time t given the expenditures m it :

  Nc j=1 (µ ij c ij ) ρ   1 ρ − → max c ij , Nc j=1 p j c ij = m i ,p j c ij − → min c ij ,   Nc j=1 (µ ij c ij ) ρ   1 ρ = c i , we obtain min   Nc j=1 p j c ij   = pi c i , pi =   Nc j=1 p j µ ij ρ ρ−1   ρ−1 ρ

.

The second step is to solve problem ( 5)-( 7) for savings as a function of time (dynamics). Rewriting ( 5)- (7) in terms of pit and c it :

1 ψ ∞ t=0 β t n it c it ψ − → max c it ,k it , p it c it + q t x it = (1 − θ it ) ω t l it + (1 − φ it ) r t k it + g it , (1 + ν it ) k i,t+1 = (1 − δ) k it + x it ,

we obtain the first-order optimality condition (Euler equation)

q t p it c it ψ−1 = β q t+1 (1 − δ) + (1 − φ i,t+1 )r t+1 p i,t+1 c i,t+1 ψ−1 .(8)

The transversality conditions

lim t− →∞ λ it k it = 0,(9)

where λ it is the Lagrange multiplier, ensures the existence of the optimal trajectory (see, e.g., [11]).

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: N C consumer goods and "investment good", or an intermediate good: N E energy types and the rest, which is called materials (N X = N C +1+N E +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

X = γ X α K (G K K) ρ X + α L (G L L) ρ X + α Ē (G Ē Ē) ρ X + α M (G M M ) ρ X 1 ρ X ,(10)

where G I is the productivity factor for I = K, L, Ē, M and γ X normalizes α I to sum to unity. The parameters α I and G I can be sector and time dependent. The producer of the good X minimizes the costs

P K K + P L L + P Ē Ē + (1 + τ M )P M 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

P K K + P L L + P Ē Ē + (1 + τ M )P M M = P X X,

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org where the price P X is1

P X = 1 γ X α 1 1−ρ X K P K G K ρ X ρ X −1 + α 1 1−ρ X L P L G L ρ X ρ X −1 + + α 1 1−ρ X Ē P Ē G Ē ρ X ρ X −1 + α 1 1−ρ X M (1 + τ M )P M G M ρ X ρ X −1 .

The input-output coefficients A I = I/X for I = K, L, Ē are given by

A I = 1 α I (γ X G I ) ρ X P I P X 1 ρ X −1 ,

and for I = M by

A M = 1 α M (γ X G M ) ρ X (1 + τ M )P M P X 1 ρ X −1 .

Similarly, production of the energy composite Ē is determined according to

Ē = γ E i α E i (G E i E i ) ρ E 1 ρ E ,(11)

where

E i , i = 1, N E are different energy types. Solving the cost-minimization problem i (1 + τ E i )P E i E i − → min E i

under the given the level of production (11), we derive the prices

P Ē = 1 γ E i α 1 1−ρ E E i (1 + τ E i )P E i G E i ρ E ρ E −1 ρ E −1 ρ E

and input-output shares A E i = E i / Ē for the energy use:

A E i = 1 α E i (γ E G E i ) ρ E (1 + τ E i )P E i P Ē 1 ρ E −1 ,

where τ E i 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,

GREV t = N d i=1 n i (φ i rk i + θ i ωl i ) + N X j=1 τ X j X j + N E s=1 τ E S j P E S j E S j + τ M j M j = = N d i=1 (φ i P K K i + θ i P L L i ) + N X j=1 τ X j + N E s=1 τ Es j P Es j A Es j A Ēj + τ M j A M j X j .

Government expenditures are the sum of purchases and transfers, GEXP t = GP t + GT t . Government purchases are assumed to be constant in the current prices: GP t /p t = (1 + ξ) t GP 0 /p 0 , where pt = (1/N d ) N d i pit . It is assumed that GP 0 is given by a Cobb-Douglas production function, so that the input-output coefficients are

A GP K = α K P K , A GP L = α L P L , A GP M = α M P M , A GP E i = α E i P E i , i = 1, N E .

Similarly, government baseline transfers are neutral. There is additional lump-sum adjustment LSA t to balance the government budget:

GT t = pt p0 (1 + ξ) t N d i=1 n it g i0 + LSA t .

Equilibrium conditions

Aggregate supply for capital and labor at each point in time are determined by summing over levels supplied by each consumer group:

K AS = N d i=1 n i k i , L AS = N d i=1 n i l i .

Aggregate demand for capital and labor are

K AD = N X j=1 A j K X j + A GP K GP, L AD = N X j=1 A j L X j + A GP L GP.

Similarly, for consumer goods and investment.

Total demand for intermediate goods z

= X AD 1 , . . . , X AD N E +1

T 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

y =   N X N E +2 A j 1 X AD j + A GP 1 GP, . . . , N X N E +2 A j N E +1 X AD j + A GP N E +1 GP   T = (Y 1 , . . . , Y N E +1 ) T .

Since the production functions are homogeneous functions, equilibrium implies

X AD i = N E +1 j=1 A j i X AD j + N X j=N E +2 A j i X AD j + A GP i GP, i = 1, N E + 1.

or in the compact form z = Az + y, where A = (A j i ) is the (N E + 1) × (N E + 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 (N E + 1) × (N E + 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 , L AD = L AS , X AD = X AS , GREV = GEXP.

Note that the consumer prices p j and q are the gross (after-output-tax) prices P X 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 P K and labor P L in industry.

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 x it , k it , l it , etc. grow exponentially with the growth rate ξ. Proceeding to the limit t − → ∞ in (8), we obtain

(1 + ξ) 1−ψ = β 1 − δ + (1 − φ)r q .(12)

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).

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). 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).

Суперкомпьютерные дни в

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 N X = 10. Then the number of variables is of the order T × N X ∼ 10 3 .

Суперкомпьютерные дни в России 2015 // Russian Supercomputing Days 2015 // RussianSCDays.org 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]).