ON THE PARALLEL SOLUTION OF PARABOLIC

EQUATIONS

E. GALLOPO ULOS

Y. $AAD

May, 1989

Research Institute for Advanced Computer Science

NASA Ames Research Center

RIACS Technical Report 89.19

NASA Cooperative Agreement Number NCC 2-387

(NASA-CR-16036_) ON THr PAPALLFL SOLUTION

OF P^R^_'_' [C:. ._ _ F_UATInN_. (Research Inst. for

Adv_anced ConlputoF Scionce) 24 p CSCL 12A

NO0-24009

G_lo4

Research Institute for Advanced Computer Science
- II
ON THE PARALLEL SOLUTION OF PARABOLIC

EQUATIONS

E. GALLOPOULOS

Y. SAAD

May, 1989

Research Institute for Advanced Computer Science

NASA Ames Research Center

RIACS Technical Report 89.19

NASA Cooperative Agreement Number NCC 2-387

ON THE PARALLEL SOLUTION OF PARABOLIC

EQUATIONS

E. GALLOPOULOS

Center for Supercomputing Research and Development

University of Illinois at Urbana-Champaign

Urbana, Illinois 61801

stratis@uicsrd, esrd .uiuc. edu

Y. SAAD

RIACS, Mail Stop _30-5

NASA Ames Research Center

Moffet Field, California 9_035

saadCriacs, edu

ABSTRACT

We propose _new parallel algorithms for the solution of linear parabolic problems.

The first of these methods is based on using polynomial approximation to the

exponential. It does not require solving any linear systems and is highly parallelizable. The two other methods proposed are based on Padd and Chebyshev

approximations to the matrix exponential. The parallelization of these methods

is achieved by using partial fraction decomposition techniques to solve the resulting systems and thus offers the potential for increased time parallelism in time

dependent problems. We also present experimental results from the All]ant FX/8

and the Cray Y-MP/832 vector multiprocessors.
1. Introduction. We consider the following linear parabolic partial

differential equation:

(1.1) au(z, t)

: Lu(_,t) + _(_), • c n

0t

u(O,_:) = uo, Yz cft

,,(t, _) = ,,(_), • c on, t >_0.

where L is a second order partial differential operator of elliptic type, acting

on functions defined on the open bounded set ft. If the method of lines is

used to solve (1.1), then this partial differential equation is first discretized

with respect to space variables, resulting in a system of ordinary differentiM

equations of the form

dw(t)

- Aw(t) + r

dt

w(O) : _o

whose solution is explicitly given by

(1.2) w(t) = A-lr + e-tA(wo -- A-'r)

which simplifies to

(1.3) w(t) : e-'Aw0

in the case of a homogeneous system (r = 0). Note that if we denote by

(v(t) =_ w(t) - A-lr and accordingly, _b0 -= w0 - A-lr, then _b(t) satisfies e_

homogeneous system and

(1.4) (o(t)-- e-tA@o

Thus, if we want to obtaip, the solution at time_ t in one sinai e stele., we

would need to operate with the matrices A -1 and e -tA on cert_n initial

vectors. This solution faces the following difficulties:

• Computing e-tAtbo may not be easy, especially for large t.

• The cost of computing A-lr is not negligible for more than one

space dimensions.

• In many problems the operator L, as well as the forcing term s,

may vary with t; If this variation is rapid the above formula is not

applicable or may be very inaccurate for large t;
Note that if we denote by f the function f(z) -- (1 - e-Z)/z we can

rewrite the solution (1.2) as follows,

w(t) = e-tAwo + f(tA)tr

(1.5) = Wo + tf(tA)[r - Awo]

This removes the term A-lr from the expression (1.2), at the expense of

dealing with the function f(z) instead of e -z [22]. The above expression

also shows more clearly the dependence of the solution with respect to the

initial condition w0 and the forcing function r.

Assume now that instead of attempting to compute the solution at time

t in one single step, we use a time-stepping procedure. At time t + At,

the solution will depend on w(t) which plays the role of w0 in the above

expressions and we get

(1.6)

from (1.4), or

w(t + At)

e(t + At) =

= e-atAw(t)+ Atf(AtA)r

(1.7) = w(t) + Atf(AtA)[r-- Aw(t)].

from (1.5).

We should observe that the use of the variable tb in formula (1.6) requires

computing A-lr only once and not at every step of the stepping procedure.

The advantage of using (1.7) over (1.6) is therefore limited, except when A

and varies with time.

In both (1.7) and (1.6), we need to compute a vector of the form

q(AtA)v, where q(z) is a known analytic function in z. The basic idea

for computing (1.6) and (1.7), is to find a suitable approximation g(A), to

the function q(A) and then substitute this approximation in (1.6) or (1.7).

This is complicated by the following facts. First, depending on the operator

L, the type of discretization performed and the boundary conditions, A may

be symmetric positive definite, or nonsymmetric. It may also be singular or

nearly singular. Moreover, in typical problems A is large and sparse making

the direct calculation of q(AtA) by usual methods prohibitive.

Note that there is no need to actually evaluate the matrix g(AtA). Instead all we need is be able to evaluate g(AtA)v inexpensively for some vector

v. In this paper we show how to do this for the specific case q(z) = e -z which

allows to solve the general problem via (1.6). For notational convenience,

we will assume that after suitable scaling At = 1.

There are two different ways of generating approximations g(A)v. The

first is by using polynomials. The resulting procedure will only require matrix by vector multiplications with the matrix A, and is therefore very easily

3
paxallelizable.In Section3, we will consider an approach of this type based

on using a Krylov subspace technique.

The second class of methods consists of taking g to be a rational function

of the form r(z) = p(z)/q(z). In this situation, a difficulty arises when computing r(A)v on parallel machines. Typically the denominator is factored

as q(z) -: I-I___l(z- )_i) and q(A)-lv is computed by solving the sequence

of linear systems (A - AiI)ui : ui-1 with u0 = v. This is a sequential

process which can be particularly damaging, especially for one-dimensional

problems where A is usually a tridiagonal matrix. Although there are many

efficient methods for solving tridiagonai systems, it is clear that if we have to

solve only one single system per step, we will very likely be under-utilizing

the computational resources. To circumvent this sequential constraint we

propose in Section 4 to use the partial fraction expansions of r(z). This will

transform the sequential solution of (A - ._I)ul = ui-1 into solving in parallel the independent linear systems (A - )_iI)ui = v and then taking a linear

combination of the results ui. The advantages for one-dimensional problems

are clear. For higher-dimensional problems, this allows to remove the need

to parallelize each of the linear systems (A)_iI)ui = ui-1. In effect it offers

