STABLE PARALLEL ALGORITHMS FOR TWO-POINT BOUNDARY VALUE PROBLEMS^{∗}

STEPHEN J. WRIGHT^{†}

SIAM J. SCI. STATIST. COMPUT. \begin{matrix}c \\ ◯ \end{matrix} 1991 Society for Industrial and Applied Mathematics Vol. 0, No. 0, pp. xxx–xxx, Month 1992 012

Abstract. Some of the most widely used algorithms for two-point boundary value ODEs, namely finite difference and collocation methods and standard multiple shooting, proceed by setting up and solving a structured system of linear equations. It is well known that the linear system can be set up efficiently in parallel; we show here that a structured orthogonal factorization technique can be used to solve this system, and hence the overall problem, in an efficient, parallel, and stable way.

Key words. parallel algorithms, two-point boundary value problems, error analysis and stability

AMS(MOS) subject classifications. 65F05, 65G05, 65L10, 65L20, 65W05

1. Introduction. Many numerical algorithms for solving the linear two-point boundary value problem

y′ = M(t)y + q(t), t ∈ [a, b], y ∈ R^{n}, (1) B_{a}y(a) + B_{b}y(b) = d, (2)

have been proposed and studied over the last 30 years. Many of these methods require a structured linear algebraic system (for example, a block-tridiagonal or staircase system) to be solved as a “core” operation, and considerable effort has been devoted to minimizing the amount of computer time and storage required during the factorization of the coefficient matrix of this system. Efficient factorization schemes, based on structured Gaussian elimination, have been implemented and are widely available (see §2, and Varah [19], Diaz, Fairweather and Keast [7], Lentini and Pereyra [13], and Keller [9].)

During the last 10 years, the question of stability of algorithms for (1),(2) has received a great deal of attention (see, for example, Mattheij [15].) It has been recognized that in a well-conditioned problem (that is, one whose solution is not too sensitive to perturbations in M, q or the boundary conditions (2)), the fundamental solution space generally contains both exponentially increasing and exponentially decreasing modes. The stability of a numerical method for (1),(2) depends on its ability to at some point perform a “decoupling” of these modes. For the standard multiple shooting and finite difference algorithms, this decoupling is performed implicitly, during the factorization of the structured linear system, through the use of a pivoting strategy that prevents element growth in the factors. Unfortunately, parallel and vectorizable algorithms that have been proposed for solving the linear system invariably place some restriction on the choice of pivots. This can lead to undesirable element growth in the factors, and such methods are in fact similar to compactification algorithms for (1),(2) which are known to be unstable. In this paper, we use instead a structured orthogonal factorization technique which is stable, has an identical serial complexity to the best-known algorithms, and can be efficiently implemented on a wide variety of parallel architectures. A variant of the algorithm vectorizes efficiently,

∗ Received by the editors September xx, 1990; accepted for publication (in revised form) January xx, 1991. This research was supported by the Applied Mathematical Sciences subprogram of the Office of Energy research, United States Department of Energy, contract W-31-109-Eng-38. A grant of computer time at the North Carolina Supercomputing Center is gratefully acknowledged.

† Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439.

in much the same way as cyclic reduction for block-tridiagonal linear systems. No separability of the boundary conditions is needed. Although we focus on matrices arising from multiple shooting and one-step differencing schemes, the technique can be applied equally well to the more general staircase matrix structures which arise in other numerical schemes, such as collocation.

We assume throughout that n is too small to allow efficient parallel or vector implementation of order-n matrix and vector operations. Instead, parallelism is sought by partitioning the domain [a, b] of the independent variable.

The remainder of the paper is organized as follows: in §2 we review the multiple shooting and finite difference algorithms, the structured linear systems that result from them, and the efficient Gaussian elimination techniques used to solve them on serial computers. In §3 we briefly review previous work on parallel and vectorizable algorithms. The new algorithm is presented in §4, together with a stability result and analysis of serial and parallel complexity. In §5 we describe parallel implementation of a standard technique for estimating the condition number of the coefficient matrix, which can be used for purposes of a posteriori error analysis. Finally in §6 we describe the results of implementation of the scheme on shared-memory and vector architectures.

2. Serial algorithms. The “standard” multiple shooting technique for (1),(2) proceeds by defining a mesh

a = t_{1} < t_{2} < … < t_{k+1} = b, (3)

and finding a fundamental and particular solution on each interval of the mesh:

Y′_{i} = M(t)Y_{i}, Y_{i}(t) ∈ R^{n×n}, t ∈ [t_{i}, t_{i+1}], Y_{i}(t_{i}) = I,

y′_{pi} = M(t)y_{pi} + q(t), y_{pi}(t) ∈ R^{n}, t ∈ [t_{i}, t_{i+1}], y_{pi}(t_{i}) = 0.

Then, we try to find s_{i} ∈ R^{n}, i = 1, …, k + 1, such that

y(t) = Y_{i}(t)s_{i} + y_{pi}(t), t ∈ [t_{i}, t_{i+1}], i = 1, …, k.

(Note in particular that s_{i} = y(t_{i}), i = 1, …, k + 1.) By applying the boundary conditions (2), and continuity at the meshpoints, we obtain the following linear system in s_{1}, …, s_{k+1}:

[ \begin{matrix}B_{a} & & & & B_{b} \\ Y_{1}(t_{2}) & −I \\ & Y_{2}(t_{3}) & −I \\ & & ... & ... \\ & & & Y_{k}(t_{k+1}) & −I \end{matrix} ] [ s_{1} ; s_{2} ; s_{3} ; ⋮ ; s_{k+1} ] = [ d ; −y_{p1}(t_{2}) ; −y_{p2}(t_{3}) ; ⋮ ; −y_{pk}(t_{k+1}) ] . (4)

In the (quite common) case in which the boundary conditions are separated, this system can be rearranged into a block-banded form. If p is the number of initial conditions and q is the number of final conditions (p + q = n), we can assume without loss of generality that the boundary data (2) can be partitioned as follows:

B_{a} = [ B̄_{a} ; 0 ] , B_{b} = [ 0 ; B̄_{b} ] , d = [ \begin{matrix}d_{a} \\ d_{b} \end{matrix} ] .

If the rows p + 1, …, n of (4) are moved to the bottom, we obtain the form

[ B̄_{a} ; Y_{1}(t_{2}) −I ; Y_{2}(t_{3}) −I ; ...  ... ; Y_{k}(t_{k+1}) −I ; B̄_{b} ] [ s_{1} ; s_{2} ; s_{3} ; ⋮ ; s_{k} ; s_{k+1} ] = [ d_{a} ; −y_{p1}(t_{2}) ; −y_{p2}(t_{3}) ; ⋮ ; −y_{pk}(t_{k+1}) ; d_{b} ] . (5)

One-step finite difference methods proceed by again choosing a mesh of the form (3) and seeking s_{1}, …, s_{k+1} to approximate y(t_{1}), …, y(t_{k+1}). On the interval [t_{i}, t_{i+1}], (1) is approximated by the linear relationship

\frac{s_{i+1} − s_{i}}{h_{i}} = Ā_{i}s_{i} + C̄_{i}s_{i+1} + f̄_{i}, (6)

where h_{i} = t_{i+1} −t_{i}. Two schemes with a local truncation error of O(h_{i}^{2}) are the “box scheme” (Ā_{i} = C̄_{i} = \frac{1}{2}M(t_{i+ \frac{1}{2}} ), f̄_{i} = q(t_{i+ \frac{1}{2}} )) and the trapezoidal rule (Ā_{i} = \frac{1}{2}M(t_{i}), C̄_{i} = \frac{1}{2}M(t_{i+1}), f̄_{i} = \frac{1}{2}(q(t_{i}) + q(t_{i+1})).) By including (2), we obtain a linear system of the form

[ \begin{matrix}B_{a} & & & & B_{b} \\ A_{1} & C_{1} \\ & A_{2} & C_{2} \\ & & ... & ... \\ & & & A_{k} & C_{k} \end{matrix} ] [ s_{1} ; s_{2} ; s_{3} ; ⋮ ; s_{k+1} ] = [ d ; f_{1} ; f_{2} ; ⋮ ; f_{k} ] , (7)

where A_{i} = −I − h_{i}Ā_{i}, C_{i} = I − h_{i}C̄_{i}, and f_{i} = h_{i}f̄_{i}. If the boundary conditions are separated, we obtain

