PARALLEL STABLE COMPACTIFICATION FOR ODE WITH PARAMETERS AND MULTIPOINT

CONDITIONS^{∗}

R. M. M. MATTHEIJ^{†}  S. J. WRIGHT^{‡}

June 30, 1993

Abstract

Many algorithms for solving ordinary differential equations with parameters and multipoint side conditions give rise to systems of linear algebraic equations in which the coefficient matrices have a bordered block diagonal structure. In this paper, we show how these problems can be solved by using parallel algorithms based on stabilized compactification.

1 Introduction

We discuss stable algorithms for solving algebraic linear systems arising from ordinary differential and difference equations. These equations may include some global parameters, and their solutions are subject to multipoint side conditions. To make the discussion in subsequent sections clearer, we group the problems to be addressed into three categories:

Two-point boundary value problems:

ẋ = A(t)x + c(t),  t ∈ [a, b],  (1)
M_{a}x(a) + M_{b}x(b) = d, x(t), c(t), d ∈ R l ^{n};

Parametrized problems with multipoint conditions (see, for example, soluton formulations discussed in Rinzel and Terman [11]):

ẋ = A(t)x + C(t)λ + c(t),  t ∈ [a, b],  (2)
∑_{j=1}^{p} M_{i}x(τ_{j}) + Nλ = d, a ≤ τ_{1} ≤ ⋯ ≤ τ_{p} ≤ b,
x(t), c(t) ∈ R l ^{n},  λ ∈ R l ^{m},  M_{i} ∈ R l ^{(m+n)×n},  N ∈ R l ^{(m+n)×m};

Parameter identification problems (see, for example, chemical rate constant identification as in Bock [3]):

min ∑_{j=1}^{p} ∥M_{j}x(τ_{j}) + N_{j}λ − d_{j}∥^{2} + ∥N_{0}λ − d_{0}∥^{2},  (3)
ẋ = A(t)x + C(t)λ + b(t), t ∈ [a, b]
x(t) ∈ R l ^{n}, λ ∈ R l ^{m}, M_{j} ∈ R l ^{n_{j}×n}, N_{j} ∈ R l ^{n_{j}×m}.

∗This research was supported by the Applied Mathematical Sciences subprogram of the Office of Energy Research, U.S.

Department of Energy, under Contract W-31-109-Eng-38.

†Department of Mathematics and Computing Science, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, the Netherlands.

‡Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439, USA.

Although these categories appear to be listed in order of increasing generality, the ordering is not strict: problems of the form (2) can be recast in the form (1) (see, for example, Ascher and Russell [2]). For purposes of efficiency, it is usually best to use the original form when solving such problems.

Each of these problem categories has a corresponding “discrete” form, in which the differential equation is replaced by a difference equation. We categorize the discrete problems similarly:

A_{i}x_{i} + B_{i}x_{i+1} = c_{i},  i = 1, ⋯, N,  (4)
M_{a}x_{1} + M_{b}x_{N+1} = d, x_{i}, c_{i}, d ∈ R l ^{n};

A_{i}x_{i} + B_{i}x_{i+1} + C_{i}λ = c_{i},  i = 1, ⋯, N,  (5)
∑_{i=1}^{N+1} M_{i}x_{i} + Nλ = d,
x_{i}, c_{i} ∈ R l ^{n},  λ ∈ R l ^{m},  M_{i} ∈ R l ^{(m+n)×n},  N ∈ R l ^{(m+n)×m};

min ∑_{i=1}^{N+1} ∥M_{i}x_{i} + N_{i}λ − d_{i}∥^{2} + ∥N_{0}λ − d_{0}∥^{2},  (6)
A_{i}x_{i} + B_{i}x_{i+1} + C_{i}λ = c_{i}, i = 1, ⋯, N,
x_{i} ∈ R l ^{n}, λ ∈ R l ^{m}, M_{i} ∈ R l ^{n_{i}×n}, N_{i} ∈ R l ^{n_{i}×m}.

The “discrete” categories can be obtained by applying algorithms such as multiple shooting or finite differencing to the corresponding “continuous” categories. The widely used collocation algorithm also gives rise to such discrete problems after a condensation step has been applied.

We will not address the special difficulties caused by stiffness and singularly perturbed problems. Neither will we specifically address nonlinear problems, though it is well known that algorithms for such problems often require the solution of linear systems like those described above as a core operation.

In the next section, we will review the concept of conditioning for the continuous problems above, and its implications for the conditioning of the linear systems to be solved in the discrete problems (4), (5), and (6). In §3, a parallel algorithm for (4) is described and is contrasted with previously proposed algorithms (e.g., Wright [13, 12]). This algorithm is extended to problems of the form (5) and (6) in §5. Finally, some numerical results are presented in §6.

2 Conditioning of the Problem and Structure of the Solution Spaces

The various problems defined in §1 relate to solution spaces of potentially different types. This variety of solution types influences the possibilities for parallelization that will be discussed in the next section. Here we characterize these solution spaces; that is, we describe certain growth properties of homogeneous modes in the fundamental solution Φ(t). Throughout, we assume that the problem is well conditioned. By this we mean that small perturbations of inhomogeneities in the ODE and the side conditions manifest themselves in the solution as quantities that are only a moderate constant larger in norm than the original perturbations, that is,

∥δx∥ ≤ κ max(|d|, ∥c∥),  (7)

where ∥.∥ is a suitable function space norm, in either the continuous or discrete setting, and |.| is any Hölder norm.

In the case of the two-point boundary value problem (1), well-conditioning implies that the underlying solution space is dichotomic; that is, there exists a projection P and a moderate constant κ̃ such that