a means for achieving parallelism in an extremely simple manner, far simpler than would be needed in optimizing the linear solves in the traditional

approach. This is achieved by using high order schemes, i.e., high degree

rational approximations. As a result an added benefit is that the overall

amount of work is also reduced. As is stated in our conclusion, it seems

that high order integration schemes in ODE methods offer a tremendous

potential in a parallel processing environment.

2. Previous work. The previous discussion underscores the direct connection that exists between the topic of this paper and that of of parallel

solution of ordinary differential equations (ODEs). As argued in [8], the most

important situation when considering parallel methods for solving systems of

ODE's is when the problem is very large, as is the case for systems resulting

from a Method of Lines semi-discretization of a partial differential equation

such as (1.1). We refer the reader to [13] for methods to approximate the

matrix exponential, to [18,21] for polynomial approximations in parabolic

problems, and to the work of Varga and co-authors for rational approximations ([22,3,2,12]). In [7,20] a method was introduced for the parallelization

of Block Cyclic Reduction (BCR). The connection between the method considered in these papers and the question addressed here, is that when using

BCR one must e;caluate a vector of the form q(A)v, where q is a rational

function. Partial fractions in a sequential context for time dependent problems were used in [12,19,26] and suggested in parallel complexity studies in

[11]. Finally recent experiments of Reusch et al. ([17]) demonstrated re-
markablegainsin efficiencyand accuracyfor the solutionof homogeneous

linearevolutionequationsbymeansof high-orderdiagonalPad6approximations. As will be arguedin Section4 suchschemesareextremelyattractive

onparallelmachines,whenproperlyimplemented.

8. Polynomial approximations. In this section we consider using

polynomial approximation to the exponential, i.e., we seek an approximation

to e-Av of the form

(3.1) e-Av '_ pm_l(A)v

where pro-1 is a polynomial of degree m - 1. The main attraction of polynomial based techniques is their explicit nature, that is the fact that they do

not require solving linear systems. In fact the only operations required with

the matrix A are multiplications with vectors, an operation that is very easy

to parallelize and vectorize. On the other hand polynomial approximation

cannot handle very stiff problems as well as rational approximations. As a

result the trade-off is a large number of matrix by vector multiplications versus no linear systems to solve. For two-dimensional and, more importantly,

for three-dlmensional problems polynomial based schemes, if implemented

with care can be very attractive.

There are several ways in which polynomial approximations can be

found. The simplest technique is to attempt to minimize some norm of

the error e -z - p,-,,_l(z) on a continuum in the complex plane that encloses

the spectrum of A. For example, Chebyshev approximation can be used.

The disadvantage of this is that it requires some approximation to the spectrum of A. In this paper we consider only approaches that do not require

any information on the spectrum of A.

The approximation (3.1) to e-'4v is an element of the Krylov subspace