[ B̄_{a} ; A_{1} C_{1} ; A_{2} C_{2} ; ... ... ; A_{k} C_{k} ; B̄_{b} ] [ s_{1} ; s_{2} ; s_{3} ; ⋮ ; s_{k} ; s_{k+1} ] = [ d_{a} ; f_{1} ; f_{2} ; ⋮ ; f_{k} ; d_{b} ] . (8)

The accuracy of finite difference schemes is often enhanced by the use of deferred correction techniques (see, for example, Pereyra [18]).

As has been observed, the two algorithms are closely related, in the sense that for a reasonable choice of the approximation (6), −C_{i}^{−1}A_{i} should be close to Y_{i}(t_{i+1}). Hence, the conditioning of the matrices in (4) and (7) is quite similar when the h_{i} are small. If we quantify the conditioning of the problem (1),(2) by a bound κ on its Green’s function (see Ascher, Mattheij, and Russell [3, §3.2]), it has been shown by Osborne [16] (and also by Mattheij [14] and Lentini, Osborne, and Russell [12]) that the inverse of A_{S} from (4) satisfies the bound

∥A_{S}^{−1}∥_{∞} ≤ kκ.

Hence if we define γ by

γ = 1 + max_{i=1,...,k} ∥Y_{i}(t_{i+1})∥_{∞}

and assume that B_{a} and B_{b} are scaled such that ∥[B_{a}, B_{b}]∥_{∞} = 1, then from (4),(5),

cond_{∞}(A_{S}) ≤ γkκ. (9)

For the finite-difference coefficient matrix A_{D} in (7), a similar analysis applies if we note that γ = 2 + O(h̄) for small h̄ = max_{i=1,...,k} h_{i}. We use this fact, together with A_{D} = A_{S} + O(h̄), to derive the bound

cond_{∞}(A_{D}) ≤ (2 + O(h̄))kκ. (10)

The slightly different form of the bounds (9) and (10) is motivated by the fact that k is typically larger in a finite difference algorithm than in a multiple shooting algorithm. Bounds for the coefficient matrices in (5) and (8) are, of course, identical to their non-separated counterparts.

As stated earlier, Gaussian elimination algorithms with various forms of pivoting have been previously proposed for solving (4),(5),(7),(8). The practical stability of such algorithms for general matrices is well known, but is also well known that the worst-case behavior can be very bad, as a result of possible growth of elements in the factors which is exponential in the dimension (k + 1)n of the system. It has been shown (see, for example, Mattheij [15]) that the presence of an exponential dichotomy in (1),(2) ensures that this worst-case behavior does not arise when elimination algorithms are applied to the matrices in (4),(5),(7),(8). In the partially separated cases (5),(8), similar results can be proved, without referring to the dichotomy at all, by using the results of Bohte [4]. Bohte shows that for banded linear systems, the bound on element growth in partial pivoting algorithms is exponential only in the bandwidth.

In the simplest case of Gaussian elimination with row partial pivoting (with coefficient matrices and right-hand sides denoted by A_{SP} and f_{SP} in (5) and A_{DP} and f_{DP} in (8)), possible fill-in of kpn elements can occur in the upper triangular factor. Element growth is, however, bounded by 2^{2n−1}, and so we obtain the following theorem:

Theorem 2.1. Let s = (s_{1}^{T} , …, s_{k+1}^{T})^{T} denote the true solution of the system (8), and suppose that ŝ is the approximate solution obtained by Gaussian elimination with row partial pivoting. Then provided that

(i) h̄ is chosen small enough that cond_{∞}(A_{DP}) ≤ 4kκ, (ii) 8c_{1}κu ≤ 1, where c_{1} = (1.12)2^{2n}n^{3}k(k + 1)(k + 8) = O(k^{3}n^{3}2^{2n}) and u is the

unit roundoff error, the following relative error bound applies:

\frac{∥ŝ − s∥_{∞}}{∥s∥_{∞}} ≤ 16c_{1}κu. (11)

Suppose that the fundamental and particular solutions that are used to construct A_{SP} and f_{SP} in (5) are calculated to a tolerance of τ, that is,

(12)  ∥A_{SP} − Ā_{SP} ∥_{∞} ≤ τ∥A_{SP} ∥_{∞}
∥f_{SP} − f̄_{SP} ∥_{∞} ≤ τ∥f_{SP} ∥_{∞}.

Let s̃ be the solution obtained by using Gaussian elimination with row partial pivoting. Then provided that

(iii) τkκγ ≤ 1/2, (iv) 6c_{1}uκγ ≤ 1, we have

(13) \frac{max_{i=1,...,k+1} ∥s̃_{i} − y(t_{i})∥_{∞}}{max_{i=1,...,k+1} ∥y(t_{i})∥_{∞}} ≤ 4γκ[kτ + 3c_{1}u + c_{2}uτ],

where c_{2} = 12c_{1}kκγ.

Proof. For the first part of the theorem, the proof follows from Theorem 2.7.2 in Golub and Van Loan [8], and §4 of Bohte [4]. The latter assumes that

u ≤ .009  and  (k + 1)nu ≤ 0.1,

which are clearly implied by assumption (ii), since κ ≥ 1, k ≥ 1 and n ≥ 1. Bohte then shows that ŝ solves the system

(A_{DP} + E)ŝ = f_{DP} ,

where E satisfies the bound

∥E∥_{∞} ≤ Λu(0.56)(4n − 1)n(k + 1)(kn + 7n + 3),

and Λ is a bound on the maximal element which arises during the elimination. By simplifying this expression, and noting from the discussion above that

Λ ≤ 2^{2n−1} max_{i,j} |(A_{DP})_{ij}| ≤ 2^{2n−1}∥A_{DP} ∥_{∞},

we have that

∥E∥_{∞} ≤ (c_{1}/k)u∥A_{DP} ∥_{∞}.

The result now follows by setting r = 1/2 in [8, Theorem 2.7.2].

For the second part, note first that

∥Ā_{SP} ∥_{∞} ≤ (1 + τ)∥A_{SP} ∥_{∞}

and

∥Ā_{SP}^{−1} ∥_{∞} ≤ ∥[I + A_{SP}^{−1} (Ā_{SP} − A_{SP})]^{−1}∥_{∞}∥A_{SP}^{−1} ∥_{∞}.

Provided

∥A_{SP}^{−1} (Ā_{SP} − A_{SP})∥_{∞} ≤ 1, (14)

we have from [8, Lemma 2.3.3] that

∥Ā_{SP}^{−1} ∥_{∞} ≤ \frac{∥A_{SP}^{−1} ∥_{∞}}{1 − ∥A_{SP}^{−1} ∥_{∞}τ} ≤ \frac{kκ}{1 − τkκ} ≤ 2kκ.

The inequality (14) is implied by (iii). Since k, γ and κ are all at least 1, assumption (iii) implies that τ ≤ 1/2, and so

cond_{∞}(Ā_{SP}) ≤ 2(1 + τ)kκγ ≤ 3kκγ. (15)

Now let s̄ be the exact solution of Ā_{SP} s̄ = f̄_{SP}. Direct application of [8, Theorem 2.7.2] with r = 1/2 shows, again using (iii), that

(16) \frac{max_{i=1,...,k+1} ∥s̄_{i} − y(t_{i})∥_{∞}}{max_{i=1,...,k+1} ∥y(t_{i})∥_{∞}} ≤ 4cond_{∞}(A_{SP})τ ≤ 4kκγτ.

Application of Bohte’s results shows that the computed solution s̃ satisfies (Ā_{SP} + E)s̃ = f̄_{SP}, where

∥E∥_{∞} ≤ \frac{c_{1}}{k} u∥Ā_{SP} ∥_{∞}.

Assumption (iv) and (15) imply that

\frac{c_{1}}{k} ucond_{∞}(Ā_{SP}) ≤ 1/2,

and so

\frac{∥s̃ − s̄∥_{∞}}{∥s̄∥_{∞}} ≤ \frac{4c_{1}u}{k}cond_{∞}(Ā_{SP}) ≤ 12c_{1}uκγ. (17)

The result follows by combining (16) and (17) in an elementary way.

A more economical scheme for solving (5) and (8), described by Varah [19] and implemented by Diaz, Fairweather, and Keast [7], is the method of alternate row and column elimination. For the first p stages of the process, we use column pivoting and elimination (involving columns 1 through n); this produces no fill-in. For stages p + 1 through n we use row pivoting and elimination for the same reason. Column and row elimination alternate in this way until a factorization of the form