\{ \begin{matrix}|Φ(t)P Φ^{−1}(s)| ≤ κ̃, & t > s, \\ |Φ(t)(I − P)Φ^{−1}(s)| ≤ κ̃, & t < s. \end{matrix}  (8)

(See de Hoog and Mattheij [5].) Hence, Φ can be properly split into one set of solution modes that do not (significantly) increase and another set of modes that do not (significantly) decrease.

If we have multipoint side conditions (m = 0 in (2)), then well-conditioning allows for a “switch” of growth behavior at any of the internal conditions. This property was referred to by de Hoog and Mattheij [5] as polychotomy. Specifically, we can choose projections P_{1}, ⋯, P_{M} (M ≤ min(p, n)) such that ∑_{j=1}^{M} P_{j} = I, P_{i}P_{j} = P_{j}P_{i} = δ_{ij}P_{j}, and

\{ \begin{matrix}|Φ(t) ∑_{j=1}^{k} P_{j}Φ^{−1}(s)| ≤ κ̃, & β_{k} < s ≤ β_{k+1},  t > s, \\ |Φ(t) ∑_{j=k+1}^{M} P_{j}Φ^{−1}(s)| ≤ κ̃, & β_{k} < s ≤ β_{k+1},  t < s, \end{matrix}  (9)

where, again, κ̃ is a moderate constant and the switch points a = β_{1}, β_{2}, ⋯, β_{M+1} = b are a subset of {a, τ_{1}, ⋯, τ_{p}, b}. Polychotomy means that Φ is dichotomic on each subinterval (β_{k}, β_{k+1}) but that the dimension of the subspace of “nonincreasing” modes may become larger — by rank(P_{k}) — at each switch point β_{k}.

A more or less complementary result to the foregoing holds for “pure” parameter problems (for which p = 2, a = τ_{1} and b = τ_{2} in (2)). Given well-conditioning, Mattheij [10] showed that there exist a moderate constant κ̃, projections P_{1}, ⋯, P_{M} with (M ≤ min(m, n)) satisfying the conditions described above, and a set of switch points α_{k}, k = 1, ⋯ M + 1 with a = α_{1} < α_{2} < ⋯ < α_{M+1} = b such that

\{ \begin{matrix}|Φ(t) ∑_{j=k+1}^{M} P_{j}Φ^{−1}(s)| & ≤ & κ̃, & α_{k} < s ≤ α_{k+1},  t > s, \\ |Φ(t) ∑_{j=1}^{k} P_{j}Φ^{−1}(s)| & ≤ & κ̃, & α_{k} < s ≤ α_{k+1},  t < s \end{matrix}  (10)

Again, Φ is dichotomic on each subinterval, but the dimension of the subspace of nonincreasing modes may now become smaller at each switch point α_{k} as we move from left to right.

In the more general case of m ≠ 0 and p > 2 in (2), one should expect that any combination of (9) and (10) can occur. We will use the term “polychotomy” to refer to the multiple splitting of the solution space in this case as well.

Thus, we may expect for (2) a total of at most min(p, n) switch points at which the dimension of the subspace of nonincreasing modes may increase, and at most min(m, n) switch points at which this dimension may decrease. Note, however, that a typical “pure” mode of this subspace may undergo only one switch in behavior; otherwise, it would violate the local well-conditioning on each subinterval (β_{k}, β_{k+1}) or (α_{k}, α_{k+1}).

It should be noted that the potential switch points τ_{k} (from which the β_{k} are chosen) are known beforehand, while the α_{k} may be anywhere in the interval (a, b).

3 The Basic Algorithm

The “stabilized compactification” algorithm has been described in Ascher, Mattheij, and Russell [1, pp. 157–161]. In this section, we outline a parallel version and discuss some aspects of its implementation in a message-passing computational environment. We describe the method for the case of problem (4); this is the simplest of our discrete problems since, when it is well conditioned, the dimensions of the nonincreasing and nondecreasing subspaces remain constant across all the stages. We use ℓ to denote the dimension of the nonincreasing subspace. The dimension of the complementary nondecreasing subspace will therefore be (n − ℓ). For the moment, we also assume that the value of ℓ is known.

Suppose that the N + 1 indices in the problem (4) are split into p partitions of approximately equal size. This is done by choosing indices k_{1}, k_{2}, k_{3}, ⋯, k_{p+1} such that

k_{1} = 0, k_{p+1} = N,
k_{i+1} ≥ k_{i} + 2, i = 1, ⋯, p.

The j-th partition of the linear system is taken to be

[ \begin{matrix}A_{k_{j}+1} & B_{k_{j}+1} \\ & A_{k_{j}+2} & B_{k_{j}+2} \\ & & ... & ... \\ & & & A_{k_{j+1}} & B_{k_{j+1}} \end{matrix} ] [ x_{k_{j}+1} ; x_{k_{j}+2} ; ⋮ ; x_{k_{j+1}+1} ] = [ b_{k_{j}+1} ; b_{k_{j}+2} ; ⋮ ; b_{k_{j+1}} ] .  (11)

The algorithm starts by transforming the recurrence (11) to one involving upper triangular blocks. Then, explicit decoupling of nonincreasing and nondecreasing modes is used to find particular and fundamental solutions for this partition. Because of the decoupling, element growth in the fundamental solution is avoided.

Suppose for the moment that we can choose an n × n orthogonal matrix R_{k_{j}+1}^{(j)} that is effective in

decoupling these two sets of modes. (The choice of R_{k_{j}+1}^{(j)} will be discussed further below.) Starting with this matrix, the transformation process performs repeated QR (orthogonal) factorizations to find the following matrices, all of which are square with dimension n:

R_{i}^{(j)}, i = k_{j} + 2, ⋯, k_{j+1},  orthogonal,
Q_{i}^{(j)}, i = k_{j} + 1, ⋯, k_{j+1} − 1, orthogonal,
U_{i}^{(j)}, i = k_{j} + 1, ⋯, k_{j+1} − 1, upper triangular,
V_{i}^{(j)}, i = k_{j} + 1, ⋯, k_{j+1} − 1, upper triangular,

such that

A_{i}R_{i}^{(j)} = Q_{i}^{(j)}U_{i}  (12a)
Q_{i}^{(j)T} B_{i} = V_{i}R_{i+1}^{(j)T} .  (12b)

The recurrence

A_{i}x_{i} + B_{i}x_{i+1} = c_{i}, i = k_{j} + 1, ⋯, k_{j+1} − 1,

has now been transformed to

U_{i}z_{i} + V_{i}z_{i+1} = Q_{i}^{T} c_{i},  i = k_{j} + 1, ⋯, k_{j+1} − 1,  (13)

where

z_{i} = R_{i}^{T} x_{i}.

(Here and subsequently, the superscripts on R_{i}^{(j)}, Q_{i}^{(j)}, etc., will be dropped when their values are clear from the context.) From an implementation point of view, these operations create little fill-in. The matrices U_{i} and V_{i} can be stored in the upper triangles of the data structures formerly occupied by A_{i} and B_{i}. Most of the information needed to reconstruct the orthogonal matrices Q_{i} and R_{i} can be stored in the lower triangles, though an extra n-vector is also needed for each matrix.

We now construct a fundamental solution to (13), denoted by

{Φ_{i}}_{k_{j+1}}^{i=k_{j}+1},  Φ_{i} ∈ R l ^{n×n},  upper triangular,  (14)

that satisfies the boundary conditions

(Φ_{k_{j}+1})_{2,.} = [0 | I],  (Φ_{k_{j+1}})_{1,.} = [I | 0].  (15)

Here, (.)_{1,.} denotes the first (n − ℓ) rows of a matrix and (.)_{2,.} denotes the last ℓ rows. The remaining components of these two sequences are computed by back substitution with the matrix formed by the recurrence (13). The (.)_{2,.} components are calculated in a forward sweep and the (.)_{1,.} components in a reverse sweep. Specifically,

for i = k_{j} + 1, ⋯, k_{j+1} − 1

(Φ_{i+1})_{2,.} = −(V_{i})_{22}^{−1}[0 | (U_{i})_{22}(Φ_{i})_{22}];  (16)

for i = k_{j+1} − 1, ⋯, k_{j} + 1

(Φ_{i})_{1,.} = −(U_{i})_{11}^{−1} {(V_{i})_{11}(Φ_{i+1})_{1,.} + (V_{i})_{12}(Φ_{i+1})_{2,.} + (U_{i})_{12}(Φ_{i})_{2,.}} .  (17)

Here, (.)_{11} is the principal (n − ℓ) × (n − ℓ) submatrix of an n × n matrix; the three blocks (.)_{22}, (.)_{12}, and (.)_{21} complete a 2 × 2 block partitioning. Examination of (16) and (17) will confirm that the Φ_{i}^{(j)} are upper triangular. The homogeneous parts of (16) and (17) have coefficients (V_{i})_{22}^{−1}(U_{i})_{22} and −(U_{i})_{11}^{−1}(V_{i})_{11}, respectively. By our original assumption of well-conditioning, the recurrences are stable; they are effectively producing components of the nonincreasing modes (forward sweep) and nondecreasing modes (backward sweep), provided that the initial matrix R_{k_{j}} has been chosen in a suitable way.

The particular solution, denoted by

{ẑ_{i}}_{i=k_{j}+1}^{k_{j+1}},

is formed in a similar way. The boundary conditions are

(ẑ_{k_{j}+1})_{2,.} = 0,  (ẑ_{k_{j+1}})_{1,.} = 0.  (18)

If one uses the notation ĉ_{i} = Q_{i}^{T} c_{i}, the forward and reverse sweeps have the form

for i = k_{j} + 1, ⋯, k_{j+1} − 1

(ẑ_{i+1})_{2} = (V_{i})_{22}^{−1} {(ĉ_{i+1})_{2} − (U_{i})_{22}(ẑ_{i})_{2}};  (19)

for i = k_{j+1} − 1, ⋯, k_{j} + 1

(ẑ_{i})_{1} = (U_{i})_{11}^{−1} {(ĉ_{i})_{1} − (V_{i})_{11}(ẑ_{i+1})_{1} − (V_{i})_{12}(ẑ_{i+1})_{2} − (U_{i})_{12}(ẑ_{i})_{2}} .  (20)

We now deal with the issue of choosing the values of ℓ and R_{k_{j}+1}. Mattheij [9, p. 329] describes a heuristic that, applied to the present situation, would find R_{k_{j}+1} so that the diagonal elements of V_{i}^{−1}U_{i} (namely, (U_{i})_{ll}/(V_{i})_{ll}, l = 1, ⋯, n) appear in descending order for most i = k_{j} + 1, ⋯, k_{j+1} − 1. We generalize this heuristic slightly for two reasons. First, we usually do not have boundary condition information to guide the choice of R_{k_{j}+1}; and second, we wish to allow possible singularity of U_{i} and V_{i} (that is, of A_{i} and B_{i}.)

To make an initial guess at R_{k_{j}+1}, we make use of the generalized singular value decomposition. In the following result, which is a direct consequence of Theorem 8.7.4 of Golub and Van Loan [6], all matrices are assumed to be n × n.

Theorem 3.1 Provided that [A, B] has full rank n, there are orthogonal matrices P, Q, and W and an upper triangular matrix T such that

W^{T} T^{−T} AP = C, W^{T} T^{−T} BQ = S,

where C = diag(c_{1}, ⋯, c_{n}) and S = diag(s_{1}, ⋯, s_{n}), with 1 ≥ c_{1} ≥ ⋯ ≥ c_{n} ≥ 0 and 0 ≤ s_{1} ≤ ⋯ ≤ s_{n} ≤ 1, and c_{i}^{2} + s_{i}^{2} = 1, i = 1, ⋯, n.

Our initial guess for R_{k_{j}+1} is obtained as follows:

(i) Set A = A_{k_{j}+1}, B = B_{k_{j}+1}, and compute W, T, P, Q, C, and S as in Theorem 3.1.

(ii) Perform a QR factorization to obtain Z (n × n orthogonal) and Y (n × n upper triangular) such that T^{T} W = ZY.

(iii) Set R_{k_{j}+1} ← P, R_{k_{j}+2} ← Q, Q_{k_{j}+1} ← Z, U_{k_{j}+1} ← Y C, V_{k_{j}+1} ← Y S.

We note that this construction is valid since, by our assumption that (4) is well conditioned, [A_{k_{j}+1}, B_{k_{j}+1}] has full rank. It is also easy to check that equations (12) hold for i = k_{j} + 1. Moreover the diagonals of U_{k_{j}+1} and V_{k_{j}+1} are ordered properly, relative to each other. By this, we mean that

• any zero diagonal elements in V_{k_{j}+1} occur in the upper left of the matrix;

• any zero diagonals in U_{k_{j}+1} occur in the lower right; and

• for the indices l for which both (V_{k_{j}+1})_{ll} and (U_{k_{j}+1})_{ll} are nonzero, the ratio

(U_{k_{j}+1})_{ll}/(V_{k_{j}+1})_{ll} = c_{l}/s_{l}

decreases as l increases.

Using this initial choice of R_{k_{j}+1}, we now compute a few more successive Q_{i}, R_{i}, U_{i}, and V_{i} by using the formulae (12) and check to see whether the diagonals of U_{i} and V_{i} continue to appear in the proper order. If not, the columns of R_{k_{j}+1} are permuted in accordance with the present ordering, and the orthogonal transformation process is begun anew.

Before computing the fundamental and particular solutions, we need to determine the dimension ℓ of the nonincreasing subspace. When the boundary conditions are separated (that is, zero rows of M_{a} correspond to nonzero rows of M_{b}, and vice versa) and the problem is well conditioned, ℓ may just be taken to be the number of nonzero rows in M_{a}. Otherwise, we require the lower right ℓ × ℓ submatrix (V_{i})_{22} of V_{i} to dominate the corresponding submatrix (U_{i})_{22} of U_{i}. For “close calls,” we can use the test suggested in [9]: If the quantity

∑_{j=1}^{p} ∑_{i=k_{j}+1}^{k_{j+1}−1} log [(U_{i}^{(j)})_{ll} / (V_{i}^{(j)})_{ll}]

is negative, the l-th diagonal element is included in the (2, 2) partition. In a multiprocessor setting, the formation of this “global” sum requires interprocesssor communication. If vendor-supplied primitives for this type of operation are not available, it can be implemented easily by making a binary tree out the processors.

In partition j, each component x_{i} of the true solution can be expressed in terms of the fundamental and particular solutions in the following way:

z_{i} = Φ_{i}^{(j)}s_{j} + ẑ_{i}^{(j)},  x_{i} = R_{i}^{(j)}z_{i},  i = k_{j} + 1, ⋯, k_{j+1},  (21)

where the values of s_{j} ∈ R l ^{n}, j = 1, ⋯, p and x_{N+1} are determined by solving a reduced system. The boundary conditions contribute one block row to this system:

M_{a}x_{1} + M_{b}x_{N+1} = d
⟹ M_{a}R_{1}^{(1)}(Φ_{1}^{(1)}s_{1} + ẑ_{1}^{(1)}) + M_{b}x_{N+1} = d.  (22)

The remaining p blocks of the reduced system are obtained by considering the equations that bind adjacent partitions together. These are obtained by taking the last block equation from each partition. For j = 2, ⋯, p, we have

A_{k_{j}}x_{k_{j}} + B_{k_{j}}x_{k_{j}+1} = c_{k_{j}}
⟹ A_{k_{j}}R_{k_{j}}^{(j−1)} [Φ_{k_{j}}^{(j−1)}s_{j−1} + ẑ_{k_{j}}^{(j−1)}] + B_{k_{j}}R_{k_{j}+1}^{(j)} [Φ_{k_{j}+1}^{(j)}s_{j} + ẑ_{k_{j}+1}^{(j)}] = c_{k_{j}}.  (23)

Table 1: Operation counts and storage requirements for five algorithms, assuming separated end conditions (N = number of stages, n = dimension of each x_{i}, ℓ = number of left-hand end conditions, R = number of right-hand sides, p = number of partitions in first level of PSC)

Algorithm Operation Count Storage
LU (row pivoting) N[\frac{5}{3}n^{3} + 3ℓn^{2} + R(4n^{2} + 2ℓn)]  3Nn^{2}
DECOMP/SOLVE  N[\frac{2}{3}n^{3} + (4R + 5ℓ)n^{2} − 2nℓ^{2}]  2Nn^{2}
Structured QR  N[\frac{46}{3}n^{3} + 11Rn^{2}]  4Nn^{2}
Structured LU  N[\frac{23}{3}n^{3} + 8Rn^{2}]  4Nn^{2}
PSC  N[\frac{22}{3}n^{3} + 6Rn^{2}] + p[\frac{46}{3}n^{3} + 11Rn^{2}] \frac{5}{2}Nn^{2} + 4pn^{2}

The remaining equation, obtained from the last partition, is

A_{N}R_{N}^{(p)} [Φ_{N}^{(p)}s_{p} + ẑ_{N}^{(p)}] + B_{N}x_{N+1} = c_{N}.  (24)

Aggregating (22), (23), and (24), and defining s_{p+1} =^{△} x_{N+1}, we see that the reduced system has the form

Ã_{j}s_{j} + B̃_{j}s_{j+1} = c̃_{j},  j = 1, ⋯, p,  (25a)
M̃_{a}s_{1} + M_{b}s_{p+1} = d̃.  (25b)

This system has the same form as the original system (4), which immediately suggests that the partitioning/reduction process can be performed recursively. Such an approach is indeed desirable when the number of processors is large, making the reduced system (25) too large to solve on a single processor. We return to this point in §6.

Tables 1 and 2 compare storage requirements and operation counts for five algorithms, including the one discussed above (partitioned stabilized compactification, or PSC). In reporting the operation counts, lower-order (n^{2}) terms are ignored. We assume that, for each algorithm, the calculation of fundamental and particular solutions was carried out in separate phases. It is therefore necessary to store enough information to reconstruct the transformations that were used to factor the matrix associated with the recurrence (4). The statistics for the structured QR algorithm differ slightly from those reported in Wright [13], since we refer here to an improved variant of the technique which is based on Givens rotations rather than Householder transformations. (This variant will be described in a forthcoming report.) The statistics for PSC assume that the reduced system (25) is solved by using structured QR. We use this approach in our implementation in §6, since it can be implemented stably, in a fashion akin to cyclic reduction, on a binary tree of processors.

For problems with separated end conditions, all five algorithms are stable, while only structured QR and PSC are guaranteed to be stable when the end conditions are coupled (although LU and structured LU are almost always stable in these circumstances). The LU and DECOMP/SOLVE algorithms are not parallelizable unless n is large.

As well as halving the computational cost of structured QR, PSC has a clear advantage in storage requirement. The Householder vectors that define the orthogonal matrices Q_{i} and R_{i} can be largely stored in the lower triangle that is vacated when A_{i} and B_{i} are transformed to U_{i} and V_{i}. The only real fill-in is due to the fundamental solution matrices Φ_{i}. When N is significantly larger than p, the factorization process requires only about 25% more storage than the cost of storing the original system.

The algorithm PSC can be viewed as a particular factorization of the matrix corresponding to the recurrence (4). We state this result as a theorem.

Theorem 3.2 Algorithm PSC is equivalent to a particular solution scheme for the linear system Ax = c, Table 2: Operation counts and storage requirements for five algorithms, assuming coupled end conditions (N = number of stages, n = dimension of each x_{i}, R = number of right-hand sides, p = number of partitions in first level of PSC)

Algorithm Operation Count Storage
LU (row pivoting) N[\frac{23}{3}n^{3} + 8Rn^{2}]  4Nn^{2}
DECOMP/SOLVE  N[\frac{14}{3}n^{3} + 4Rn^{2}]  3Nn^{2}
Structured QR  N[\frac{46}{3}n^{3} + 11Rn^{2}]  4Nn^{2}
Structured LU  N[\frac{23}{3}n^{3} + 8Rn^{2}]  4Nn^{2}
PSC  N[\frac{22}{3}n^{3} + 6Rn^{2}] + p[\frac{46}{3}n^{3} + 11Rn^{2}] \frac{5}{2}Nn^{2} + 4pn^{2}

where

A = [ \begin{matrix}A_{1} & B_{1} \\ & A_{2} & B_{2} \\ & & ... & ... \\ & & & A_{N} & B_{N} \\ M_{a} & & & & M_{b} \end{matrix} ] , c = [ c_{1} ; c_{2} ; ⋮ ; c_{N} ; d ] .

In this scheme, A is factored in the form

P_{L}QARP_{R} = [ \begin{matrix}L̂_{1} & 0 \\ L̂_{2} & Ã \end{matrix} ] [ \begin{matrix}I & Û_{2} \\ 0 & I \end{matrix} ] = L̃Ũ,  (26)

where P_{L} and P_{R} are permutation matrices, Q and R are block diagonal matrices with orthogonal n × n blocks, and L̂_{1} is lower triangular. The block Ã in (26) has the form

Ã = [ \begin{matrix}Ã_{1} & B̃_{1} \\ & Ã_{2} & B̃_{2} \\ & & ... & ... \\ & & & Ã_{p} & B̃_{p} \\ M̃_{a} & & & & M_{b} \end{matrix} ]

(that is, the coeficient matrix associated with the recurrence (25)). The remaining blocks in (26) are defined in the proof.

Proof. We prove only the case of p > 1 (trivial modifications are required for p = 1). The matrices Q and R come from the initial orthogonal transformation phase. They are defined as

Q^{T} = diag(Q_{1}^{(1)}, Q_{2}^{(1)}, ⋯, Q_{k_{2}−1}^{(1)}, I, Q_{k_{2}+1}^{(2)}, ⋯, Q_{k_{2}−1}^{(2)}, I, ⋯, Q_{k_{p+1}−1}^{(p)}, I, I),
R = diag(R_{1}^{(1)}, R_{2}^{(1)}, ⋯, R_{k_{2}}^{(1)}, R_{k_{2}+1}^{(2)}, ⋯, R_{k_{3}}^{(2)}, R_{k_{3}+1}^{(3)}, ⋯, R_{k_{p+1}}^{(p)}, I).

Then QAR is

[ \begin{matrix}U_{1} & V_{1} \\ & ... & ... \\ & & U_{k_{2}−1} & V_{k_{2}−1} \\ & & & (A_{k_{2}}R_{k_{2}}^{(1)}) & (B_{k_{2}}R_{k_{2}+1}^{(2)}) \\ & & & & U_{k_{2}+1} & V_{k_{2}+1} \\ & & & & & ... & ... \\ & & & & & & A_{N}R_{N−1}^{(p)} & B_{N} \\ (M_{a}R_{1}^{(1)}) & & & & & & & M_{b} \end{matrix} ] .

Column pivoting is now performed on this matrix. The columns of QAR that contain the last (n−ℓ) columns of U_{k_{j}+1}, j = 1, ⋯, p, are shifted to the right of the matrix, as are the columns of QAR that contain the first ℓ columns of V_{k_{j}−1}, j = 2, ⋯, p + 1. Denoting the overall pivot matrix by P_{c}, we write the pivoted form as

QARP_{c} = [ \begin{matrix}T_{1} & & & & & S_{1} \\ Â_{1}^{(2)} & B̂_{1}^{(1)} & & & & Â_{1}^{(1)} & B̂_{1}^{(2)} \\ & T_{2} & & & & & S_{2} \\ & Â_{2}^{(2)} & B̂_{2}^{(1)} & & & & Â_{2}^{(1)} & B̂_{2}^{(2)} \\ & & T_{3} & & & & & S_{3} \\ & & Â_{3}^{(2)} & B̂_{3}^{(1)} & & & & Â_{3}^{(1)} & B̂_{3}^{(2)} \\ & & & ⋮ & & & & & ⋮ \\ & & & & T_{p} & & & & & S_{p} \\ 0 & 0 & ⋯ & 0 & Â_{p}^{(2)} & 0 & 0 & ⋯ & 0 & Â_{p}^{(1)} & B_{N} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 & 0 & M̂_{a}^{(2)} & 0 & ⋯ & 0 & 0 & M_{b} \end{matrix} ] .  (27)

The component blocks of this matrix can be described by using the following notation: N_{j} = (k_{j} −k_{j−1} −1)n denotes the dimension of the square matrix T_{j}, (.)_{.,1} denotes the first n − ℓ columns of an n × n matrix, and (.)_{.,2} denotes the last ℓ columns of an n × n matrix. Then

T_{j} = [ \begin{matrix}(U_{k_{j}+1})_{.,1} & V_{k_{j}+1} \\ & U_{k_{j}+2} & V_{k_{j}+2} \\ & & ... & ... \\ & & & U_{k_{j+1}−1} & (V_{k_{j+1}−1})_{.,2} \end{matrix} ] ∈ R l ^{N_{j}×N_{j}},

S_{j} = [ 0 (U_{k_{j}+1})_{.,2} ; 0 ; ⋮ ; 0 ; (V_{k_{j+1}−1})_{.,1} 0 ] ∈ R l ^{N_{j}×n},

B̂_{j}^{(1)} = [(B_{k_{j+1}}R_{k_{j+1}−1}^{(j+1)})_{.,1} 0 ⋯ 0] ∈ R l ^{n×N_{j+1}},  B̂_{j}^{(2)} = [0 (B_{k_{j+1}}R_{k_{j+1}−1}^{(j+1)})_{.,2}] ∈ R l ^{n×n},
Â_{j}^{(1)} = [(A_{k_{j+1}}R_{k_{j+1}}^{(j)})_{.,1} 0] ∈ R l ^{n×n},  Â_{j}^{(2)} = [0 ⋯ 0 (A_{k_{j+1}}R_{k_{j+1}}^{(j)})_{.,2}] ∈ R l ^{n×N_{j}},
M̂_{a}^{(1)} = [(M_{a}R_{1}^{(1)})_{.,1} 0 ⋯ 0] ∈ R l ^{n×N_{1}},  M̂_{a}^{(2)} = [0 (M_{a}R_{1}^{(1)})_{.,2}] ∈ R l ^{n×n}.

We can now perform a block LU factorization on the remaining matrix. By restating (16) and (17) as a system of linear equations, we find that

T_{j}Φ̂_{j} = −S_{j}, j = 1, ⋯, p,

where

Φ̂_{j} = [ (Φ_{k_{j}+1}^{(j)})_{1,.} ; Φ_{k_{j}+2}^{(j)} ; ⋮ ; Φ_{k_{j+1}−1}^{(j)} ; (Φ_{k_{j+1}}^{(j)})_{2,.} ] .

STABLE PARALLEL COMPACTIFICATION FOR ODEs 10 Therefore, (27) can be written as the following product:

QARP_{c} = [ \begin{matrix}L^{(1)} & 0 \\ 0 & I \end{matrix} ] [ \begin{matrix}U_{1}^{(1)} & U_{2}^{(1)} \\ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} & \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(1)} & B_{N} \\ M̂_{a}^{(2)} & 0 & ⋯ & 0 & M_{b} \end{matrix} \end{matrix} ] ,

where

L^{(1)} = [ \begin{matrix}T_{1} \\ & I \\ & & T_{2} \\ & & & I \\ & & & & . & . & . \\ & & & & & & & I \\ & & & & & & & & T_{p} \end{matrix} ]

(where each I in L^{(1)} is n × n),

U_{1}^{(1)} = [ I ; Â_{1}^{(2)} B̂_{1}^{(1)} ; I ; Â_{2}^{(2)} B̂_{2}^{(1)} ; I ; Â_{3}^{(2)} ; ... ; I ] , U_{2}^{(1)} = [ −Φ̂_{1} ; Â_{1}^{(1)} B̂_{1}^{(2)} ; −Φ̂_{2} ; Â_{2}^{(1)} B̂_{2}^{(2)} ; −Φ̂_{3} ; Â_{3}^{(1)} ; ... ; −Φ̂_{p} 0 ] .

Now

U_{1}^{(1)} = [ I ; Â_{1}^{(2)} I B̂_{1}^{(1)} ; I ; Â_{2}^{(2)} I B̂_{2}^{(1)} ; I ; Â_{3}^{(2)} ; ... ; I ] [ \begin{matrix}I \\ 0 & 0 \\ & I \\ & 0 & 0 \\ & & I \\ & & & . & . & . \\ & & & & & & I \end{matrix} ]
=^{△} L^{(2)}U_{1}^{(2)}.

A little manipulation using the definitions of the coefficient matrices for the reduced system (25), the definitions of Â_{j}^{(1)}, Â_{j}^{(2)}, B̂_{j}^{(1)}, and B̂_{j}^{(2)}, and the fundamental boundary conditions (15) yields

−Â_{j}^{(2)}Φ̂_{j} + Ã_{j} = Â_{j}^{(1)}, −B̂_{j}^{(1)}Φ̂_{j+1} + B̃_{j} = B̂_{j}^{(2)}.

Hence

U_{2}^{(1)} = [ \begin{matrix}I \\ Â_{1}^{(2)} & I & B̂_{1}^{(1)} \\ & & I \\ & & Â_{2}^{(2)} & I & B̂_{2}^{1} \\ & & & & I \\ & & & & & . & . & . \\ & & & & & & & & I \end{matrix} ] [ \begin{matrix}−Φ̂_{1} \\ Ã_{1} & B̃_{1} \\ & −Φ̂_{2} \\ & Ã_{2} & B̃_{2} \\ & & −Φ̂_{3} \\ & & & . & . & . \\ & & & & & & −Φ̂_{p} & 0 \end{matrix} ]

=^{△} L^{(2)}U_{2}^{(2)}.

Using the definitions of U_{2}^{(2)} and U_{1}^{(2)} and of M̂_{a}^{(1)}, M̂_{a}^{(2)}, M̂_{b}^{(1)}, and M̂_{b}^{(2)}, we have that

[ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} ] U_{1}^{(2)} = [ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} ]