Km= span{v, Av,..., A'_-Iv}.

In order to manipulate vectors in Km it is convenient to generate an orthonormal basis V,_ = [vl,v2, v3,...,vm]. We will take as initial vector

vl = v/[[vii2 and generate the basis V,_ with the well-known Arnoldi algorithm, described below.

5
Algor|thm- Arnoldi

1. Initialize:

Compute vl := v/llvll_.

2. Iterate: Do j = 1, 2, ..., m

1. Compute w := Avj

2. Compute a set of j coefficients h 0 so that

J

(3.2) W := W -- E hijvi

i=1

is orthogonal to all previous vi's.

3. Compute hj+,j -- Ilwll2and vj+l = w/hj+xd.

By construction the above algorithm produces an orthonorma! basis

V,_ = [vl, v2,...,vm], of the Krylov subspace Kin. If we denote the m x m

upper Hessenberg matrix consisting of the coefficients hij computed by the

algorithm by Hm we have the relation

(3.3) AVm = VmHm + hrn+l,rnVm+le T

from which we get Hm = V_AVm. Therefore Hm represents the projection

of the linear transformation A to the subspace Kin, with respect to the basis

Ym.

We can write the desired solution x = pm_l(A)v as x = V,_y where y,

is an m-vector. Ideally we would like to minimize the norm lie-Up - VmyH2.

The solution to this optimization problem is known to be

(3.4) uo ,: vge-Av

Unfortunately, this is not computable because it involves the unknown vector

e-av. However, if we assume that vx = j0v then we have y_t = flVTe-AVmel

and it is natural to approximate VTe-AVm by e -n'', leading to the approximation,

(3.5) e-av ._ _Vme-H"el

This immediately raises a question concerning the quality of this approximation. A first observation is that the above approximation is exact

for m -- n. This is because in this situation vm+l = 0 and (3.3) becomes

AVm = VmH_, where V,, is an n × n orthogonal matrix. In fact, similarly to

the conjugate gradient method and the ArnoIdi process, the approximation

will be exact for m whenever m is larger or equal to the degree of the minimal polynomial of vt with respect to A. This however, is unlikely to take

6
placebeforem = n. More generally, the following theorem provides a rough

bound on the error and establishes convergence when m increases.

THEOREM 3.1. Let A be any matriz and let p : IIAII2. Then the error

of the approzimation (3.5) is such that

(3.6) lie-% ll2<_2 P-7 p.

The proof of this result is onfitted. This result as well as sharper bounds

will be fully discussed in a forthcoming paper.

The theorem shows convergence of the approximation (3.5). It can also

serve as a guide to choosing the step size in a time-stepping procedure.

Indeed, if we were to replace A by the scaled matrix rA, then the Krylov

subspace will remain the same, i.e., V,n will not change, and Hm will be

scaled to rH,_. As a result the bound (3.6) becomes,

(3.7) IIe- av - II<-

The consequence of {3.7) is that by reducing the step-size one can always

make the scheme accurate enough, without changing the dimension rn.

We note that the idea of exploiting the Lanczos algorithm to evaluate

terms of the exponential of Hamiltonian operators has been extensively used

in Chemistry ([25]). The work in [15] is also related wherein systems of

ODEs are solved by first projecting into Krylov subspaces and then solving

reduced tridiagonal systems of ODEs.

So far we have not considered the important particular case where A is

symmetric. As is well-known in this situation Arnoldi's algorithm simplifies

into the Lanczos process, which entails a three-term recurrence. This is a

result of the fact that the matrix Hm = VTAVm must be symmetric and

therefore tridiagonal symmetric, and so all hi,j = 0 for i = 1,2,..,j- 2.

However, the resulting vectors which are in theory orthogonal to each other,

tend to loose their orthogonality rapidly.

From the practical point of view several problems must be addressed.

For example we can mention the following issues:

1. How should one compute the vector eH"el ?

2. In the case where A is symmetric, should orthogonality be enforced?

We note that for (1) we can use the methods described in the next sections

efficiently since Hm is either tridiagonal or Hessenberg. If m is small enough

as is the case in practical situations_ then the cost of computing eH'_ex will

be negligible. An important observation here is that since Vm is orthogonal,

the integration scheme based on the formula (3.5) is likely to inherit the

stability properties of the scheme used in approximating e-lt"'e 1. For this

reason it is crucial to use rational approximation.
For (2), if the matrix is nonsymmetric it is recommended to perform a

modified Grala-Schmidt process with partial reorthogonalization. Selective

or partial reorthogonalization can be used when A is symmetric [16].

4. Rational approximations.

4.1. Overview. As was mentioned earlier, a popular way of computing

approximations to e-Av is via a rational approximation to the exponential,

i.e.,

(4.1) e-'4v ._ q,.(A)-lP,_(A)v

The simplest approximation of this type, referred to as Pad_ approximation

can be found by matching the Taylor series expansion of the left-hand-side

and right-hand-side of (4.1) at the origin. This approximation is local, i.e.,

it is very accurate near the origin but may be inaccurate far away. For

this reason schemes based on more global approximation have been devised

[22,2]. Thus, for typical parabolic problems that involve a second order

partial differential equation L that is self-adjoint elliptic, the eigenvalues of L

are located in the interval ( -o¢_, 0) and it is therefore natural to seek the best

rational approximation to the function e z on this interval, or equivalently to

the function e -_ on the interval (0, _).

One of the main reasons why rational approximations have been preferred to polynomial approximations, is the better stability properties of the

corresponding integration schemes. Thus, it is known that the Pad_ approximations to the matrix exponential give rise to unconditionally stable

methods if and only if the degree m of the numerator is at most equal to the

degree of the denominator ([22]).

By far the best known rational approximation to the exponential is the

( 1, 1) Pad_ approximation e z _ ( ! _- 1⁄2z) / ( 1 - 1⁄2z) which when used in conjunction with (1.6) leads to the well-known Crank-Nicolson scheme. However,

because of its modest accuracy, there are limitations as to how large the step

si_e At can be and there might be, some large number, say, rnl, applications

of formula (1.6) before the solution at the final time T is found. At the other

extreme assume that one can find a highly accurate rational approximation

that allows to compute w(T) in just one application of (1.6). If the rational

approximation is of the type (m2, m2) then it is very likely that m2 << ml,

meaning that the total amount of work is far less with the more accurate

scheme. Thus, the more accurate schemes have tremendous potential in the

context of parallel processing precisely because of this feature. By their very

nature, low order schemes do not allow for much work to be shared at every

step while, as will be seen in the next section, high order schemes are easily

and safely parallelizab!e.

8
4.2. The use of partial fraction expansions. Let the rational approximation to the exponential of e -z be of the form

(4.2) qr(z)

where we assume that m < r. Then, at each application of the scheme

corresponding to (1.6) we would have to evaluate the vector

(4.3) z =p,_(A).q,(A)-lv

There are several ways in which one can compute qr(A)-lv. One economical

procedure involves factoring the polynomial q, (z) as

(4.4)

i=1

and then solving the successive 1/near systems

(4.5) (A - AJ)ui : ui_l,i = 1,2,...,r

with u0 = v. The final result is the desired solution. One then needs to

multiply the result by p,n(A). In fact, several modifications to the aforementioned scheme have been proposed in the literature to avoid this last extra

step [4]. Incidentally, we should mention that partial fraction expansion formulations can be used to explain many of these efficient implementations of

time stepping procedures. This is discussed in the note [6].

Clearly, a significant difficulty with the use of (4.5) is that it is a sequential process. System number i must be solved before system i + I since its

solution will be the right hand side of the next system.

An alternative approach used in [7] in a different context is to resort to

the partial fraction expansion of (4.2), namely,

(4.6) R,,,,_(z) = mrm + _ ai

where

w( 4 (4*) -

and 7r,,_ and r_ are the leading coefficients of the polynomials Pm and qr

respectively.

With this expansion the algorithm for computing (4.3) becomes:
Algorithm:

1. For i --1,2,...,rsolve(A - Ail)xi--v in parallel.

2. Compute z = -_v + _r=l Cq_i.

Note that forthe usualapproach, schemes with repeated poles(e.g.the

restrictedschemes [14])have oftenbeen preferredbecause they involvefewer

factorizatlonswith the standard techniqueswhen directmethods are used.

These factorizationsare very expensive for 2-D and 3-D problems. These

schemes cannot be used with our approach sincewe need to have distinct

poles.

As we can see,the method isvery wellsuitedforsystems offeringhie1"-

archical parallelism realized with multicluster architectures ([10]). For problems in two and three space dimensions, we could use a rational approximation generating as many independent systems as there are clusters. The

solution of each system could then proceed independently in each cluster.

We thus see the interesting phenomenon that not only numerical considerations but also the amount of available parallelism will drive the choice of

the order of the approximation.

4.3. Padd Approximation. in this section we outline the procedure

using Padd approximation to the exponential. Given the degrees of the

numerator and denominator, it is easy to automatically generate the rational

function as a pair of two polynomials both given in power form. More

precisely, the coefficients rrj,j = 0, ..., m of the polynomial p,, and _j,i =

0, ..., r of the polynomial qr are explicitly given by [22]:

= (-1)# (r+m-j)!m! (,. + - j)!'J. 0,.. @ m

(m + r - j)!r!

(4.8) = (m + ,)!2(,- J)!'J = 0,...r.

Then we need to compute the roots of the denominator. This we do by

some standard polynomial rootfinder. Once the roots are computed one

then needs to compute the coefficients ai of the partial fraction expansion

(4.6) using formula (4.7). For high degree polynomials, numerical difficulties

both in evaluating accurately the roots and in computing oi by formula (4.7)

are to be expected.

4.4. Chebyshev Approx;mation. When using Chebyshev approximation, one must first decide on which region the best rational function

nmst be computed. In this paper we only consider the best uniform approximation to e -z for z E (0, oo). Unlike the Padd case, the coefficients of the

numerator and denominator polynomials are not available analytically and

must be computed as the solution of an optimization problem. This can lead

10
to a fairly involved procedure, requiring the use of a Remez type algorithm.

We preferred instead to use directly the very accurate values from [1], where

the polynomial coefficients are provided for up to tile degree 30. Once these

coefficients are input, we proceed as before, calling a polynomial rootfinder

and evaluating the partial fraction coefficients.

The big advantage of Chebyshev methods, stressed in the work by Varga,

is the ability to use large step size. In fact, when A is Hermitian, a relation

of the form

[Iw(t)- w,,,,lh <_ Am,_Nw0 - A-lr[l_ Vt >_ 0

holds, where win,, is the solution computed with an (m, r) order Chebyshev

approximation, and Am,_ are constants converging to zero geometrically (see

[3] and [1] for a list of A,,,'s).

Although such a technique will have its limitations for time-varying coefficients and boundary conditions, it can often produce excellent results as

is demonstrated in our experiments. Thus by combining the large time step,

together with the problem decoupling for parallelism, we obtain a very efficient procedure in the sense that fast convergence, low error and efficient

exploitation of the parallel resources are achieved.

4.5. Handling complex poles. The partial fractions in the decomposition could involve complex shifts of the operator A. These complex shifts

come in conjugate pairs and correspond to complex poles of the rational function under consideration, that is the roots of the (real) denominator. Similar

problems occur elsewhere in linear algebra, e.g. in the course of the QR algorithm [24, Section 41]. Since the coefficients of the corresponding principal

parts in the partial fraction expansion also have conjugate coefficients, complex arithmetic can be avoided completely by writing the rational function as

a sum of fractions whose denominators can contain quadratic factors. Each

quadratic factor corresponds to a product of the form (z - _i)(z - )_i) for all

roots A of the denominator having non-zero imaginary parts. This is just a

case of incomplete partial fraction decomposition (see [9, Section 7.1]).

Some drawbacks to this technique are the need to form A 2, the extra

computations due to the first order expansion coefficients and the need to

store the quadratic factors.

In case A is banded, however, the formation of A s means a doubling of

the bandwidth, and will result in an increase of data locality when working

with A. If the corresponding matrix operations are designed carefully (e.g.

using blocking as is done in BLAS3) increased efficiency will result on architectures with hierarchical memories. A numerical drawback is due to the

squaring of the (possibly already large) condition of the matrix factors with

the ensuing drawbacks in the application of iterative methods.

11
A simplerway of dealing with complex poles is to observe that the

expansion coefficients oi associated with two complex conjugate pairs must

be conjugate. Then we can write

(4.9) oti(A - AiI)-I_: + _i(A - _iI)-l:r : 2_[ai(A - )_iI)-:z]

This requires solving one complex system as opposed to two. It has the

advantage of requiring less storage and fewer arithmetic operations than with

the squaring approach. Moreover, data locality is also preserved through the

use of complex arithmetic.

The above discussion addresses only the use of direct solvers. For 2-D

and, more importantly, for 3-D problems, iterative procedures become attractive and we would like next to discuss how complex poles can be handled

in this case. The first observation to be made is that we can again exploit the

fact that the poles usually come in conjugate pairs. Thus, we can use the conjugate gradient technique to solve a system of the form (A-AI)(A-_I);r = f,

which will involve no complex arithmetic in the CG iteration. Indeed, the

only operations that are needed with this matrix are matrix by vector mnltiplications of the form w = (A - )_I)(A - _I)v which can be performed inexpensively in real arithmetic when v is real as w = [AJZv + A(A - 2_(A)I)v.

Moreover, the storage requirement is also not affected since only A is needed.

Note that this is not equivalent to the normal equations approach. Preconditioning can also be easily retrofitted in this scheme. Indeed, the ICC(0)

preconditioners require only an extra diagonal of data. This extra diagonal

is complex and can be easily constructed. Using extra fill-in is, however_

troublesome since all of L and U matrices must be treated as complex. Once

the preconditioning M = LU has been built then the CG iteration can be

performed with the matrix M-I(A - )_I)lYI-I(A - _I) in real arithmetic.

We should note that the scheme described here represents the simplest, certainly not the best, of a mmlber of possible options. In particular, there are

methods which will not be described here, that do not involve the matrix

l)_12v + A(A - 2_()_)I) but the original matrix A. We also mention that the

problem of solving complex linear systems of the form (A - )_I)z = f has

been addressed by Freund who devised special iteration schemes [5].

5. Numerical experiments. In this section we will describe a few

tests to illustrate the behavior of the schemes described in this paper. In

particular we do not compare here the polynomial approach with the rational approximation approach. Further experiments will be presented in a

forthcoming report. Except where noted, our tests were conducted on an

Alliant FX/8 vector multiprocessor using 64 bit floating-point arithmetic.

The test problem is issued from the discretization of

ut -- uzx, x E (0,1)

12
o) = u(t, 1) = o