A_{DP} = P LBUΠ (18)

is produced, where P and Π are permutation matrices, L and U are lower and upper triangular matrices of multipliers, and B has the form

B = [ \begin{matrix}L \\ X & R & X & X \\ X & & L \\ & & X & R & X & X \\ & & X & & L \\ & & & & ... & ... & ... \\ & & & & & & X & 0 & L \\ & & & & & & & & X & R \end{matrix} ] .

(L denotes a lower triangular p×p block, R denotes an upper triangular (n−p)×(n−p) block, and X denotes a dense block.) It is easy to show that element growth in the B factor is bounded by 2^{n−1} and, hence, that the scheme is very stable. We can prove the following:

Theorem 2.2. Let s = (s_{1}^{T} , …, s_{k+1}^{T})^{T} denote the true solution of the system (8), and let ŝ be the approximate solution calculated by alternate row and column elimination. Suppose further that

(i) h̄ is chosen small enough that cond_{∞}(A_{DP}) ≤ 4kκ, and (ii) 8c_{4}κu ≤ 1,

where

c_{4} = k(k + 2)n(5 + 6n)(1 + n^{2}2^{n−1}) = O(k^{2}n^{4}2^{n}).

Then

\frac{∥ŝ − s∥_{∞}}{∥s∥_{∞}} ≤ 16c_{4}κu.

Let s̃ be obtained by multiple shooting with alternate row and column elimination, and suppose that (12) holds, with (iii) τkκγ ≤ 1/2, (iv) 6c_{4}uκγ ≤ 1.

Then

\frac{max_{i=1,...,k+1} ∥s̃_{i} − y(t_{i})∥_{∞}}{max_{i=1,...,k+1} ∥y(t_{i})∥_{∞}} ≤ γkκ[kτ + 3c_{4}u + c_{5}uτ],

where c_{5} = 12c_{4}kκγ.

Proof. The assumptions of Lemma A.1 are consequences of (ii) and (iv) above. The first part of the result follows directly from Lemma A.2 and [8, Theorem 2.7.2]. The second part is analogous to the second part of Theorem 2.1, with c_{4} replacing c_{1} and c_{5} replacing c_{2}.

3. Parallel elimination algorithms. Other parallel and vectorizable algorithms, based on Gaussian elimination, have recently been proposed for (4),(5),(7),(8). In general, they suffer either from poor stability properties or from limitations in the amount of parallelism which is possible.

Wright and Pereyra [21] describe variants of a block factorization algorithm applied to (7). In the most highly vectorized variant, a factorization of the form

[ Ã_{1} ; Ã_{2} ; ... ; Ã_{k} ; Z_{1} Z_{2} … Z_{k} Ã_{k+1} ] [ \begin{matrix}I & W_{1} \\ & I & W_{2} \\ & & ... & ... \\ & & & I & W_{k} \\ & & & & I \end{matrix} ] =

[ \begin{matrix}P_{1} \\ & P_{2} \\ & & . & . & . \\ & & & & & P_{k} \\ & & & & & & P_{k+1} \end{matrix} ] [ \begin{matrix}A_{1} & C_{1} \\ & A_{2} & C_{2} \\ & & ... & ... \\ & & & A_{k} & C_{k} \\ B_{a} & & & & B_{b} \end{matrix} ]

(where the P_{i} are permutation matrices) is produced. It is easy to show that for h̄ sufficiently small, this factorization exists. It is equivalent to compactification [3, page 153] which, because of its similarity to single shooting, is known to be potentially unstable. However, this instability can usually be recognized by the presence of large elements in the Z_{i} blocks. The strategy described in [21, §5] is to use this factorization where possible, and discard it in favor of a more stable method if the ∥Z_{i}∥ are too large. In many applications (such as the one described in [21]) the lack of stability is not a problem.

Paprzycki and Gladwell [17] describe a partitioned elimination algorithm in which (8) is torn into P submatrices, each of which has the same form as the original A_{DP}, and alternate row and column elimination is applied to each piece. This corresponds to partitioning the interval [a, b] and solving a number of sub-BVPs, each of which has p initial and q final conditions. Although the number of initial and final conditions is correct, this alone is not enough to ensure well-conditioning of the sub-BVPs. In a linear algebra sense, well-conditioning of the whole matrix A_{DP} does not guarantee well-conditioning of each of the P submatrices. Ascher and Mattheij [2] develop a “theoretical multiple shooting” framework in which they show how boundary values for the sub-BVPs should be chosen to ensure well-conditioning. Ascher and Chan [1] suggest how to implement this in a parallel environment.

Another possibility, which leads to near-perfect speedup on two processors (but cannot be generalized to a higher order of parallelism) is the “burn at both ends” or “twisted” factorization. Here, some form of pivoted Gaussian elimination is applied simultaneously from both ends of the matrix (either A_{DP} or A_{SP}). When the factorizations meet in the center, a small reduced system is formed and factored. This is analogous to the approach of Lentini [11].

Finally, we mention the approach of Ascher and Chan [1], who form the normal equations for (5) and (8), and factorize the resulting symmetric, positive definite, block-tridiagonal system using cyclic reduction. In exact arithmetic, this scheme produces a triangular factor which is identical to that given by the “cyclic reduction” variant of the algorithm to be described in the next section. The difference lies in the fact that by explicitly forming the normal equations, the condition number of a linear system is squared, an effect that is avoided when the algorithm from §4 is used.

4. Structured orthogonal factorization. We now describe the structured QR factorization algorithm, as applied to the system (7). It can of course be applied equally well to the systems (8),(4),(5), since it is indifferent to separability of the boundary conditions.

The first step is to partition the system into, say, P pieces of approximately equal size. We choose indices k_{1}, k_{2}, …, k_{P +1} to satisfy

0 = k_{1} < k_{2} < … < k_{P +1} = k, k_{j+1} ≥ k_{j} + 2, j = 1, …, P.

(P could for instance be the number of processors on the physical machine being used to solve the problem.) Partition j then consists of rows (k_{j} + 1)n + 1, …, (k_{j+1} + 1)n of (7). Each partition is now processed independently; in effect, the variables s_{i} for i ≠ k_{j} + 1, j = 0, …, P, are eliminated from the problem.

We describe this process in detail for the first partition, which, if augmented with its right-hand side, has the form

[ \begin{matrix}A_{1} & C_{1} & & & & f_{1} \\ & A_{2} & C_{2} & & & f_{2} \\ & & ... & ... & & ⋮ \\ & & & A_{k_{2}} & C_{k_{2}} & f_{k_{2}} \end{matrix} ] . (19)

We first find an orthogonal matrix Q_{1} ∈ R^{2n×2n} such that

Q_{1}^{T} [ \begin{matrix}C_{1} \\ A_{2} \end{matrix} ] = [ \begin{matrix}R_{1} \\ 0 \end{matrix} ] ,

where R_{1} ∈ R^{n×n} is upper triangular. If we form Q_{1} as a product of n Householder transformations, the information needed to reconstruct Q_{1} could be stored in the space formerly occupied by the “zeroed” elements of [C_{1}^{T} , A_{2}^{T}]^{T}, plus an additional n locations. We also need to apply Q_{1}^{T} to the other columns of the matrix (19), and the overall effect is

[ \begin{matrix}Q_{1}^{T} & 0 \\ 0 & I \end{matrix} ] [ \begin{matrix}A_{1} & C_{1} & & & & f_{1} \\ & A_{2} & C_{2} & & & f_{2} \\ & & ... & ... & & ⋮ \\ & & & A_{k_{2}} & C_{k_{2}} & f_{k_{2}} \end{matrix} ] =

[ \begin{matrix}G_{1} & R_{1} & E_{1} & & & & g_{1} \\ Ḡ_{2} & 0 & C̄_{2} & & & & f̄_{2} \\ & & A_{3} & C_{3} & & & f_{3} \\ & & & ... & ... & & ⋮ \\ & & & & A_{k_{2}} & C_{k_{2}} & f_{k_{2}} \end{matrix} ] . (20)

The next step is to find an orthogonal Q_{2} ∈ R^{2n×2n} such that

Q_{2}^{T} [ C̄_{2} ; A_{3} ] = [ \begin{matrix}R_{2} \\ 0 \end{matrix} ] ,