(U_{1}^{(2)} has the effect of deleting a number of the zero columns) and

[ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} ] U_{2}^{(2)} + [ \begin{matrix}0 & ⋯ & 0 & Ã_{p} & B_{N} \\ M̃_{a} & 0 & ⋯ & 0 & M_{b} \end{matrix} ]

= [ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(1)} & B_{N} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 & M_{b} \end{matrix} ] .

Hence, by modifying the factorization (27), we find that

QARP_{c} = [ \begin{matrix}L^{(1)}L^{(2)} & 0 \\ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} & I \end{matrix} ] [ \begin{matrix}U_{1}^{(2)} & U_{2}^{(2)} \\ 0 & \begin{matrix}0 & ⋯ & 0 & Ã_{p} & B_{N} \\ M̃_{a} & 0 & ⋯ & 0 & M_{b} \end{matrix} \end{matrix} ] .

Now, let P_{r} be a permutation matrix that pivots all zero rows of U_{1}^{(2)} to the bottom:

P_{r}U_{1}^{(2)} = P_{r} [ \begin{matrix}I \\ 0 & 0 \\ & I \\ & 0 & 0 \\ & & I \\ & & & . & . & . \\ & & & & & & I \end{matrix} ] = [ \frac{I}{0} ] .

Then, augmenting P_{r} with a 2n × 2n identity matrix to obtain