using 22 grid points, yielding a matrix of size n=20. The initial conditions

are chosen once the matrix is discretized, in such a way that the solution is

known for all t. More precisely,

_ I . jk_r (5.1) u(0, zj) = _ sm

k:l '/2 + 1

Note that the vector {sin jk'ra n+l Sj=l,..,n is an eigenvector of the discretized

operator. In order to decouple from the influence of errors due to spatial

discretization, for all experiments in this section we take the solution of the

semi-discrete problem ut = -Au to be the true solution.

We begin by illustrating the behavior of the polynomial approximation

to the exponential. First, we would like to show how the accuracy of the

polynomial scheme varies as m varies but At is fixed. We take At = 0.01

and let m vary from 1 to 20. The infinity-norm of the error between the

exact result e-aatWo and the approximation obtained from using the Arnoldi

process as described in Section 3 is computed and scaled by the infinity norm

of the exact solution. These relative error norms are then plotted in Figure 5

versus the subspaee dimensions m. Notice that as is expected for m = 20

the error is zero up to the machine accuracy and the errors induced by

the computation of e-He1, since the approximation (3.5) is exact in this

situation. In the tests dealing with polynomial approximation, the vector

e-AtH"ex is computed to very high accuracy, by compounding Taylor series

expansions of degree 10. The composition is done by scaling H,, by a scalar

6 in such a way that the spectral radius of 6Hm is less than 1/2. Thus, the

evaluation of e-AtFt'el may require a large number of successive evaluations

of vectors of the form e-6atH'e. A more elaborate implementation using

rational approximations as suggested earlier is under way.

Next we fix the dimension of the Krylov subspace to m = 10 and let At

vary from At = 0.005 to At = 0.1 with an increment of 0.005. The relative

error norms are plotted in Figure 2. Notice how the accuracy deteriorates

at once instead of progressively.

We now consider a three-dimensional version of the previous test problem, i.e., we discretize the problem

ut = u:_+u_u+u_, z,y, zC (0,1)

u = 0 on the boundary

using 17 grid points in each direction, yielding a matrix of size N = 153 =

3375. The experiment we now describe was performed on a Cray Y-MP/832

13
_3

S

-2

-4

-6

-8

dO

-12

-14

-160-_ _ -

r

x

o

i

x

x

o

o

x

x

o

o

x

z

x

o

o

o

x

z

][

l

o I

I

o

1

o

o

o

o

o

4 6 8 lO 12 14 16 18 20

Subspaccorder rn

FIG. !. Behavior of polynomial approximation for At = 0.005 ('o) and At = 0.01 (z) as

the degree m varies.

(8 processors). Again the initial conditions are chosen once the matrix is

discretized, in such a way that the solution is known for all t. We take

(5.2) zk) :

1 ii'rr . jjtx , kk'Tr i' + j' + k' sin _sm _sm n+ 1

iJ,j_,kr=l

The above expression is simply an explicit linear combination of the eigenvectors of the discretized operator.

The purpose of this experiment is to iliustrate the efficiency of using

high accuracy schemes versus low accuracy schemes. This point was stressed

earlier and constitutes one of the main motivations for this paper. As will be

seen later the same conclusions also hold for the methods based on rational

approximations to the exponential.

Assume that we want to integrate the above equation between t=0 and

t=0.1, and achieve art error-norm at t = 0.1 wtfich is less than e = 10 -l°.

Here by error-norm we mean the 2-norm of the absolute error. We Can

14
"d

-0.5

-1

-1.5

-2

-2.5

-3

-3.5

-4

-4.5

o

-5

0

o

o

o

o

i--7 i i 1

oooo o°°

o

o

i i.__. _ 1 .t ..x___ .... i i

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Timestep size

FIG. 2. Behavior of polynomial approximation for ra = 10 as timestep At varies.

vary both the degree m and the time-step At. In a normal procedure we

would first choose a degree m and then try to determine the maximum At

allowed to achieve the desirable error level. However, for convenience, we

proceed in the opposite way: we first select a step-size At and then determine

the minimum m that is needed to achieve the desirable error level. Here

the vector e-ar'_el was computed via the diagonalization of the tridiagonal

matrix Hm using EISPACK's routine IMTQL2. This is clearly not the most

efficient technique since/arm is tridiagonal. What is shown in Table 1 is the

various time steps chosen (column 1) and the minimum values ofm (column

2) to achieve an absolute error less than e z 10 -l° at t=0.1. We show in the

third column the total number of matrix-by-vector multiplications required

to complete the integration. The times required to complete the integration

on a Cray Y-MP are shown in two parts in columns 4 and 5. Since we used

an inefficient algorithm to compute e-H'_el we showed the total time for

performing this operation (denoted by Timem) separately in column 5. The

remaining time denoted by Timely and shown in column 4, represents the

time for performing the Arnoldi process and the linear combinations of the

15
TABLE 1

Performance of the polynomial scheme with varying accuracy on the Cray Y-MP/832.

At rn I

0.5000E-04 ,_ ,

0.1000E-03

0.5000E-03 10 I

0.1000E-02 _._ ,'

0.5000E-02 20 I

--- 4

0.1000E-01 26 I

0.2000E-01 33 [

0.3000E-01 391

0.4000E-01:44 ]

0.5000E-01 49 ]