and to apply Q_{2}^{T} to rows n + 1, …, 3n of the reduced system (20). This process is repeated a total of k_{2} − 1 times, until finally we obtain a system equivalent to (19) which has the form

[ \begin{matrix}G_{1} & R_{1} & E_{1} & & & & g_{1} \\ G_{2} & & R_{2} & E_{2} & & & g_{2} \\ ⋮ & & & ... & ... & & ⋮ \\ G_{k_{2}−1} & & & & R_{k_{2}−1} & E_{k_{2}−1} & g_{k_{2}−1} \\ Ã_{1} & & & & & C̃_{1} & f̃_{1} \end{matrix} ] . (21)

Formally, the reduction process for partition j can be specified as follows:

C̄_{k_{j}+1} = C_{k_{j}+1}, Ḡ_{k_{j}+1} = A_{k_{j}+1}, f̄_{k_{j}+1} = f_{k_{j}+1}.

for i = k_{j} + 1, …, k_{j+1} − 1 Find orthogonal Q_{i} such that

Q_{i}^{T} [ C̄_{i} ; A_{i+1} ] = [ \begin{matrix}R_{i} \\ 0 \end{matrix} ] , (22)

where R_{i} ∈ R^{n×n} is upper triangular; Set

[ E_{i} ; C̄_{i+1} ] ← Q_{i}^{T} [ \begin{matrix}0 \\ C_{i+1} \end{matrix} ] , [ G_{i} ; Ḡ_{i+1} ] ← Q_{i}^{T} [ Ḡ_{i} ; 0 ] ,

[ g_{i} ; f̄_{i+1} ] ← Q_{i}^{T} [ f̄_{i} ; f_{i+1} ] .

end (for)

Set Ã_{j} = Ḡ_{k_{j+1}}, C̃_{j} = C̄_{k_{j+1}}, f̃_{j} = f̄_{k_{j+1}}.

Clearly, if we knew the values of s_{k_{j}+1} and s_{k_{j+1}+1}, the values of s_{k_{j}+2}, …, s_{k_{j+1}} could be recovered by the simple back-substitution:

for i = k_{j+1} − 1, k_{j+1} − 2, …, k_{j} + 1

Solve R_{i}s_{i+1} = g_{i} − G_{i}s_{k_{j}} − E_{i}s_{i+2} (23)

In fact, these “separator” variables can be found by forming and solving a reduced system by taking the last n rows of each (processed) partition, and the boundary conditions. This reduced system has the form

[ \begin{matrix}B_{a} & & & & B_{b} \\ Ã_{1} & C̃_{1} \\ & Ã_{2} & C̃_{2} \\ & & ... & ... \\ & & & Ã_{P} & C̃_{P} \end{matrix} ] [ s_{k_{1}+1} ; s_{k_{2}+1} ; s_{k_{3}+1} ; ⋮ ; s_{k_{P +1}+1} ] = [ d ; f̃_{1} ; f̃_{2} ; ⋮ ; f̃_{P} ] . (24)

The form of (24) is obviously identical to the form of the original system (7), so this immediately suggests that it may be possible to apply the whole process recursively, that is, (24) could be partitioned into P_{2} ≤ P/2 pieces, and the algorithm described above could be applied to obtain a smaller reduced system. This process could be repeated for, say, L levels, so at the innermost level a (P_{L} + 1)n-dimensional system would remain.

In a parallel implementation of the algorithm, the number of levels L to be used and the number of processors to be used at each level depend on the number of available processors on the physical machine, and, on a distributed-memory machine, on the cost of interprocessor communication. The “extreme” cases are as follows:

• A one-level version (P = 1, k_{2} = k), in which (7) is reduced to the 2n × 2n system

[ \begin{matrix}B_{a} & B_{b} \\ Ã_{1} & C̃_{1} \end{matrix} ] [ \begin{matrix}s_{1} \\ s_{k+1} \end{matrix} ] = [ d ; f̃_{1} ] , (25)

which is then solved by QR factorization. This is the “serial” version of the algorithm.

• A two-level version, in which P processors are used to do the reduction and back-substitution (23), and the system (24) is solved by the one-level algorithm just described. It is easy to see that, assuming that each partition contains about k/P rows of blocks, the time required for the reduction phase is proportional to (k/P − 1)n^{3}. The time required to solve (24) on a single processor is proportional to (P + 1)n^{3}. The latter operation will not be a bottleneck for the whole process, provided that k ≫ P^{2}.

• A “cyclic reduction” version, in which P is equal to the number of available processors (assuming that P ≤ k/2), and then P_{2} = P/2, P_{3} = P/4, …, P_{L} = 1. Here clearly L = log_{2} P + 1. The total computation time required will be proportional to (k/P − 1 + log_{2} P)n^{3}, while the communication time will be proportional to n^{2} log_{2} P. With regard to complexity, the algorithm is optimal: when k = P log_{2} P, the execution time is O(log_{2} P) on P processors, while the serial time is O(P log_{2} P).

The cyclic reduction variant is also appropriate for implementation on a single vector processor. In this environment, we could choose P = k/2, P_{2} = k/4, P_{3} = k/8, …, and so L = log_{2} k + 1. At level l, the reduction and back-substitution processes can be coded so that most of the arithmetic involves vectors of length k/2^{l}. At low levels, a reversal of the loop nesting can be used to ensure that vector lengths do not fall below n.

It is not difficult to see that the factorization and solution phases can be separated, provided that the Householder vectors are stored (in the locations vacated by the zeroed elements). This feature is useful when quasi-linearization is used to solve a firstorder nonlinear BVP. The same coefficient matrix (and its factorization) can be used for a number of consecutive iterations to produce “chord method” approximations to Newton steps.

It is important to note that the scheme proposed above is simply standard Householder QR factorization applied to an initial row- and column-permuted form of the original matrix. Thus, we can apply the standard stability analysis for (k + 1)ndimensional matrices from Lawson and Hanson [10] to obtain error bounds for the computed solutions. We have refined these bounds to take into account the structure of the matrix. In addition, to allay any possible concerns about instability at the level of the O(u^{2}) terms, we have removed these terms using the style of analysis in Wilkinson [20]. The relevant results are stated and proved in Appendix B. Here, we summarize the analysis in the following Lemma and Theorem:

Lemma 4.1. Let s = (s_{1}^{T} , s_{2}^{T} , …, s_{k+1}^{T})^{T} denote the true solution of the system (7), and suppose that ŝ is the approximate solution obtained by using the structured QR factorization and back-substitution. Suppose that the orthogonal matrices used in the factorization are all constructed from Householder reflectors, and that n, k and u satisfy

(12n + 51)n(k + 1)u ≤ 0.1. (26)

Then ŝ is the exact solution of a perturbed system

(A_{D} + δA_{D})ŝ = f_{D} + δf_{D}, (27)

where

(28)  ∥δA_{D}∥_{F} ≤ (1.106)(12n + 51)(k + 2)nu∥A_{D}∥_{F} ,
(29)  ∥δf_{D}∥_{2} ≤ (1.106)(12n + 51)(k + 1)nu∥f_{D}∥_{2}.

Proof. See Appendix B.

Theorem 4.2. Suppose that the conditions of Lemma 4.1 are satisfied. Then provided that

(i) h̄ is chosen small enough that cond_{∞}(A_{DP}) ≤ 4kκ,

(ii) 8c_{6}κu ≤ 1, where c_{6} = (1.106)k(k + 2)^{2}n^{2}(12n + 51),

the following relative error bound applies:

\frac{∥ŝ − s∥_{∞}}{∥s∥_{∞}} ≤ 16c_{6}κu.

Suppose that the fundamental and particular solutions that are used to construct A_{S} and f_{S} in (4) are calculated to a tolerance of τ, that is,

∥A_{S} − Ā_{S}∥_{∞} ≤ τ∥A_{S}∥_{∞}
∥f_{S} − f̄_{S}∥_{∞} ≤ τ∥f_{S}∥_{∞}.

Let s̃ be the solution obtained by using structured QR algorithm, where (26) holds. Then provided that

(iii) τkκγ ≤ 1/2, (iv) 6c_{6}uκγ ≤ 1, we have

\frac{max_{i=1,...,k+1} ∥s̃_{i} − y(t_{i})∥_{∞}}{max_{i=1,...,k+1} ∥y(t_{i})∥_{∞}} ≤ 4γκ[kτ + 3c_{6}u + c_{7}uτ],