P̃_{r} = [ \begin{matrix}P_{r} \\ & I \end{matrix} ] ,

we have that

P̃_{r}QARP_{c}
= P̃_{r} [ \begin{matrix}L^{(1)}L^{(2)} & 0 \\ \begin{matrix}0 & ⋯ & 0 & Â_{p}^{(2)} \\ M̂_{a}^{(1)} & 0 & ⋯ & 0 \end{matrix} & I \end{matrix} ] P̃_{r}^{T} P̃_{r} [ \frac{U_{1}^{(2)}  U_{2}^{(2)}}{\begin{matrix}0 & ⋯ & 0 & Ã_{p} & B_{N} \\ M̃_{a} & 0 & ⋯ & 0 & M_{b} \end{matrix}} ]

= [ \begin{matrix}L^{(3)} & 0 \\ L^{(4)} & I \end{matrix} ] [ \begin{matrix}I & U_{2}^{(3)} \\ 0 & Ã \end{matrix} ] ,

where

L^{(3)} = [ \begin{matrix}T_{1} \\ & T_{2} \\ & & . & . & . \\ & & & & & T_{p} \end{matrix} ] ,  L^{(4)} = [ Â_{1}^{(2)} B̂_{1}^{(1)} ; Â_{2}^{(2)} B̂_{2}^{(1)} ; ... ; Â_{p}^{(2)} ; M̂_{a}^{(1)}  0 ] ,

U_{2}^{(3)} = [ \begin{matrix}−Φ̂_{1} \\ & −Φ̂_{2} \\ & & . & . & . \\ & & & & & −Φ̂_{p} & 0 \end{matrix} ] ,  Ã = [ \begin{matrix}Ã_{1} & B̃_{1} \\ & Ã_{2} & B̃_{2} \\ & & ... & ... \\ & & & Ã_{p} & B_{N} \\ M̃_{a} & & & & M_{b} \end{matrix} ] .

Finally, since linear systems involving each matrix T_{j} (and hence L^{(3)}) can be solved by a triangular substitution process, it follows that there are permutation matrices P_{0} and P_{1} such that P_{0}L^{(3)}P_{1} is lower triangular. Hence

[ \begin{matrix}P_{0} \\ & I \end{matrix} ] P̃_{r}QARP_{c} [ \begin{matrix}P_{1} \\ & I \end{matrix} ]

= [ \begin{matrix}P_{0} \\ & I \end{matrix} ] [ \begin{matrix}L^{(3)} & 0 \\ L^{(4)} & I \end{matrix} ] [ \begin{matrix}P_{1} \\ & I \end{matrix} ] [ \begin{matrix}P_{1}^{T} \\ & I \end{matrix} ] [ \begin{matrix}I & U_{2}^{(3)} \\ 0 & Ã \end{matrix} ] [ \begin{matrix}P_{1} \\ & I \end{matrix} ]  (28)
= [ \begin{matrix}P_{0}L^{(3)}P_{1} & 0 \\ L^{(4)}P_{1} & I \end{matrix} ] [ \begin{matrix}I & P_{1}^{T} U_{2}^{(3)} \\ 0 & Ã \end{matrix} ] .

By making the obvious identifications

P_{L} = [ \begin{matrix}P_{0} \\ & I \end{matrix} ] P̃_{r}, P_{R} = P_{c} [ \begin{matrix}P_{1} \\ & I \end{matrix} ] ,

L̂_{1} = P_{0}L^{(3)}P_{1}, L̂_{2} = L^{(4)}P_{1}, Û_{2} = P_{1}^{T} U_{2}^{(3)},

we obtain (26).

To complete the proof, we note that the quantities in both factors on the right-hand side of (26) are quantities that are actually computed by Algorithm PSC. The remainder of the algorithm consists of

1. performing the orthogonal “preprocessing” of A with Q and R, and operating with Q on the right-hand side c;

2. permuting the rows and columns of QAR and the rows of Qc in a predetermined way;

3. performing an LU factorization with no further pivoting, stopping the elimination at stage (N +1−p)n;

4. doing a “forward substitution” with L̂_{2} to get the right-hand side of the reduced system;

5. continuing the forward substitution with L̂_{1} to get the particular solution ẑ_{i};

6. solving the reduced system (coefficient matrix Ã) to obtain the s_{j}; and

7. doing a back substitution, followed by orthogonal transformations, to recover all the x_{i}.

Stability of the algorithm can now be proved by using error analysis techniques from numerical linear algebra. We have the following result. Theorem 3.3 Suppose that PSC is used to solve (4) in finite-precision floating-point arithmetic, with unit roundoff u ≪ 1. Assume that

(i) the growth in the fundamental solution is not excessive; that is, there is a moderate constant γ_{1} such

that

max_{j=1,···,p} max_{i=k_{j−1}+1,···,k_{j}} ∥Φ_{i}^{(j)}∥_{2} ≤ γ_{1},  max_{j=1,···,p} ∥T_{j}^{−1}∥_{2} ≤ γ_{1};

(ii) the reduced system Ãs = c̃ is solved in a stable way; that is, there is a moderate constant γ_{2} such that the computed solution s̄ satisfies