O.lOOOE+OOL691.

M-vec's TimeN Wimern

12006 0.8264E+01

7007 0.4779E+01

2010 0.1329E_01

1200

llE,'ro ll2

0.3840E+00 0.3808E-10

0.2459E+00 0.1298E-10

0.9062E-01 0.1338E-10

0.7968E+00 0.6331E-01 0.1946E-10

400 0.2593E+00 0.3479E-01 0.7336E-10

260 0.1646E+O0 0.2939E-01 0.5304E-10

165 O.103OE+00 0.2487E-01 0.9857E-10

156 0.9504E-01

!32 0.8154E-01

98 0.5879E-01

69 0.4134E-01

O.2856E-01 0.6247E-10

0.2847E-01 0.4098E-10

0.2446E-01 0.5787E-10

0.3170E-01 0.7494E-10

Arnoldi vectors. Since for all tests m _< 69, and Hm is tridiagonal symmetric,

one can expect that with an efficient algorithm the total time for evaluating

the vector e-R"el will represent a small portion of the total execution time,

given the size of this problem. However, as is indicated by the last entry

of the table, this timc may become nonnegligible compared with the time

Timeg as m increases, if an inefficient algorithm is used, even though the

total number of operations involved is much smaller than that in the rest

of the computation. Another point is that the matrix is symmetric, so we