where c_{7} = 12c_{6}kκγ.

Proof. The proof is identical to that of Theorem 2.1, once we convert the Frobenius norm bounds (28) and (29) into ∞-norm bounds. This is done by noting that for an N × N matrix A,

\frac{1}{N^{1/2}} ∥A∥_{∞} ≤ ∥A∥_{F} ≤ N^{1/2}∥A∥_{∞}.

Hence

∥δA_{D}∥_{∞} ≤ (1.106)(k + 2)^{2}n^{2}(12n + 51)u∥A_{D}∥_{∞}.

5. Condition estimation. For purposes of assessing the reliability of the computed solution, it is useful to have an estimate of the conditioning of the discrete system. Such an estimate can be obtained, simultaneously with the factorization and solution process, by adapting the procedure described in Cline et al. [6] to our situation.

We aim to compute an estimate of the quantity

κ̂ = ∥A∥_{∞}∥R^{−1}∥_{∞},

where A is one of the coefficient matrices from (4),(7),(5),(8), and R is the upper triangular factor produced by the procedure just described. It is easy to show that

\frac{1}{\sqrt{(k + 1)n}}κ̂ ≤ cond_{∞}(A) ≤ \sqrt{(k + 1)n}κ̂.

∥A∥_{∞} can of course be calculated directly, and so computation of the estimate of ∥R^{−1}∥_{∞} is the major part the task of finding κ̂. Following [6], we do this by first finding vectors z and v such that

R^{T} v = z,

where the components of z are all ±1 and are chosen by a heuristic that attempts to maximize ∥v∥_{∞}. This is done during the factorization of A; as reduction of each partition into a single row of blocks is performed, the heuristic can be applied to find the components of z and v corresponding to the rows and columns of the original system that are eliminated. Next, the solution of

Rw = v

is found concurrently with the solution of the original linear system Rs = Q^{T} f. We then use the estimate

∥R^{−1}∥_{∞} ≈ \frac{∥w∥_{∞}}{∥v∥_{∞}} .

The operation count for calculation of κ̂ is approximately four times that of doing a single backsolve with R, and the parallel complexity is the same as that of the factorization and solution. Comparisons with the LINPACK condition number estimator of cond_{1}(A) (which is based on an LU factorization of A but uses similar heuristics) show that the two estimates are usually within a factor of 3 of each other.

Table 1 Operation counts and storage requirements for four algorithms, assuming separated end con-

ditions. k = number of meshpoints, n = dimension of y, p = number of left-hand end conditions, R = number of right-hand sides

Algorithm  Operation count  Storage
LU (row pivoting) k[\frac{5}{3}n^{3} + 3pn^{2} + R(4n^{2} + 2pn)] kn(2n + p)
Structured QR k[\frac{46}{3}n^{3} + (15R + 30)n^{2}] 4kn^{2}
Normal equations k[\frac{38}{3}n^{3} + 12Rn^{2}]  4kn^{2}
DECOMP/SOLVE  k[\frac{2}{3}n^{3} + (4R + 5p)n^{2} − 2np^{2}]  2kn^{2}

Table 2 Operation counts and storage requirements for four algorithms, assuming non-separated end

conditions. k = number of meshpoints, n = dimension of y, R = number of right-hand sides

Algorithm  Operation count  Storage
LU (row pivoting) k[\frac{23}{3}n^{3} + 8Rn^{2}]  4kn^{2}
Structured QR  k[\frac{46}{3}n^{3} + (15R + 30)n^{2}] 4kn^{2}
Normal equations k[\frac{62}{3}n^{3} + 12Rn^{2}]  5kn^{2}
DECOMP/SOLVE  k[\frac{14}{3}n^{3} + 4Rn^{2}]  3kn^{2}

6. Computational results. Versions of the structured QR algorithm have been implemented on the Alliant FX/8 vector multiprocessor at Argonne National Laboratory and on the CRAY Y-MP at the North Carolina Supercomputing Center. Performance comparisons were made with

• a row partial pivoting code (two versions, for separated and non-separated boundary conditions), and

• the DECOMP and SOLVE routines from the PASVA codes [13]. The DECOMP routine uses alternate row and column pivoting (as does the algorithm described in Varah [19] and in §2) but always eliminates by rows.

Approximate operation counts and storage requirements for these algorithms and the normal equations method of Ascher and Chan [1] are given in Tables 1 and 2, for separated and non-separated boundary conditions, respectively. In tabulating storage requirements, it has been assumed that intermediate information generated during the factorization — namely, the multipliers and Householder vectors — is stored for possible later use with a different right-hand side. In general, the structured QR algorithms require the most operations. This result is not surprising, since it is well known that orthogonal factorization of dense matrices requires about twice as much work as Gaussian elimination. Moreover, in the case of separated end conditions, structured QR generates fill-in that is avoided by the elimination-based methods, and this adds further to the operation count. On the other hand, the operation counts and storage requirements for strutured QR are not affected if the end conditions are non-separated rather than separated, while for the other methods, they increase substantially. The method based on normal equations requires about the same amount of work as structured QR; however, as we noted earlier, it is less stable. We add the caveat that operation counts are notoriously bad predictors of run time for factorizations of narrow-banded matrices. For small n, much of the CPU time is taken up with non-numerical operations. This is borne out by our results, which show that the QR algorithm does not do as badly as predicted in serial mode.

The linear system solvers described above have been incorporated into both a multiple shooting and finite difference (box method) framework. The dverk code from netlib was used to solve the IVPs on each interval of multiple shooting, with the global error tolerance parameter set to 10^{−10}. A user-specified number k of equallyspaced intervals is used for both methods. For practical codes, the choice of the number of intervals and their lengths (and efficient parallel implementation of this) are important issues, but we focus here on the linear algebra, which typically is the most computationally intensive part of a finite difference-based code. In accord with the theoretical results of §2 and §4, virtually no difference was noted between the stability properties of structured QR and row-pivoted LU. As evidence of this, we quote results for the following test problem:

Problem 1 (Ascher and Chan [1]) a = 0, b = 1, n = 2,

y′(t) = [ \begin{matrix}−λ cos 2ωt & ω + λ sin 2ωt \\ −ω + λ sin 2ωt & λ cos 2ωt \end{matrix} ] y(t) + f(t),

B_{a} = [ \begin{matrix}1 & 0 \\ 0 & 0 \end{matrix} ] , B_{b} = [ \begin{matrix}0 & 0 \\ 1 & 0 \end{matrix} ] ,

with f(t) and d chosen so that y(t) = e^{t}(1, 1)^{T}.

A fundamental solution is

Y (t) = [ \begin{matrix}cos ωt & sin ωt \\ − sin ωt & cos ωt \end{matrix} ] [ \begin{matrix}e^{−λt} & 0 \\ 0 & e^{λt} \end{matrix} ] ,

so clearly there exist one growing and one decaying mode. The problem was solved using both multiple shooting and the box method for λ = 200 and ω = 1. Four different algorithms were used to solve the linear system, namely,

• SQR-1 — one-level structured QR;

• ROWPP — LU factorization with row partial pivoting;

• DECOMP — the DECOMP and SOLVE routines from PASVA; and

• COMPACT — compactification, as implemented in the codes D4/S4 described

in [21]. Tables 3 and 4 show the maximum error in the first component of the computed result. Because of its failure to decouple the fundamental solution modes, compactification performs poorly. The accuracy of box method solutions is limited by discretization error, while the multiple shooting solutions are accurate up to the conditioning of the discrete system and the tolerance imposed on the IVP solver. The DECOMP code gives accurate results here because the end conditions are separated. When this is not the case, as in problem 3 below, DECOMP is known to be unstable.

To test the relative speed of the linear solvers, two further problems from the literature were used in addition to problem 1. These were

Problem 2 (Brown and Lorenz [5]) a = −1, b = 1, n = 4,

−ϵy′′ − \frac{t}{2}y′ + \frac{t}{2}z′ + z = ϵπ^{2} cos πt + \frac{1}{2}πt sin πt,
ϵz′′ = z
y(−1) = −1  y(1) = e^{−2/\sqrt{ϵ}}
z(−1) = 1  z(1) = e^{−2/\sqrt{ϵ}}.

Table 3 Box method, error in first component of computed solution for Problem 1 for four different

linear system solvers