(Ã + EÃ)s̄ = (c̃ + ec̃),  (29)

where

∥EÃ∥_{2} ≤ γ_{2}u∥Ã∥_{2},  ∥ec̃∥_{2} ≤ γ_{2}u∥c̃∥_{2}.  (30)

Then the computed solution x̄ of Ax = b satisfies

(A + E_{A})x̄ = (b + e_{b}),  (31)

where

∥E_{A}∥_{2} ≤ γ_{3}u∥A∥_{2},  ∥e_{b}∥_{2} ≤ γ_{4}(1 + ∥A∥_{2})u∥b∥_{2},  (32)

where γ_{3} = O(Nn^{5/2} + N^{3}np^{−1/2}γ_{1}^{2}γ_{2}) and γ_{4} = O(n^{2}p^{1/2}Nγ_{1}γ_{2}).

Proof. We assume throughout the proof that ∥.∥ denotes the Euclidean norm ∥.∥_{2}. Our result depends on standard analysis for Householder QR factorizations (Lawson and Hanson, [8, pp. 86–89]) and LU factorizations (Golub and Van Loan [6, §3.3]). The key to the stability argument is the fact that element growth in the L and U factors of P_{L}QARP_{R} is bounded, because of (i). We stress that our assumption (i) is reasonable: when (4) is well conditioned, a dichotomy exists, and so our use of decoupling will ensure that (i) holds provided that the starting matrices R_{k_{j}+1}^{(j)} for each partition are chosen appropriately.

First, we take account of the errors arising from the initial orthogonal transformation. By the argument of Lawson and Hanson [8, p. 87], we have that if C ∈ R l ^{n×n} is a general matrix and Q ∈ R l ^{n×n} is a product of n Householder transformations, then

comp(QC) = Q(C + H),

where comp(.) denotes the computed value of its argument, taking roundoff error into account, and

∥H∥ ≤ (3n + 40)n^{3/2}u∥C∥ + O(u^{2}).

For the case of a vector c ∈ R l ^{n},

comp(Qc) = Q(c + h),

where

∥h∥ ≤ (3n + 40)nu∥c∥ + O(u^{2}).

Applying these results to the initial orthogonal transformation of A, we find that

U_{i} = comp(Q_{i}^{T} comp(A_{i}R_{i}))
= Q_{i}^{T} [(A_{i} + H_{i}^{1})R_{i} + H_{i}^{2}]
= Q_{i}^{T} [A_{i} + H_{i}^{3}]R_{i},

where

∥H_{i}^{3}∥ ≤ ∥H_{i}^{1}∥ + ∥H_{i}^{2}∥ ≤ 2(3n + 40)n^{3/2}∥A_{i}∥u + O(u^{2}).

Similarly,

V_{i} = Q_{i}^{T} [B_{i} + H_{i}^{4}]R_{i+1},

where

∥H_{i}^{4}∥ ≤ 2(3n + 40)n^{3/2}∥B_{i}∥u + O(u^{2}).

Similar analysis can be applied to comp(M_{a}R_{1}^{(1)}) and to the block rows k_{j} that are not multiplied from the left by an orthogonal transformation. Since the pivoting operations do not incur any roundoff error,

comp(P_{L}QARP_{R}) = P_{L}Q(A + H̃)RP_{R},  (33)

STABLE PARALLEL COMPACTIFICATION FOR ODEs 14 where

∥H̃∥ ≤ (N + 1) [sup_{i} ∥H_{i}^{3}∥ + sup_{i} ∥H_{i}^{4}∥]
≤ 4(N + 1)n^{3/2}(3n + 40)∥A∥u + O(u^{2}).

Similarly,

comp(Qb) = Q(b + h),

where

∥h∥ ≤ (N + 1)^{1/2}(3n + 40)n∥b∥u + O(u^{2}).

We now examine the LU factorization of comp(P_{L}QARP_{R}). Although the Gaussian elimination process is terminated prematurely, we can use the proof of Theorem 3.3.1 of Golub and Van Loan [6] to show that

L̃Ũ = comp(P_{L}QARP_{R}) + Ĥ,  (34)

where

|Ĥ| ≤ 3(N + 1)nu (|comp(P_{L}QARP_{R})| + |L̃||Ũ|) + O(u^{2}).

Hence,

∥Ĥ∥_{2} ≤ 3(N + 1)nu [∥A∥ + ∥L̃∥∥Ũ∥] + O(u^{2}).  (35)

Now

∥L̃∥ ≤ ∥L̂_{1}∥ + ∥L̂_{2}∥ + ∥Ã∥ = ∥L^{(3)}∥ + ∥L^{(4)}∥ + ∥Ã∥,  (36)

and so, by using the definitions that appear during the proof of Theorem 3.2,

∥L^{(3)}∥ ≤ ∥comp(QAR)∥ ≤ ∥A∥ + O(u),
∥L^{(4)}∥ ≤ 2(p + 1)^{1/2} max (max_{i}(∥A_{i}R_{i}∥, ∥B_{i}R_{i+1}∥), ∥M_{a}R_{1}∥) + O(u)
≤ 2(p + 1)^{1/2}∥A∥ + O(u),
∥Ã∥ ≤ 2(p + 1)^{1/2} max (max_{j}(∥Ã_{j}∥, ∥B̃_{j}∥), ∥M̃_{a}∥, ∥M_{b}∥) .

In the last of these inequalities,

∥Ã_{j}∥ ≤ ∥A_{k_{j}+1}∥(1 + ∥Φ̂_{j}∥ + O(u)) ≤ ∥A∥(1 + (N/p)γ_{1} + O(u)).

Similar inequalities can be derived for ∥B̃_{j}∥ and ∥M̃_{a}∥, so

∥Ã∥ ≤ 2(p + 1)^{1/2}(1 + (N/p)γ_{1})∥A∥ + O(u).

For the upper triangular factor, we have that

∥Ũ∥ ≤ 2 + ∥Û_{2}∥ = 2 + ∥U_{2}^{(3)}∥

≤ 2 + max_{j} ∥Φ̂_{j}∥ + O(u) ≤ 2 + (N/p)γ_{1} + O(u).  (37)

By using some elementary inequalities, we then have from (35), (36), and (37) that

∥Ĥ∥ ≤ 6(N + 1)(p + 1)^{1/2}n(3 + (N/p)γ_{1})^{2}∥A∥u + O(u).

Consider now the forward and back substitution process involving L̃ and Ũ, that is,

L̃Ũy = comp(P_{L}Qb) = P_{L}Q(b + h̃).  (38)

Using (29), (30), and a simple backwards error argument, we can show that the computed solution ȳ of (38) will satisfy

[ \begin{matrix}L̂_{1} + Ĥ_{1} & 0 \\ L̂_{2} + Ĥ_{2} & Ã + EÃ \end{matrix} ] [ \begin{matrix}I & Û_{2} + Ĥ_{3} \\ 0 & I \end{matrix} ] [ ȳ_{1} ; ȳ_{2} ] = [ \begin{matrix}Q̄_{1}(b + h̃) \\ Q̄_{2}(b + h̃) + ẽ \end{matrix} ]  (39)
⟹ (L̃ + EL̃) (Ũ + EŨ) ȳ = P_{L}Q(b + h̃ + eb̃),

where Q̄_{1} consists of the first (N − p)n rows of P_{L}Q, Q̄_{2} contains the last (p + 1)n rows, and EL̃, EŨ, and eb̃ are defined in obvious ways. Now (29) and (30) imply that

∥ẽ∥ ≤ γ_{2}u ∥Q_{1}(b + h̃) − (L̂_{2} + Ĥ_{2})(L̂_{1} + Ĥ_{1})^{−1}Q_{2}(b + h̃)∥
≤ γ_{2}u (1 + ∥L̂_{2}∥∥L̂_{1}^{−1}∥ + O(u)) (∥b∥ + O(u))
≤ γ̄_{4}u∥A∥∥b∥ + O(u^{2}),

where γ̄_{4} = O(γ_{1}γ_{2}p^{1/2}), and

∥EÃ∥ ≤ γ_{2}u∥Ã∥ ≤ 2γ_{2}(p + 1)^{1/2}(1 + (N/p)γ_{1})u∥A∥ + O(u^{2}).

Following Golub and Van Loan [6, p. 106], we find that the remaining error terms in (39) satisfy

∥Ĥ_{1}∥ ≤ nu∥L̂_{1}∥ ≤ nu∥A∥ + O(u^{2}),
∥Ĥ_{2}∥ ≤ nu∥L̂_{2}∥ ≤ 2nu(p + 1)^{1/2}∥A∥ + O(u^{2}),
∥Ĥ_{3}∥ ≤ nu∥Û_{2}∥ ≤ 2nu(N/p)γ_{1}.

Hence

∥EL̃∥ ≤ ∥Ĥ_{1}∥ + ∥Ĥ_{2}∥ + ∥EÃ∥ ≤ γ̄_{3}u∥A∥ + O(u^{2}),
∥EŨ∥ ≤ 2nNγ_{1}u,
∥eb̃∥ ≤ γ̄_{4}u∥A∥∥b∥,

where γ̄_{3} = O(np^{−1/2}Nγ_{1}γ_{2}). From (33), (34), and (39), we have that

(L̃Ũ + EL̃Ũ + L̃EŨ + O(u^{2})) ȳ = P_{L}Q(b + h̃ + eb̃)
⟹ (P_{L}QARP_{R} + H̄) ȳ = P_{L}Q(b + h̃ + eb̃),

where

H̄ = P_{L}QH̃RP_{R} + Ĥ + EL̃Ũ + L̃EŨ + O(u^{2}).

Hence

∥H̄∥ ≤ ∥H̃∥ + ∥Ĥ∥ + ∥EL̃∥∥Ũ∥ + ∥L̃∥∥EŨ∥ ≤ γ_{3},

where γ_{3} is as defined in the statement of the theorem. Also,

∥h̃ + eb̃∥ ≤ γ_{4}∥b∥.

The result follows by making the identifications

E_{A} = Q^{T} P_{L}^{T} H̄P_{R}^{T} R^{T} ,  e_{b} = h̃ + eb̃,  ȳ = RP_{R}ȳ.

If the algorithm of this section is applied recursively (that is, if the reduced system is itself solved on multiple processors by using the same technique), the theorem indicates the amount of deterioration in accuracy that can be expected in moving between levels of the recursion.

4 A No-fill-in Variant