have used a Lanczos algorithm to generate the v_s instead of the full Arnoldl

algorithm. No reorthogonalization of any sort was performed. The matrix

consists of 7 diagonals, so the matrix by vector products are performed by

diagonals resulting in a very effective use of the vector capabilities of the

YMP. Based on the time Timeg for the last entry of the table, we have

estimated that the average Mflops rate reached (excluding the calculation

of exp{-H,_}ei) was around t61. This is achieved with virtually no code

optimization.

Note the very rapid decrease in the total number of matrix by vector

products required. The ratio between the lowest degree m --- 6 and the

highest degree m = 69 is 174. The corresponding ratio between the two

times is roughly 200. The case m = 69 can achieve the desired accuracy in

just one step, i.e., with At = 0.1. On the other hand for m = 6 a timestep of At : 0.0(}005 must be taken resulting in a total of 2000 steps. We

should point Out that we are restricting ourselves to a constant time-step,

but more efficient variable time stepping procedures are likely to reduce

the total nmnber of steps needed. From the result of Theorem 3.1, these

16
observations come with no surprise. In effect, increasing the dimension of

the Krylov subspace, will increase the accuracy in such a way that a much

larger p, i.e., a larger At, can quickly be afforded.

We next test the rational approximations described in Section 4. The

program asks for the type (Pad6 or Chebyshev) and order of the diagonal rational approximation. In the Pad6 case, the coefficients of the numerator and