k = 16 k = 64 k = 1024
κ̂  .17(+2) .51(+1) .21(+1)
ROWPP .21(-2) .10(-3) .32(-6)
DECOMP .22(-2) .10(-3)  .32(-6)
SQR-1 .21(-2) .10(-3) .32(-6)
COMPACT .22(-2) .93(+27) .16(+72)

Table 4 Multiple shooting, error in first component of computed solution for Problem 1 for four different

linear system solvers

k = 16  k = 32 k = 128
κ̂  .85(+6) .17(+4) .13(+2)
ROWPP  .45(-3)  .67(-6)  .64(-7)
DECOMP  .45(-3)  .67(-6)  .64(-7)
SQR-1  .45(-3)  .67(-6)  .64(-7)
COMPACT .51(+72) .21(+72) .11(+72)

(We use ϵ = .001.)

Problem 3 (Mattheij [15]) a = 0, b = π, n = 3,

y′(t) = [ \begin{matrix}1 − 19 cos 2t & 0 & 1 + 19 sin 2t \\ 0 & 19 & 0 \\ −1 + 19 sin 2t & 0 & 1 + 19 cos 2t \end{matrix} ] y(t)

+ e^{t} [ \begin{matrix}−1 + 19(cos 2t − sin 2t) \\ −18 \\ 1 − 19(cos 2t + sin 2t) \end{matrix} ] ,

y_{1}(0) = 1
y_{3}(0) + y_{3}(π) = 1 + e^{π}
y_{2}(0) + y_{2}(π) = 1 + e^{π}.

The solution is y(t) = e^{t}(1, 1, 1)^{T}.

Problems 1 and 2 have separated end conditions, while two of the three end conditions for problem 3 are non-separated. We report on five cases (two different values of k were tried for problems 1 and 2, and values of λ = 1 and ω = 50 were used in problem 1). Table 5 gives condition estimates for the multiple shooting and finite difference matrices.

Results from “scalar” implementations on one processor of the Alliant FX/8 are shown in Table 6. We have tabulated the times required to solve the linear systems. The -Og compiler option was used with each code, so the vector processing capabilities of the Alliant were not used. In addition to the linear solvers already mentioned, we tested SQR-CR, which was the cyclic-reduction variant of structured QR. Note that the SQR codes typically take two to three times as long as ROWPP, though the penalty is

Table 5

Dimensions of the five test cases, and conditioning of the multiple shooting and finite difference matrices

Problem n k ms conditioning fd conditioning
1a  2 64  .41(+2)  .17(+2)
1b  2 1024  .30(+2)  .30(+2)
2a  4 64  .37(+9)  .83(+3)
2b  4 1024  .26(+2)  .17(+2)
3  3 1024  .20(+1)  .34(+1)

Table 6

Alliant FX/8, one-processor timings for linear system solvers (times in seconds)

Problem ROWPP  DECOMP  SQR-1  SQR-CR

1a  .041  .081  .071  .098
1b  .439  .660  .813  1.04
2a  .094  .181  .220  .312
2b  1.32  1.95  3.22  4.55
3  1.23  1.44  1.72  2.34

Table 7

CRAY Y-MP, one-processor timings for linear system solvers. Vectorized code (times in milliseconds)

Problem ROWPP  DECOMP  SQR-1  SQR-CR

1a  2.19  3.05  8.63  1.45
1b  34.6  48.3  139.  10.9
2a  5.32  8.38  21.1  6.40
2b  84.4  133.  341.  51.6
3  136.  116.  232.  26.3

Table 8

Alliant FX/8, eight-processor timings for linear system solvers (times in seconds)

  Problem ROWPP  DECOMP  SQR-2

1a  .029  .072  .031
1b  .367  .697  .136
2a  .053  .150  .067
2b  .739  1.67  .463
3  .685  1.39  .262

Table 9

Alliant FX/8, Ratio of times for ROWPP (one-processor) to times for SQR-2 (eight processors)

Problem Speedup
1a  1.3
1b  3.2
2a  1.4
2b  2.9
3  4.7

much smaller when the end conditions are not separated (as in problem 3). In either case, the overhead for using structured QR is not as great as the operation counts in Tables 1 and 2 would suggest.

Timings for a vectorized implementation on one processor of a CRAY Y-MP are shown in Table 7. In general, the SQR-CR code becomes very competitive, particularly when n = 2 or 3, k is large, and/or the end conditions are not separated. This code performs extremely well on problems 1b and 3. When n = 4 (problem 2), the small amount of vectorization that occurs in the other codes lessens the advantage of SQR-CR, while in problems 1a and 2a the value of k makes the overall computational task to small to benefit from vectorization.

Table 8 gives results for an eight-processor parallel implementation on the Alliant FX/8. The -Ogc option was used during compilation. Here, SQR-2 refers to the two-level version of structured QR, in which the original system is broken into eight partitions of equal size, which are factorized concurrently. On the largest problem, the parallel efficiency of structured QR (measured by comparing serial SQR-1 to parallel SQR-2) is 87% — quite acceptable, given that the solution of the reduced system is an unavoidable bottleneck. The efficiency improves further for still larger problems. Defining speedup as the ratio of the one-processor time for the best serial algorithm to the eight-processor time for the best parallel algorithm, we see, from Table 9, that in three of the five cases good parallel efficiency is attained. The remaining two problems were too small for parallelism to have much effect.

Comparing Tables 6 and 8, it can be seen that ROWPP and DECOMP also speed up a little when extra processors are available. This is because the Alliant is a sharedmemory machine. It is important to note that on the current generation of messagepassing machines, these algorithms will not benefit from multiprocessing unless n is large enough that rows or columns within each block can profitably be distributed around the processor array. This is unlikely to happen until n is at least 50 or 100. On the other hand, efficient implementations of multilevel SQR on these machines will be possible for much more typical problem sizes.

To summarize, we conclude that the structured QR codes are useful in the following circumstances:

• when the computational task of solving the linear equations is substantial enough to benefit from vectorization or parallelism;

• especially, when the end conditions are not separated;

• on a vector processor, when the value of n is too small (say, only 2 or 3) to allow efficient vectorized factorization of n × n blocks;

• on the current generation of distributed-memory multiprocessors, unless n is very large and the number of processors is very small;

• on a shared-memory multiprocessor, unless n is quite large (say, greater than 8) and there are fewer than four processors.

A. Appendix A. We start with a result which is similar to [8, Theorem 3.3.1]:

Lemma A.1. Suppose that alternate row and column elimination, without pivoting, is applied to an N × N matrix A with bandwidth b_{w} to produced computed factors L̂, B̂, Û. Assume that N, b_{w} and the unit roundoff error u satisfy

(30)  Nu < 0.1
(31)  Nu(2 + 1.06b_{w} + 2.12b_{w}u) < 0.5.

18 STEPHEN J. WRIGHT

Then

L̂B̂Û = A + H,

where |H| ≤ c_{3}(N − 1)u{|A| + |L̂||B̂||Û|} (32) and

c_{3} = 5 + 3b_{w}.

Proof. The result is trivially true for N = 1. Suppose for induction that (30) holds for matrices of size up to N − 1. Let

A = [ \begin{matrix}α & w^{T} \\ v & A_{1} \end{matrix} ] ,

where α ∈ R, A_{1} ∈ R^{(N−1)×(N−1)}, etc., and suppose that row elimination will be used to eliminate v. We compute ẑ = \frac{1}{α}v + f, |f| ≤ \frac{|v|}{α}u, (33) Â_{1} = A_{1} − ẑw^{T} + F, |F | ≤ 2u(|A_{1}| + |ẑ||w|^{T}), (34) It follows immediately from (34) that |Â_{1}| ≤ (1 + 2u)(|A_{1}| + |ẑ||w|^{T}). (35) An LBU factorization of Â_{1} is then performed, yielding

L̂_{1}B̂_{1}Û_{1} = Â_{1} + H_{1},

with |H_{1}| ≤ c_{3}(N − 2)u{|Â_{1}| + |L̂_{1}||B̂_{1}||Û_{1}|}. (36) The calculated factors of A are therefore

L̂ = [ \begin{matrix}1 & 0 \\ ẑ & L̂_{1} \end{matrix} ] ,  B̂ = [ \begin{matrix}α & ŵ^{T} \\ 0 & B̂_{1} \end{matrix} ] ,  Û = [ \begin{matrix}1 & 0 \\ 0 & Û_{1} \end{matrix} ] ,