In the preceding section, we showed how the fundamental and particular solutions for the problem (4) could be calculated by using the blocks U_{i} and V_{i} from the transformed coefficient matrix. As described there, the algorithm assumed that the data for the initial orthogonal reduction (that is, the U_{i}, V_{i}, Q_{i}, and R_{i} matrices) is not overwritten in storage. This assumption means that new storage must be used to store the fundamental solution blocks Φ_{i}, resulting in a fill-in of approximately 25%.

In this section, we outline a scheme in which the fill-in is reduced by storing the components of Φ_{i} in storage formerly occupied by U_{i} and V_{i}. At the same time, we wish to allow for the possibility that the particular solution is calculated at some later time than the fundamental solution, that is, the right-hand sides c_{i}, i = 1, ⋯, N and d are not known at the time the fundamental solution is calculated. This situation may arise when a “chord method” approach is applied to a nonlinear version of (4). In this method, approximate Newton iterations are calculated by using the same Jacobian information for a number of successive iterations, while the right-hand side changes from iteration to iteration. In the preceding section, particular solutions are calculated by using (19) and (20). Hence, if we plan to overwrite components of U_{i} and V_{i}, we need to devise new formulae for finding the ẑ_{i}s.

To satisfy both requirements of the preceding paragraph, we need to define a “general solution” Ψ_{i}^{(j)}, i = k_{j} + 1, ⋯, k_{j+1} on partition j in such a way that each Ψ_{i}^{(j)} ∈ R l ^{n×n} is nonsingular. Moreover, Ψ_{i}^{(j)} will be identical to the fundamental solution Φ_{i}^{(j)} whenever all U_{i} and V_{i} in partition j are nonsingular. Computation of the sequence Ψ_{i} involves a simple modification of the forward sweep/reverse sweep process (16) and (17). The only enhancement is that when some (U_{i})_{22} or (V_{i})_{11} are singular, diagonal terms are added to (U_{i})_{22}(Ψ_{i})_{22} or (V_{i})_{11}(Ψ_{i+1})_{11} to prevent singularity in the next general solution matrix in the sequence. Note that we are taking advantage here of the fact that U_{i} and V_{i} (and also Φ_{i} and Ψ_{i}) are upper triangular. This means that singularity manifests itself as zero elements on the diagonal and is easily remedied by replacing these zeros by, say, ones, which is what we do.

For i = k_{j} + 1, ⋯, k_{j+1} − 1, we define

Z_{i}^{F} = [e_{l}^{ℓ}]_{l∈A_{i}^{F}} ,

where e_{l}^{ℓ} is the R l ^{ℓ} unit vector, with zeros everywhere except for a 1 in position l, and

A_{i}^{F} = {l | l-th diagonal element of (U_{i})_{22} is zero}.

Note that Z_{i}^{F} (Z_{i}^{F})^{T} is an ℓ × ℓ matrix that is zero everywhere except for ones in the diagonal positions in which (U_{i})_{22} has zeros. Similarly, define

Z_{i}^{R} = [e_{l}^{n−ℓ}]_{l∈A_{i}^{R}},

where e_{l}^{n−ℓ} is the R l ^{n−ℓ} unit vector with a 1 in the l position, and

A_{i}^{R} = {l | l-th diagonal element of (V_{i})_{11} is zero}.

The boundary conditions for Ψ_{i} coincide with (15); that is,

(Ψ_{k_{j}+1})_{2,.} = [0 | I],  (Ψ_{k_{j+1}})_{1,.} = [I | 0].  (40)

The sweeps are defined as follows:

for i = k_{j} + 1, ⋯, k_{j+1} − 1

(Ψ_{i+1})_{2,.} = −(V_{i})_{22}^{−1} \{[0 | (U_{i})_{22}(Ψ_{i})_{22}] + Z_{i}^{F} [ \begin{matrix}0 \\ Z_{i}^{F} \end{matrix} ]^{T} \};  (41)

for i = k_{j+1} − 1, ⋯, k_{j} + 1

(Ψ_{i})_{1,.} =
−(U_{i})_{11}^{−1} \{(V_{i})_{11}(Ψ_{i+1})_{1,.} + (V_{i})_{12}(Ψ_{i+1})_{2,.} + (U_{i})_{12}(Ψ_{i})_{2,.} + Z_{i}^{R} [ \begin{matrix}Z_{i}^{R} \\ 0 \end{matrix} ]^{T} \} .  (42)

Note that the extra terms in (41) and (42) do not disturb the upper triangularity of each Ψ_{i}.

By making use of the relationship between the fundamental and general solution, we can store the general solution in a compressed format. We show how this is done for the (2, 2) block of Ψ_{i}; for the other two nonzero blocks the technique is similar.

First, we show that we can write

(Ψ_{i})_{22} = (Φ_{i})_{22} + W̄_{i}^{F} (Z̄_{i}^{F})^{T} ,  (43)

where Z̄_{i}^{F} contains all e_{l}^{ℓ} such that the l-th diagonal of (U_{k})_{22} is zero for some, and possibly more than one, k = k_{j} + 1, ⋯, i − 1. For i = k_{j} + 1, W̄_{i}^{F} and Z̄_{i}^{F} are null. Assuming that (43) is true for k = k_{j} + 1, …, i, we have from (41) that

(Ψ_{i+1})_{22} = −(V_{i})_{22}^{−1} \{(U_{i})_{22} [(Φ_{i})_{22} + W̄_{i}^{F} (Z̄_{i}^{F})^{T} ] + Z_{i}^{F} [ \begin{matrix}0 \\ Z_{i}^{F} \end{matrix} ]^{T} \}
= (Φ_{i+1})_{22} − (V_{i})_{22}^{−1}(U_{i})_{22}W̄_{i}^{F} (Z̄_{i}^{F})^{T} − (V_{i})_{22}^{−1}Z_{i}^{F} (Z_{i}^{F})^{T} .

We obtain Z̄_{i+1}^{F} by merging the columns of Z̄_{i}^{F} and Z_{i}^{F}, while W̄_{i+1}^{F} is obtained by merging −(V_{i})_{22}^{−1}(U_{i})_{22}W̄_{i}^{F} and −(V_{i})_{22}^{−1}Z_{i}^{F}. We “merge” rather than simply append Z_{i}^{F} to Z̄_{i}^{F} since there is no need for Z̄_{i+1}^{F} to have two copies of the same column. (The same effect can be obtained by adding the two corresponding columns of W̄_{i+1}^{F}.) The corresponding formula to (43) for the (1, 1) and (1, 2) blocks is

(Ψ_{i})_{1,.} = (Φ_{i})_{1,.} + W̄_{i}^{R}(Z̄_{i}^{R})^{T} .  (44)

In each case, the substantive additional storage requirements are for the matrices W̄_{i}^{F} and W̄_{i}^{R} (the Z̄_{i} matrices can be stored in a few integer locations.) In the worst case, this will require the same amount of storage as the general solution itself, but we usually expect it to be much less. For example, if the (ℓ, ℓ) element of one or more of the (U_{i})_{22} matrices is zero, then we need about (k_{j+1} − k_{j} − 1)n locations to store the W̄_{i}, or about 2/n of the space required by the entire general solution.

As we noted earlier, a particular solution of the recurrence (13) can be calculated by performing the forward and backward sweeps (19) and (20). Since we would like to overwrite some components of the U_{i} and V_{i} matrices by components of Φ_{i}, we now describe an alternative method for calculating the ẑ_{i} which makes use of the Ψ_{i} but only of selected components of U_{i} and V_{i}.

The boundary conditions (18) are used for ẑ_{i}, i = k_{j} + 1, ⋯, k_{j+1}, as before. Defining the change of variables

ẑ_{i} = Ψ_{i}v_{i},  (45)

and using the boundary conditions (40), we obtain

(v_{k_{j}+1})_{2,.} = 0,  (v_{k_{j+1}})_{1,.} = 0.  (46)

We now substitute in (13) and define ĉ_{i} = Q_{i}^{T} c_{i} to obtain

U_{i}Ψ_{i}v_{i} + V_{i}Ψ_{i+1}v_{i+1} = ĉ_{i}.  (47)

Note from (41) and (42) that

− U_{i}Ψ_{i} = V_{i}Ψ_{i+1} + [ \begin{matrix}Z_{i}^{R}(Z_{i}^{R})^{T} & 0 \\ 0 & Z_{i}^{F} (Z_{i}^{F})^{T} \end{matrix} ] .  (48)

By isolating the last ℓ rows of (47), and by considering the (2, 2) block of (48), we obtain

−[(V_{i})_{22}(Ψ_{i+1})_{22} + Z_{i}^{F} (Z_{i}^{F})^{T}](v_{i})_{2} + (V_{i})_{22}(Ψ_{i+1})_{22}(v_{i+1})_{2} = (ĉ_{i})_{2}
⇒ (v_{i+1})_{2} = (v_{i})_{2} + [(V_{i})_{22}(Ψ_{i+1})_{22}]^{−1}[(ĉ_{i})_{2} + Z_{i}^{F} (Z_{i}^{F})^{T} (v_{i})_{2}].  (49)

Here, (.)_{2} denotes the last ℓ rows of a vector in R l ^{n}. By isolating the first n − ℓ rows of the expressions (47) and (48), we find that

(v_{i})_{1} = (v_{i+1})_{1} + [(U_{i})_{11}(Ψ_{i})_{11}]^{−1} {Z_{i}^{R}(Z_{i}^{R})^{T} (v_{i+1})_{1} + (ĉ_{i})_{1}
+[(U_{i})_{11}(Ψ_{i})_{12} + (U_{i})_{12}(Ψ_{i})_{22}][(v_{i+1})_{2} − (v_{i})_{2}]} .  (50)

We conclude that v_{i}, i = k_{j} + 1, ⋯, k_{j+1}, can be found by doing a forward sweep using (49) followed by a reverse sweep using (50). Note that we need to solve linear systems with coefficient matrices (V_{i})_{22}, (U_{i})_{11}, (Ψ_{i})_{11}, and (Ψ_{i})_{22} in order to obtain ẑ_{i} from (45) and (49),(50). If we assume that (4) is well conditioned and that the partitioning is done correctly, these matrices are invertible.

The blocks (U_{l})_{22}, (V_{l})_{11}, and (V_{l})_{12} are not used in (49) and (50). We can therefore overwrite these blocks in storage to avoid fill-in. During the forward sweep (16), (U_{i})_{22} can be progressively overwritten by (V_{i})_{22} in memory, and (Φ_{i+1})_{22} can be stored in the space vacated by (V_{i})_{22}. During the reverse sweep (17), (Φ_{i})_{11} and (Φ_{i})_{12} can progressively overwrite (V_{i})_{11} and (V_{i})_{12}. The matrices W̄_{i}^{F} and W̄_{i}^{R} which are needed to recover Ψ_{i} from Φ_{i} will need to be stored in new locations.