denominator polynomials in the rational function are numerically evaluated

from (4.8). In the Chebyshev case the coefficients are taken directly from [1].

Subsequently, the IMSL routine ZRPOLY (based on the Jenkins-Traub algorithm) is used to compute the roots and poles of the rational approximation.

The partial fraction coefficients are then computed from (4.7). The new solution is found by solving the independent complex tridiagonal systems in

parallel and combining the results. In our algorithm we take advantage of

the feature of the decomposition alluded to in (4.9). In this fashion, for an

approximation of degree r we only need to solve [_J complex systems and

[_] - L_J real systems (corresponding to the possible real root).

Figure 3 shows the behavior of the error in one time-step for Pad6

and Chebyshev diagonal approximations of orders 1, 10, 4 and 14 as the

steplength At varies. Thus, it is the rational approximation analogue of

Figure 2. The case of diagonal Pad6 order 1 is of course identical to the

Crank-Nicolson scheme, so that figure indicates us the behavior of a standard method in comparison to the high-order methods we are proposing.

Judging from the errors, it would seem that one could discard the Pad6

schemes as inferior. It is known however that Chebyshev approximation

reaches maximum error at 0, exactly where Pad6 does best. An example of

the better behavior of Pad6 for small steps can be seen in Figure 4, where

At = 0.0005 to 0.01. We note that there exist techniques to a,'oid this error

behavior of Chebyshev approximation [23].

To test the performance of rational approximation method we used the 1-

dimensional problem but with u(0, z j) taken to be the jth componen t sin

of the eigenvector of A corresponding to the eigenvalue of smallest modulus.

The dimension n was chosen to be 98. Diagonal approximation was used

throughout (m = r). When m = 1 real arithmetic was used. The objective was to integrate from 0 to T C [1, 1 + At] so that the maximtun error

(llerrorl!oo)at T is less than e = 10 -9. The optimal Atop is the maximum

At which achieves error tolerance E at T. This is difficult to compute exactly

and we determined numerically an approximation Atov so that an underestimation will add only a minimal amount of iterations. We only show the

results for the Pad6 case with degrees 1, 2, 4, 6, 8. From Section 4.5 there

are 1, 1, 2, 3, and 4 tridiagonal systems to be solved per step (using LU).

Our results are summarized in Table 2 and Figure 5.

17
2 r 1 i E r r _ J J

-2

-4

-8

-10

p] -----------

/

CIO ...............

C14 ":5__" .................................................................................

"120 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.0_ 0.09 O.1

Timestepsize

FIG. 3. Behavior of rational approximation: Pro, Crn are the curves ]or rn th order Padd

and Chebyshev respectively as At varies.

The experiment demonstrates the following Crucial facts:

1. Crank-Nico!son on 8 CEs achieves a speedup of 1.43 over its 1 CE

run.

2. The 8th order scheme achieves a speedup of 3.5 (for 4 or more CEs)

over its 1 CE run.

3. The 8th order scheme achieves a speedup of 167 over Crank-Nicolson,

ltem (1) shows the difficulty of the low order scheme to profit from parallel

processing. Item (2) shows the considerably better behavior of the higher

order schemes due to the use of partial fractions. Item (3) shows the excellent

behavior and potential of high degree methods compared with standard low

order schemes.

6. Concluding remarks. We have proposed three parallel techniques

for solving parabolic equations. The first method based on Krylov subspaees

is the easiest to implement. It has the advantage of not requiring any solution

of linear systems. On the other hand it is basically an explicit method and

may be inefficient for very stiff problems. The other two methods rely on a

18
a_

-1

-2

-3

-4

-5

-6

-7

-8

-9

-1o

o

i i i i 1

C4

i i _t_._____ i t i i i ++_

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Timestep size

FIG. 4. Behavior of fourth (rn = 4) order Pad_ (P_) and Chebyshev (C4) rational approzimations as At varies.

rational approximation to the exponential. The basic idea of their parallel

implementations is to resort to partial fraction expansions. This transforms

the basic problem of solving a linear system with a product of matrices into

that of solving independent linear systems.

We would like to conclude with two comments, placing ourselves in the

more general framework of the parallel solution of systems of Ordinary Differential Equations. First, it is becoming apparent that explicit methods will

regain interest with parallel processing. These methods are particularly appealing for three-dimensional problems, especially in conjunction with highly

accurate schemes. Second, high order integration methods seem to be important in ODE methods, in order to achieve parallelism. In one-dimensional

problems they are mandatory since each step requires solving a tridlagonal

system, with little room for parallelization. For two or three dimensional

problems, the use of the techniques based on partial fraction expansions described in this paper, allow us to bypass the need to parallelize the sparse

linear system solvers which are difficult to optimize on supereomputers.

19
TABLE 2

Performance of different degree Pad_ schemes on the Alliant FX/8.

At

0.4910E-03

0.1950E-01

0.1600E+00

0.4000E+00

0.5000E+00

m Steps Time

1 2037 0.1838E+01

2 52 0.1740E+00

4 7 0.2590E-01

6 3 ' 0.1320E-01

8 2 0.1110E-01

Acknowledgement. The research of the firstauthor was supported by

the National Science Foundation under Grants No. US NSF-MIP-8410110,