where ŵ is the computed solution of the system Û_{1}^{T} w̃ = w. Defining b̄_{w} = 1.06b_{w}, and noting that Û_{1}^{T} has lower bandwidth b_{w}, it is easy to show that ŵ exactly solves

(Û_{1}^{T} + δ_{U})ŵ = w,

where

|δ_{U}| ≤ b̄_{w}u|Û_{1}|^{T} .

Hence |Û_{1}^{T} ŵ − w| = |δ_{U} ŵ| ≤ |δ_{U}||ŵ| ≤ b̄_{w}u|Û_{1}|^{T} |ŵ|, (37)

PARALLEL ALGORITHMS FOR BVP 19

and so |w| ≤ (1 + b̄_{w}u)|Û_{1}|^{T} |ŵ|. (38) Now

Â_{1} = A_{1} − ẑw^{T} + F
⇒ ẑŵ^{T} Û_{1} + L̂_{1}B̂_{1}Û_{1} = A_{1} + F + H_{1} − ẑ[w − Û_{1}^{T} ŵ]^{T} .

Combining this with (33),

L̂B̂Û = A + [ \begin{matrix}0 & ŵ^{T} Û_{1} − w^{T} \\ αf & F + H_{1} − ẑ[w − Û_{1}^{T} ŵ]^{T} \end{matrix} ] . (39)

Now, combining (34), (35), (36) and (38), we find that

|F + H_{1} − ẑ[w − Û_{1}^{T} ŵ]^{T} |

≤  |F | + |H_{1}| + |ẑ||w − Û_{1}^{T} ŵ|^{T}
≤  [c_{3}(N − 2)u(1 + 2u) + 2u]|A_{1}|
+[c_{3}(N − 2)u(1 + 2u)(1 + b̄_{w}u) + 2u(1 + b̄_{w}u) + b̄_{w}u]|ẑ||ŵ|^{T} |Û_{1}|

+c_{3}(N − 2)u|L̂_{1}||B̂_{1}||Û_{1}|. (40) We now show that c_{3}(N − 2)u(1 + 2u)(1 + b̄_{w}u) + 2u(1 + b̄_{w}u) + b̄_{w}u ≤ c_{3}(N − 1)u. (41) This is equivalent to

c_{3}(N − 2)u(2u + b̄_{w}u + 2b̄_{w}u^{2}) + (2u + b̄_{w}u + 2b̄_{w}u^{2}) ≤ c_{3}u
⇔ c_{3}[1 − (N − 2)u(2 + b̄_{w} + 2b̄_{w}u)] ≥ (2 + b̄_{w} + 2b̄_{w}u).

Since by assumption (31),

1 − (N − 2)u(2 + b̄_{w} + 2b̄_{w}u) ≥ 1/2,

inequality (41) will hold provided that c_{3} ≥ 2(2 + b̄_{w} + 2b̄_{w}u). (42) From (30), and using b_{w} < N, it is clear that 2b̄_{w}u ≤ 0.5, so (42) follows trivially.

Since the left-hand side of (41) is the largest of the three coefficients in (40), we can combine (40) and (41) to obtain |F + H_{1} − ẑ[w − Û_{1}^{T} ŵ]^{T} | ≤ c_{3}(N − 1)u{|A_{1}| + |ẑ||ŵ|^{T} |Û_{1}| + |L̂_{1}||B̂_{1}||Û_{1}|}. (43) Combining (39), (33), (37) and (43), we therefore find that

|A − L̂B̂Û|

= |[ \begin{matrix}0 & ŵ^{T} Û_{1} − w^{T} \\ αf & F + H_{1} − ẑ[w − Û_{1}^{T} ŵ]^{T} \end{matrix} ]|
≤ c_{3}(N − 1)u \{[ \begin{matrix}|α| & |w|^{T} \\ |v| & |A_{1}| \end{matrix} ] + [ \begin{matrix}1 & 0 \\ |ẑ| & |L̂_{1}| \end{matrix} ] [ \begin{matrix}|α| & |ŵ|^{T} \\ 0 & |B̂_{1}| \end{matrix} ] [ \begin{matrix}1 & 0 \\ 0 & |Û_{1}| \end{matrix} ]\}
= c_{3}(N − 1)u{|A| + |L̂||B̂||Û|}.

This proves the result for the case in which row elimination is used at stage N. When column elimination is used instead, the proof is analogous.

Lemma A.2. If the alternate row and column elimination is used to solve (5) or (8), and the assumptions of Lemma A.1 hold, with N = (k + 1)n and b_{w} = 2n, then the computed solution ŝ is the exact solution of the perturbed system

(A + E)ŝ = f,

where

∥E∥_{∞} ≤ (k + 2)n(5 + 6n)(1 + n^{2}2^{n−1})∥A_{DP} ∥_{∞}u.

Proof. Note first that the pivoting does not alter the sparsity structure of A. We can, therefore, view alternate row and column elimination as being applied to P^{T} AΠ^{T} (where P and Π are permutation matrices) to produce a computed factorization

P^{T} AΠ^{T} + H = L̂B̂Û.

It is easy to show that the procedure leading to ŝ results in the following sequence of perturbed problems:

(L̂ + δ_{L})ŵ = P^{T} f,  |δ_{L}| ≤ (1.06)nu|L̂|,
(B̂ + δ_{B})v̂ = ŵ,  |δ_{B}| ≤ 2(1.06)nu|B̂|,
(Û + δ_{U})ẑ = v̂,  |δ_{U}| ≤ (1.06)nu|Û|,
ŝ = Π^{T} ẑ.

(The bounds on |δ_{L}|, |δ_{B}|, |δ_{U}| are a consequence of the maximum number of nonzeros in each row of L̂, B̂, Û, respectively.) Hence

P^{T} EΠ^{T} = H + (L̂ + δ_{L})(B̂ + δ_{B})(Û + δ_{U}) − L̂B̂Û
⇒ ∥E∥_{∞} ≤ ∥H∥_{∞} + [4(1.06)nu + 5(1.06)^{2}n^{2}u^{2} + 2(1.06)^{3}n^{3}u^{3}]∥L̂∥_{∞}∥B̂∥_{∞}∥Û∥_{∞}.

It follows from (30) that nu ≤ 0.1, so the coefficient of ∥L̂∥_{∞}∥B̂∥_{∞}∥Û∥_{∞} can be bounded above by 5nu. Since element growth in B is bounded by 2^{n−1}, and since all entries in L̂ and Û are bounded by 1, we have

∥B̂∥_{∞} ≤ 2^{n−1}∥P^{T} AΠ^{T} ∥_{∞} ≤ 2^{n−1}∥A∥_{∞}, ∥L̂∥_{∞} ≤ n, ∥Û∥_{∞} ≤ n.

Combining these observations with the result of Lemma A.1, we obtain

∥E∥_{∞} ≤ (5 + 6n)(k + 1)nu{∥A∥_{∞} + ∥L̂∥_{∞}∥B̂∥_{∞}∥Û∥_{∞}} + 5nu∥L̂∥_{∞}∥B̂∥_{∞}∥Û∥_{∞}
≤ {(5 + 6n)(k + 1)n + [(5 + 6n)n(k + 1) + 5n]n^{2}2^{n−1}}∥A∥_{∞}u
≤ (k + 2)n(5 + 6n)(1 + n^{2}2^{n−1})∥A∥_{∞}u.

as required.

B. Appendix B. We start by stating two results on the rounding error due to Householder reduction. These results are similar to those in Lawson and Hanson [10, pp. 85–89] and Wilkinson [20, pp. 157–162]. They differ from Lawson and Hanson’s results in that the O(u^{2}) term is explicitly accounted for at every stage, and from Wilkinson’s in that we do not assume double-precision accumulation of inner products. Since the proofs are tedious and do not offer any new insight, they are omitted.

Lemma B.1. Suppose that an m_{1} × m_{2} matrix X is multiplied by an m_{1} × m_{1} Householder reflector Q. Then provided that

(6m_{1} + 18)u ≤ 0.1,

the computed result Y satisfies

Y = Q(X + E),

where

∥E∥_{F} ≤ (7m_{1} + 42)u∥X∥_{F} .

Lemma B.2. If Q is a product of r Householder reflectors whose effect is to introduce zeros into the subdiagonals into the first r columns of the m_{1} × m_{2} matrix X, then provided that

(7m_{1} + 42)ru ≤ 0.1,

the computed result Y satisfies

Y = Q(X + E),

where