5 Extension to Problems with Multipoint Conditions or Parameters

We now show how the algorithm of §3 can be modified to handle problems of the forms (5) and (6). As mentioned in §2, the fundamental solution modes in these problems exhibits not dichotomy but, in general, polychotomy and skew-polychotomy. The practical consequence for the Algorithm PSC is that the value of ℓ (the number of “decreasing” fundamental modes) is no longer constant across all partitions; in fact, it may increase or decrease repeatedly within each partition. It is even possible for the submatrix [A_{i} B_{i}] to be rank deficient. In order to adapt PSC to these circumstances, we need to be able to recognize when ℓ has changed, and to modify the compactification strategy accordingly.

A change in ℓ is recognized by periodically examining the diagonal elements of V_{i}^{−1}U_{i} during the initial orthogonal factorization process (12). As mentioned earlier, these are (U_{i})_{ll}/(V_{i})_{ll}, l = 1, …, n. If the number of diagonals that are less than 1 is either greater than or less than the current value of ℓ for a few successive stages, then we deem ℓ to have changed. The most general way of handling such a dichotomy change is simply to break off the present partition at the current stage point, and start a new one. Specifically, suppose that at stage i during the orthogonal preprocessing of partition j, we decide that ℓ has changed. We then set k_{j+} = i + 1 and terminate the preprocessing of the current partition after calculating R_{k_{j+}}^{(j)}.

Next, we skip a row and start a new partition by choosing R_{k_{j+}+1}^{(j+)} as described in §3. The skipped row must now be added to the reduced system. If s_{j+} ∈ R l ^{n} is the “reduced” variable for the new partition, the extra equation is

A_{k_{j+}}R_{k_{j+}}^{(j)} [Φ_{k_{j+}}^{(j)}s_{j} + ẑ_{k_{j+}}^{(j)}] + B_{k_{j+}}R_{k_{j+}+1}^{(j+)} [Φ_{k_{j+}+1}^{(j+)}s_{j+} + ẑ_{k_{j+}+1}^{(j+)}] = c_{k_{j+}}.

To avoid creating a new partition in the case in which ℓ increases, we can pick up an extra component during the forward sweep calculation of the “decreasing” part of Φ_{i}, at the point at which the dichotomy change occurs. Correspondingly, a row is dropped from the first part of Φ_{i} during the backward sweep. This is essentially the strategy used by de Hoog and Mattheij [4, §4].

A circumstance that causes more immediate failure of the algorithm occurs when the submatrix [A_{i} B_{i}] fails to have full rank for some i. When this occurs, [U_{i} V_{i}] is also rank deficient, and it is easy to show that for at least one index l = 1, ⋯, n, the l-th diagonals of U_{i} and V_{i} are both zero. Hence at least one of (U_{i})_{11} and (V_{i})_{22} are singular, so either the forward sweep (16) or the reverse sweep (17) will break down at stage i. Since this difficulty can be detected during the initial orthogonal factorization, the fix is the same as for a dichotomy change — we create a partition break at stage i.

We assume from this point on that the number of partitions p and the separator indices k_{j}, j = 1, ⋯, p, have been altered where necessary during the partitioned compactification, to reflect the number of new dichotomy switches that were encountered.

Construction of the reduced system is slightly different for problems of the form (5) than it is for problems of the form (4). Rather than (21), each z_{i} is now expressible as

z_{i} = Φ_{i}^{(j)}s_{j} + Λ_{i}^{(j)}λ + ẑ_{i}^{(j)},  x_{i} = R_{i}^{(j)}z_{i},  i = k_{j} + 1, ⋯, k_{j+1}.  (51)

Since

U_{i}Φ_{i}^{(j)} + V_{i}Φ_{i+1}^{(j)} = 0,  i = k_{j} + 1, ⋯, k_{j+1},
U_{i}ẑ_{i}^{(j)} + V_{i}ẑ_{i+1}^{(j)} = Q_{i}^{T} c_{i},  i = k_{j} + 1, ⋯, k_{j+1},
U_{i}z_{i} + V_{i}z_{i+1} + Q_{i}^{T} C_{i}λ = Q_{i}^{T} c_{i},  i = k_{j} + 1, ⋯, k_{j+1},

we have by substitution in (51) that

[U_{i}Λ_{i}^{(j)} + V_{i}Λ_{i+1}^{(j)} + Q_{i}^{T} C_{i}] λ = 0, i = k_{j} + 1, ⋯, k_{j+1}.

We choose Λ_{i}^{(j)}, i = k_{j} + 1, ⋯, k_{j+1}, to satisfy the recurrence suggested by this formula, namely,

U_{i}Λ_{i}^{(j)} + V_{i}Λ_{i+1}^{(j)} = −Q_{i}^{T} C_{i},  i = k_{j} + 1, ⋯, k_{j+1}.  (52)

Clearly, this recurrence has the same form as the one that is solved for the particular solution ẑ_{i}^{(j)} (except that it has m columns instead of one), and we can solve it in exactly the same way. The choice of boundary conditions is also the same as for ẑ_{i}^{(j)}. When no dichotomy switch is encountered within a partition, we set

(Λ_{k_{j}+1}^{(j)})_{n−ℓ+1:n,.} = 0,  (Λ_{k_{j+1}}^{(j)})_{1:n−ℓ,.} = 0.  (53)

In deriving the scheme (52) and (53), we are essentially treating the parameter term as a forcing term. We can now use (51) to construct the reduced system. For the side conditions, we have

∑_{i=1}^{N+1} M_{i}x_{i} + Nλ = d

⇒ ∑_{j=1}^{p} ∑_{i=k_{j}+1}^{k_{j+1}} M_{i}R_{i}^{(j)} [Φ_{i}^{(j)}s_{j} + Λ_{i}^{(j)}λ + ẑ_{i}^{(j)}] + Nλ = d  (54)

⇒ ∑_{j=1}^{p} M̃_{j}s_{j} + Ñλ = d̃,

where

M̃_{j} = ∑_{i=k_{j}+1}^{k_{j+1}} M_{i}R_{i}^{(j)}Φ_{i}^{(j)},  Ñ = N + ∑_{j=1}^{p} ∑_{i=k_{j}+1}^{k_{j+1}} M_{i}R_{i}^{(j)}Λ_{i}^{(j)},
d̃= d − ∑_{j=1}^{p} ∑_{i=k_{j}+1}^{k_{j+1}} M_{i}R_{i}^{(j)}ẑ_{i}^{(j)}.

For the remaining equations, we obtain

Ã_{j}s_{j} + B̃_{j}s_{j+1} + C̃_{j}λ = c̃_{j},

where Ã_{j}, B̃_{j} and c̃_{j} are as defined in (25a), and

C̃_{j} = C_{k_{j+1}} + A_{k_{j+1}}R_{k_{j+1}}^{(j)}Λ_{k_{j+1}}^{(j)} + B_{k_{j+1}}R_{k_{j+1}+1}^{(j+1)}Λ_{k_{j+1}+1}^{(j+1)}.

Again, recursive application of the compactification to the reduced system is possible, but only up to a point. We would expect to have an intrinsic lower bound on the number of stages in the smallest possible reduced system, namely, the total number of dichotomy switches. In other words, the smallest reduced system is one in which a dichotomy switch occurs at each stage.

For problems of the form (6), substitution like that described in (54) can be performed in the least squares objective function. We obtain

∑_{i=1}^{N+1} ∥M_{i}x_{i} + N_{i}λ − d_{i}∥^{2} + ∥N_{0}λ − d_{0}∥^{2}
= ∑_{j=1}^{p} ∥M̂_{j}s_{j} + N̂_{j}λ − d̂_{j}∥^{2} + ∥N_{0}λ − d_{0}∥^{2},  (55)

where

M̂_{j} = [M_{i}R_{i}^{(j)}Φ_{i}^{(j)}]_{k_{j+1}}^{i=k_{j}+1} ,
N̂_{j} = [N_{i} + M_{i}R_{i}^{(j)}Λ_{i}^{(j)}]_{k_{j+1}}^{i=k_{j}+1} ,
d̂_{j} = [d_{i} − M_{i}R_{i}^{(j)}ẑ_{i}^{(j)}]_{k_{j+1}}^{i=k_{j}+1} .

6 Numerical Results

We implemented the PSC and SQR algorithms on the Intel Touchstone Delta at the Concurrent Supercomputer Consortium. This machine has 513 computational nodes, each with an Intel i860 microprocessor and 16 Megabytes of memory. The nodes are arranged in a two-dimensional mesh configuration. Underlying the mesh is a high-speed bus which allows data to be ported between any two nodes in a fashion that is transparent to the user. Provided the total amount of data being transferred is not too great, the time to send a message across this bus is not strongly dependent on the physical locations of the source and destination nodes in the mesh. We chose this architecture since it seems typical of the new generation of massively parallel machines. The forthcoming Intel Paragon will be a commercial version of the Touchstone Delta, and the Thinking Machines CM-5 is architecturally similar.

Our codes handle two classes of problem — two-point boundary value problems (1) and (4), and two-point problems with parameters ((2) and (5) with τ_{1} = a and τ_{2} = b). As discussed earlier, the SQR algorithm is similar to the one described in [13], except that Givens rotations are used in addition to Householder transformations during the QR factorizations. The extra columns in the shooting matrix that are due to the presence of parameters are handled as described in Wright [12, §5]. Each processor compresses its “slice” of the shooting matrix into a single block row, leaving a reduced system with (p + 1)n + m rows and columns. This system is solved by using a “cyclic reduction” variant of SQR, which we now briefly describe. Assume that the processors are numbered 0, 1, 2, ⋯, p − 1, where p = 2^{d} for some integer d ≥ 0. At the first level of cyclic reduction, the odd-numbered processors 2i + 1, i = 0, 1, ⋯, p/2 − 1 pass their piece of the reduced system (consisting of a single block row) to the neighboring even-numbered processor 2i. Each processor 2i then compresses this block row with its own block row to produce a single block row. At the end of this first level, the size of the reduced system has been approximately halved — it now has dimension (p/2 + 1)n + m. At the second level, the processors 4i + 2 pass their data to processors 4i. After d levels, processor 0 is left with a reduced system of dimension 2n + m. It solves this system to produce x_{1}, x_{N+1}, and λ. We then backtrack along the binary tree traversed during the process of compression to recover all the intermediate solution components.