US NSF DCR85-09970, US NSF CCR-8717942 and by AT&T Grant AT&T

AFFL67Sarneh. The research of the second author was supported by NASA

under USRA Grant No. NCC 2-387.

REFERENCES

[I] A. J. CARPENTER, A. RUTTAN, AND R. S. VARGA. Extended numerical computations on the I/9 conjecture in rational approximation theory. In P. R. GravesMorris, E. B. Saff, and R. S. Varga, editors, Rational Approzimation and Interpolation, pages 383-411, Sprlnger-Verlag, Berlin, i984.

[2] J. C. CAVENDISH, W. E. CULIIAM, AND R. S. VARGA. A comparison of CrankNicolson and Chebyshev rational methods for numerically solving linear parabolic

equations. J. Comput. Phys., 10:354 368, 1972.

[3] W. J. CODY, G. MEINARDUS, AND R. S. VARGA. Chebyshev rational approximations to e -_ in (0, +or) and appilcat_ons to heat-conductlon problems. J. Approt.

Theory, 2(1):50-65, March 1969.

[4] G. FAIRWEATHER. A note on the efficient implementation of certain Padd methods

for linear parabolic problems. BIT, 18:106-I09, 1978.

[5] R. FREUND. On conjugate gradient type methods and polynomial preconditioners for

a class of complez non-Hermitian matrices. Technical Report 88-44, RIACS,

NASA Ames Research Center, 1989.

[6] E. GALLOPOULOS. A partial fraction decomposition approach to improved efficiency

of some parabolic solvers. Technical Report 874, Center for Supercomputing

Research and Development, May 1989.

[7] E. GALLOPOULOS AND Y. SAAD. Parallel block cyclic reduction algorithm for the fast

solution of elliptic equations. Technical Report 659, Center for Supercomputing

Research and Development, April 1987.

[8] C. W. G_AR. The potential for parallelism in ordinary differential equations. Technical Report R-86-1246, Department of Computer Science, University of Illinois

at Urbana Champaign, Feb. 1986.

[9] P. HENRICI. Applied and Computational Complez Analysis. Volume 1, Wiley, 1974.

[10] D. J. KUCK, E. S. DAVIDSON, D. L. LAWRIE, AND A. H SAMEH. Paralle[ supercomputing today and the Cedar approach. Science, 231:967-974, February 1986.

[11] tI. T. KUNG. New algorithms and lower bounds for the parallel evaluation of certain

rational expressions and recurrences. J. ASsoc. Comput. Mach., 23(2):252-26i,

20
0.5

-0.5

-1.5

-2

m=l

m=2

iiiiiiiii ...........................................

m_

I I I I t I l

"2'50 i 2 3 4 5 6 7 8 9

Pf0c_g0rs

FIG. 5. Performance of rational Pad( approximation of degree m on the Alliant FX/8

varying the number of processors.

April 1976.

[12] J. D. LAWSON AND D. A. SWAYNE. High-order near best uniform approximations

to the solution of heat conduction problems. In Proc. IFIP Congress 80 - In/or.

mation Processing 80, pages 741-746, North Holland, New York, 1980.

[13] C. MOLER AND C. VAN LOAN. Nineteen dubious ways to compute the exponential

of a matrix. SIAM Rev., 20(4):801-836, October 1978.

[14] S. P. NORSETT. Restricted Pad_ approximations to the exponential function. SIAM

J. Numer. Anal., 15(5):1008-1029, Oct, 1978.

[15] B. NOUR-OM/D. Applications of the Lanczos algorithm. Comp. Phy. Comm., 53,

1989.

[16] B. N. PARLE'IT. The Symmetric Eigenvalue Problem. Prentice Hall, Englewood

Cliffs, 1980.

[17] M. F. REUSCH, L. RAT'ZAN, N. POMPHREY, AND W. PARK. Diagonal Pad_ approximations for initial value problems. SIAM J. Sci. Statist. Comput., 9(5):829-838,

September 1988.

21
[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

M. J. SCHAEFER. A polynomial based iterative method/or linear parabolic equations.

Technical Report 661, Center for Supercomputlng Research and Development,

May 1987.

D. A. SWAYNE. Matrix operations with rational functions. In Proc. 7th Manitoba

Con/. Numerical Mathematics, pages 581-589, Utilitas Mathematica, Winnipeg,

Manitoba, 1977.

R. A. SWEET. A parallel and vector cyclic reduction algorithm. S[AM J. Sci. Statist.

Comput., 9(4):761-765, July 1988.

1[. TAL-EZER. Spectral methods in time for parabolic problems. SIAM J. Nurner.

Anal., 26(1):1-]1, Feb. ]989.

R. S. VARGA. On higher order stable implicit methods for solving parabolic partial

differential equations. J. Math. Phys., 40:220-231, 1961.

R. S. VARGA. Some results in approximation theory with applications to numerical

analysis, in B. E. 1[ubbard, editor, Numerical Solution o/Partial Differential

Equations - [I, pages 623-649, Academic Press, New York, 197I.

J. H. WILKINSON. The Algebraic Eigenvalue Problem. Oxford University Press,

1965.

R. E. WYATT. The recursive residue generation method. Adv. Chem. Phy. Chap. 5,

1988.

V. ZAKIAN. Properties of IMN and JMN approximant_and applications to numerical inversion of Laplace transforms and initial value problems. J. Math. Anat.

Applic., 50:191-222, i975.

22