∥E∥_{F} ≤ (8m_{1} − 4r + 51)ru∥X∥_{F} .

For the purpose of this Appendix, it is simplest to view the structured factorization process as the application of k − 1 orthogonal transformation matrices Q_{1}, Q_{2}, …, Q_{k−1} to a row- and column-reordered version of A_{D} (and the right-hand side f_{D}), followed by the application of another two matrices Q_{k} and Q_{k+1} to effect the final reduction of (25). (Q_{k} and Q_{k+1} reduce the first and last n columns of the coefficient matrix in (25), respectively.) Each of the Q_{j}, j = 1, …, k − 1 are products of n Householder reflectors, and each operates on only a small part of the matrix that it multiplies: to be precise, a 2n × 3n submatrix. Since we wish to reduce (7) to (25), it follows that exactly k − 1 such transformations are needed.

Proof. (Lemma 4.1) Let Â_{D,i+1} be the transformed version of A_{D} after i stages of

the structured factorization, and Â_{D,1} = A_{D}. Let ˆ Â_{D,i+1} be the submatrix which is actually affected at stage i, by the matrix Q_{i} from (22), and let Q̂_{i} be the orthogonal matrix which is obtained by embedding Q_{i} into a (k+1)n-dimensional identity matrix. For i = 1, …, k + 1, we have from Lemma B.2 that Â_{D,i} and Â_{D,i+1} are related by

Â_{D,i+1} = Q̂_{i}(Â_{D,i} + E_{i}),

where E_{i} consists of the 2n × 3n “error” submatrix corresponding to the factorization of ˆ Â_{D,i}, padded out with zeros to dimension (k + 1)n. Hence

∥E_{i}∥_{F} ≤ (8(2n) − 4n + 51)nu∥ˆÂ_{D,i}∥_{F} ≤ (12n + 51)nu∥Â_{D,i}∥_{F} .

Hence

∥Â_{D,i+1}∥_{F} ≤ [1 + (12n + 51)nu]∥Â_{D,i}∥_{F} ≤ [1 + (12n + 51)nu]^{i}∥A_{D}∥_{F} .

The errors made in stages k and k + 1 are bounded in the same way, since the submatrices affected at these stages are no larger than those affected at the earlier stages. Under the assumption (12n + 51)n(k + 1)u ≤ 0.1, we therefore obtain

∥Â_{D,i+1}∥_{F} ≤ [1 + (1.06)(12n + 51)n(k + 1)u]∥A_{D}∥_{F} ≤ (1.106)∥A_{D}∥_{F} ,

for i = 1, …, k + 1. Hence

∥E_{i}∥_{F} ≤ (1.106)(12n + 51)nu∥A_{D}∥_{F} .

The final (upper triangular) matrix is Â_{D,k+2}, which satisfies

Â_{D,k+2} = Q̂_{k+1}(Â_{D,k+1} + E_{k+1})
= Q̂_{k+1}Q̂_{k}(Â_{D,k} + E_{k}) + Q̂_{k+1}E_{k+1}
⋮
= Q̂_{k+1} … Q̂_{1}A_{D} + Q̂_{k+1} … Q̂_{1}E_{1} + Q̂_{k+1} … Q̂_{2}E_{2} + … + Q̂_{k+1}E_{k+1}
= Q̂[A_{D} + E_{D}].

Here Q̂ = Q̂_{k+1}Q̂_{k} … Q̂_{1} is orthogonal, and E_{D} satisfies the bound

∥E_{D}∥_{F} ≤ ∑_{i=1}^{k+1} ∥E_{i}∥_{F} = (1.106)(12n + 51)n(k + 1)u∥A_{D}∥_{F} .

Similarly, application of Q̂ to the right-hand side f_{D} results in a computed vector f̂_{D} which satisfies

f̂_{D} = Q̂[f_{D} + δf_{D}],

where

∥δf_{D}∥_{2} ≤ (1.106)(12n + 51)n(k + 1)u∥f_{D}∥_{2}

(since ∥v∥_{F} = ∥v∥_{2} when v is a vector.) Finally, back-substitution is used on the system with coefficient matrix Â_{D,k+2} and right-hand side f̂_{D}. The computed solution satisfies

(Â_{D,k+2} + E_{S})ŝ = f̂_{D}, (44)

where, since Â_{D,k+2} has at most 3n nonzeros per row,

∥E_{S}∥_{F} ≤ 3(1.06)nu∥Â_{D,k+2}∥_{F}
≤ 3(1.06)nu[1 + (1.106)(12n + 51)n(k + 1)u]∥A_{D}∥_{F}
≤ 3(1.06)(1.1106)nu∥A_{D}∥_{F} .

Substituting in (44),

[Q̂[A_{D} + E_{D}] + E_{S}] ŝ = Q̂[f_{D} + δf_{D}]
⇒ (A_{D} + E_{D} + Q̂^{T} E_{S})ŝ = (f_{D} + δf_{D}).

Defining

δA_{D} = E_{D} + Q̂^{T} E_{S},

we have

∥δA_{D}∥_{F} ≤ ∥E_{D}∥_{F} + ∥E_{S}∥_{F}
≤ [(1.106)(12n + 51)n(k + 1) + 3(1.06)(1.1106)n] u∥A_{D}∥_{F}
≤ (1.106)(12n + 51)n(k + 2)u∥A_{D}∥_{F} ,

as required.

REFERENCES

[1] U. M. Ascher and P. S. Y. Chan, On parallel methods for boundary value odes, Computing, 46 (1991), pp. 1–17.

[2] U. M. Ascher and R. M. M. Mattheij, General framework, stability and error analysis for numerical stiff boundary value problems, Numerische Mathematik, 54 (1988), pp. 355–372.

[3] U. M. Ascher, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, 1988.

[4] Z. Bohte, Bounds for rounding errors in the Gaussian elimination for band systems, Journal of the Institute of Mathematics and its Applications, 16 (1975), pp. 133–142.

[5] D. L. Brown and J. Lorenz, A high order method for stiff boundary value problems with turning points, SIAM Journal on Scientific and Statistical Computing, 8 (1987), pp. 790– 805.

[6] A. K. Cline, C. B. Moler, G. W. Stewart, and J. H. Wilkinson, An estimate for the condition number of a matrix, SIAM Journal on Numerical Analysis, 16 (1979), pp. 368– 375.

[7] J. C. Diaz, A. Fairweather, and P. Keast, FORTRAN packages for solving certain almost block diagonal linear systems by modified alternate row and column elimination, ACM Transactions on Mathematical Software, 9 (1983), pp. 358–375.

[8] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, MD, second ed., 1989.

[9] H. B. Keller, Accurate difference methods for two-point boundary value problems, SIAM Journal on Numerical Analysis, 11 (1974), pp. 305–320.

[10] C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, NJ, 1974.

[11] M. Lentini, Parallel solution of special large block tridiagonal systems: Tpbvp. Manuscript, 1989.

[12] M. Lentini, M. R. Osborne, and R. D. Russell, The close relationships between methods for solving two-point boundary value problems, SIAM Journal on Numerical Analysis, 22 (1985), pp. 280–309.

[13] M. Lentini and V. Pereyra, An adaptive finite difference solver for nonlinear two-point boundary value problems with mild boundary layers, SIAM Journal on Numerical Analysis, 14 (1977), pp. 91–111.

[14] R. M. M. Mattheij, The conditioning of linear boundary value problems, SIAM Journal on Numerical Analysis, 19 (1982), pp. 963–978.

[15] , Decoupling and stability of algorithms for boundary value problems, SIAM Review, 27 (1985), pp. 1–44.

[16] M. R. Osborne, Aspects of the numerical solution of boundary value problems with separated boundary conditions. Manuscript, 1978.

[17] M. Paprzycki and I. Gladwell, Solving almost block diagonal systems on parallel computers, Parallel Computing, 17 (1991), pp. 133–153.

[18] V. Pereyra, Iterated deferred corrections for nonlinear boundary value problems, Numerische Mathematik, 8 (1968), pp. 111–125.

[19] J. M. Varah, Alternate row and column elimination for solving certain linear systems, SIAM Journal on Numerical Analysis, 13 (1976), pp. 71–75.

[20] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, London, 1965.

[21] S. J. Wright and V. Pereyra, Adaptation of a two-point boundary value problem solver to a vector-multiprocessor environment, SIAM Journal on Scientific and Statistical Computing, 11 (1990), pp. 425–449.