The implementation of PSC is similar. Each processor is assigned an equal-sized slice of the original shooting matrix to which it applies the stabilized compactification algorithm described in §§3 and 5. At the end of this stage, each processor contains one or, in the case of dichotomy changes, more than one block row of the reduced system. If a processor contains more than one block row, orthogonal compression like that used by SQR is applied to obtain a single block row. The cyclic reduction SQR compression scheme just described is now applied to solve the reduced system that remains.

In our tests, we form discrete problems by applying the “box scheme” to our continuous problems. We seek x_{i} ∈ R l ^{n}, i = 1, ⋯, N + 1 such that x_{i} ≈ x(t_{i}), where t_{i} = a + (i − 1)h and h = (b − a)/N. The parametrized ODE (2) is approximated by

(x_{i+1} − x_{i})/h = A(t_{i+1/2})(x_{i} + x_{i+1})/2 + C(t_{i+1/2})λ + c(t_{i+1/2}).

Our first test problem has a change of dichotomy near the point t = 1/3.

Example 1 n = 3, m = 1.

ẋ = Q(t) [ \begin{matrix}20 \\ & 10(t − 1/3) \\ & & −20 \end{matrix} ] R(t)x + C(t)λ + c(t),  t ∈ [0, 1],  (56)

where, using the shorthand c = cos t, s = sin t,

Q(t) = [ \begin{matrix}c & s & 0 \\ −s & c & 0 \\ c − s & c + s & 1 \end{matrix} ] ,  R(t) = [ \begin{matrix}1 & & 1 \\ & 1 & 1 \\ & & 1 \end{matrix} ] ,  C(t) = [ \begin{matrix}3 \\ 0 \\ 5t \end{matrix} ] .

The vector function c(t) is chosen so that x(t) = e^{t}(1, 1, 1)^{T} and λ = 1 is a solution of the parametrized ODE. The boundary conditions are

[ \begin{matrix}1 \\ & 1 \\ & & 1 \\ 2 & 3 & 4 \end{matrix} ] x(0) + [ \begin{matrix}& & 1 \\ & 1 \\ 1 \\ −2 & −3 & −4 \end{matrix} ] x(1) + [ \begin{matrix}1 \\ 0 \\ −1 \\ 0 \end{matrix} ] λ = [ \begin{matrix}2 + e \\ 1 + e \\ e \\ 9 − 9e \end{matrix} ] .

PSC is able to detect the change in dichotomy near t = 1/3. On a single processor with N = 1000, a change from two decreasing modes to one decreasing mode is reported at t = .406. The “rightward shift” is due to the fact that our heuristic for detecting dichotomy changes is rather conservative; it reports a switch only when the behavior is consistently different over a significant number of consecutive intervals. When we use four processors, each with 1000 intervals (a total of N = 4000), the situation is a little more complicated. The algorithm for calculating the fundamental solution detects a number of points at which a dichotomy switch appears to occur and reports them as such. However, in all but one case, a restart at that point (as described in §5) indicates that the number of increasing and decreasing modes has not changed. The exception is, of course, the point at which the one true dichotomy switch occurs. The spurious switches do not affect the stability of the algorithm and have only a marginal effect on the computation time, since they add just a single row to the reduced system. The locations of the “breaks” are summarized in Table 3.

Tables 4 and 5 show timings for SQR and PSC, respectively. In both cases, the time for initial reduction of each partition remains essentially constant and is somewhat less for PSC than for SQR. The time to solve the reduced system by the cyclic reduction algorithm tends to increase as the depth of the tree increases, though not in a smooth way. Nevertheless, solution of the reduced system is such a small part of the overall computation that near-perfect speedup is attained.

Our second example arises from transport theory. Consider the single-group one-dimensional transport equation as defined by Jin and Levermore [7, Example 2]:

μ∂_{x}Φ(x, μ) + σ_{T} (x)Φ(x, μ) = \frac{1}{2}w(x)σ_{T} (x) ∫_{−1}^{1} Φ(x, μ′) dμ′,  x ∈ (0, L),  (57)

Table 3: Example 1, four processors, N = 4000. Restart points for stabilized compactification (excluding the breakpoints due to partitioning into four subintervals) and number of decreasing modes on the subinterval between this breakpoint and the preceding one.

t  .378 .403 .433 .961 ; No. of decreasing modes 2  1  1  1

Table 4: Example 1. PSC: Timings with number of intervals per processor fixed at 1000. All times are in seconds.

Number of processors  1  2  4  8  16 32 64 128
Time for initial reduction  .347 .348 .347 .347 .347 .347 .347 .347
Time for reduced system solution .001 .001 .001 .014 .013 .012 .013 .014
Total time  .371 .372 .371 .384 .382 .382 .383 .384

with boundary conditions

Φ(0, μ) = F_{1}(μ) (μ > 0),  Φ(L, μ) = F_{2}(μ) (μ < 0).  (58)

One widely accepted method for solving this equation, the discrete ordinates method, proceeds by replacing the integral term in (57) by a quadrature approximation. The resulting two-point boundary value problem in the remaining independent variable x can then be solved to obtain an approximation to Φ. We use Gaussian quadrature with an even number of abscissae, in which specification of weights ω_{i}, i = 1, ⋯, n, and abscissae μ_{i}, i = 1, ⋯, n with −1 < μ_{1} < ⋯ < μ_{n} < 1, μ_{i} = −μ_{n−i+1}, ω_{i} = ω_{n−i+1}, and ω_{i} > 0 leads to the approximation

∫_{−1}^{1} f(μ) dμ ≈ ∑_{i=1}^{n} ω_{i}f(μ_{i}).

Using the notation

Φ_{i}(x) =^{△} Φ(x, μ_{i}),

and discretizing the boundary conditions in an obvious (though not optimal) way, we can write the boundary value problem arising from Example 2 of Jin and Levermore [7] as

μ_{i}Φ′_{i}(x) + σ_{T} (x)Φ_{i}(x) = \frac{1}{2}w(x)σ_{T} (x) ∑_{j=1}^{n} ω_{j}Φ_{j}(x),  i = 1, ⋯, n,  (59)

Φ_{i}(0) = 5, i = \frac{1}{2}n + 1, ⋯, n,  Φ_{i}(11) = 0, i = 1, ⋯, \frac{1}{2}n.  (60)

Because of separation of the boundary conditions and the well-posed nature of this problem, we would expect exactly half of the fundamental modes to be nonincreasing and half to be nondecreasing across the entire interval. Our implementation of PSC indicates that this is indeed the case. No “spurious” dichotomy changes are flagged on any of the examples we tried.

Table 5: Example 1. SQR: Timings with number of intervals per processor fixed at 1000. All times are in seconds.

Number of processors  1  2  4  8  16 32 64 128
Time for initial reduction  .613 .613 .614 .614 .614 .613 .614 .614
Time for reduced system solution .001 .002 .002 .012 .013 .014 .014 .014
Total time  .642 .643 .644 .654 .655 .655 .656 .656

Table 6: Example 2. PSC: Timings with n = 20 and number of intervals per processor fixed at 200. All times are in seconds.

Number of processors  1  2  4  8  16
Time for initial reduction  2.18 2.22 2.21 2.19 2.19
Time for reduced system solution .028 .044 .088 .113 .093
Total time  2.26 2.31 2.36 2.35 2.33

Table 7: Example 2. SQR: Timings with n = 20 and number of intervals per processor fixed at 200. All times are in seconds.

Number of processors  1  2  4  8  16
Time for initial reduction  6.61 6.38 6.37 6.37 6.38
Time for reduced system solution .016 .049 .140 .190 .206
Total time  6.72 6.53 6.61 6.66 6.68

Results for the two algorithms are given in Tables 6 through 9 for the cases n = 20 and n = 10. A box discretization was used with constant interval length and a fixed number of intervals per processor. As in Example 1, the time taken to solve the reduced system tends to increase with the number of levels in the tree (though not smoothly), but high efficiency is obtained even on a large number of processors.

References

[1] 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.

[2] U. M. Ascher and R. D. Russell, Reformulation of boundary value problems into “standard” form, SIAM Review, 23 (1981), pp. 238–254.

[3] H.-G. Bock, Recent advances in parameter identification techniques for O.D.E., in Numerical Treatment of Inverse Problems in Differential and Integral Equations, P. Deuflhard and E. Hairer, eds., vol. 2 of Progress in Scientific Computing, Birkhauser, Boston, 1983, pp. 95–121.

[4] F. R. de Hoog and R. M. M. Mattheij, An algorithm for solving multi-point boundary value problems, Computing, 38 (1987), pp. 219–234.

[5] , On the conditioning of multipoint and integral boundary value problems, SIAM Journal of Mathematical Analysis and Applications, 20 (1989), pp. 200–214.

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

Table 8: Example 2. PSC: Timings with n = 10 and number of intervals per processor fixed at 200. All times are in seconds.

Number of processors  1  2  4  8  16 32 64 128
Time for initial reduction  .449 .450 .450 .450 .451 .451 .452 .452
Time for reduced system solution .006 .009 .012 .034 .028 .041 .043 .034
Total time  .473 .477 .480 .502 .497 .510 .512 .504

Table 9: Example 2. SQR: Timings with n = 10 and number of intervals per processor fixed at 200. All times are in seconds.

Number of processors  1  2  4  8  16 32 64 128
Time for initial reduction  1.27 1.27 1.27 1.27 1.27 1.26 1.26 1.27
Time for reduced system solution .004 .010 .016 .031 .029 .034 .040 .054
Total time  1.30 1.31 1.32 1.33 1.32 1.33 1.33 1.35

[7] S. Jin and D. Levermore, The discrete-ordinate method in diffusive regimes, Transport Theory and Statistical Physics, 20 (1991), pp. 413–439.

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

[9] R. M. M. Mattheij, Stability of block LU-decompositions of matrices arising from BVP, SIAM Journal on Algebraic and Discrete Methods, 5 (1984), pp. 314–331.

[10] , On boundary value problems for ordinary differential equations with parameters, in Differential Equations, C. M. Dafermos, ed., Marcel Dekker, New York, 1989, pp. 481–489.

[11] J. Rinzel and D. Terman, Propogation phenomena in a bistable reaction-diffusion system, SIAM Journal on Applied Mathematics, 42 (1982), pp. 1111–1137.

[12] S. J. Wright, Stable parallel elimination for boundary value ODEs, Tech. Rep. MCS–P229–0491, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Illinois, April 1991.

[13] , Stable parallel algorithms for two-point boundary value problems, SIAM Journal on Scientific and Statistical Computing, 13 (1992), pp. 742–764.