STABLE PARALLEL ELIMINATION FOR BOUNDARY VALUE ODES^{∗}

STEPHEN J. WRIGHT^{†}

PREPRINT MCS–P229–0491, MCS DIVISION, ARGONNE NATIONAL LABORATORY

Abstract. A parallelizable and vectorizable algorithm for solving linear algebraic systems arising from two-point boundary value ODEs is described. The method is equivalent to Gaussian elimination, with row partial pivoting, applied to a certain row- and column- reordered version of the usual almost-block-diagonal coefficient matrix. Analytical and numerical evidence is presented to show that the algorithm is stable. Results from implementation on a shared-memory multiprocessor and a vector processor are given. The approach can be extended to handle problems with multipoint and integral conditions, or algebraic parameters.

Key words. Parallel algorithms, two-point boundary value problems, Gaussian elimination, stability

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

1. Introduction. In a recent paper [11], we described a technique for solving the almost block-diagonal system of linear equations that arise in various algorithms for solving two-point boundary value problems, such as collocation, multiple shooting, and finite differencing. The advantage of the approach is that it can be adapted to make efficient use of parallel and vector computer architectures, while retaining the stability characteristics that are essential in robust software. The technique makes use of orthogonal transformations that take advantage of the special structure of the matrices.

In this report, we discuss a similar technique based on Gauss transformations. The new “structured elimination” technique and the structured orthogonal technique make identical use of the matrix structure, but the former requires substantially fewer floating-point operations. Timing results are given for some computational experiments on two advanced computer architectures. In Section 4, we investigate the stability of this structured elimination algorithm by relating it to a block factorization algorithm, which effectively decouples the increasing and decreasing fundamental modes when an exponential dichotomy is present and the problem is well conditioned. We show in Section 5 how the approach can be extended to handle problems with algebraic parameters or multipoint and integral side conditions.

Some of the details of the problem and the structured factorization approach have been discussed in [11] and hence are omitted here.

2. The Algorithm. Consider the linear first-order two-point boundary value problem

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

Multiple shooting and finite differencing proceed by choosing a grid

a = t_{1} < t_{2} < … < t_{k+1} = b

and seeking vectors s_{i} that approximate the true solution y(t_{i}) for i = 1, …, k + 1. In both methods, linear systems of the following form must be solved for s_{1}, s_{2}, …, s_{k+1}:

[ \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} ] . (3)

∗ 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. A grant of computer time at the North Carolina Supercomputing Center is gratefully acknowledged.

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

Here, A_{i}, C_{i} ∈ R^{n×n}, f_{i} ∈ R^{n}. Collocation techniques give rise to similar systems, after condensation has been used to remove some of the parameters. In the case of separated end conditions, we can assume that B_{a}, B_{b}, and d have the form

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

where B̄_{a} ∈ R^{l×n} and B̄_{b} ∈ R^{(n−l)×n}. In this case, the system is usually reordered as

[ 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} ] . (4)

Our algorithm starts, as in [11], by choosing an integer P ≤ k/2, and a set of “separator” indices

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

Using these indices, we divide the coefficient matrix into P “slices,” where slice j consists of rows (k_{j−1} + 1)n + 1, …, (k_{j} + 1)n, with the (structural) zero columns deleted. For example, the first slice is

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

We now perform Gaussian elimination, with row partial pivoting, on the submatrix obtained by omitting the first and last block column for each slice. In (5), the first n steps of this process yield the factorization

[ \begin{matrix}U_{1} \\ 0 \end{matrix} ] = L̃_{1,n}P_{1,n} … L̃_{1,2}P_{1,2}L̃_{1,1}P_{1,1} [ \begin{matrix}C_{1} \\ A_{2} \end{matrix} ] , (6)

where P_{1,i} ∈ R^{2n×2n}, i = 1, …, n, are permutation matrices; L̃_{1,i} ∈ R^{2n×2n}, i = 1, …, n, are Gauss transformations; and U_{1} ∈ R^{n×n} is upper triangular. (The information needed to construct each L̃_{1,i} can be stored in the locations formerly occupied by the eliminated elements of [C_{1}^{T} A_{2}^{T}]^{T}, while the permutation data for the whole process requires an integer array with dimension equal to that of the right-hand side of (3).) Applying the transformations across the first 2n rows of the matrix (5), we obtain an equivalent system whose augmented matrix is

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

where

[ \begin{matrix}G_{1} & U_{1} & E_{1} & g_{1} \\ Ḡ_{2} & 0 & C̄_{2} & f̄_{2} \end{matrix} ] = L̃_{1,n}P_{1,n} … L̃_{1,2}P_{1,2}L̃_{1,1}P_{1,1} [ \begin{matrix}A_{1} & C_{1} & 0 & f_{1} \\ 0 & A_{2} & C_{2} & f_{2} \end{matrix} ] .

Repeating the process, for rows n + 1, …, 3n, we find permutation matrices P_{2,i} and Gauss transformations L̃_{2,i}, i = 1, …, n, such that

[ \begin{matrix}G_{2} & R_{2} & E_{2} & g_{2} \\ Ḡ_{3} & 0 & C̄_{3} & f̄_{3} \end{matrix} ] = L̃_{2,n}P_{2,n} … L̃_{2,2}P_{2,2}L̃_{2,1}P_{2,1} [ \begin{matrix}Ḡ_{2} & C̄_{2} & 0 & f̄_{2} \\ 0 & A_{3} & C_{3} & f_{3} \end{matrix} ] .

Proceeding in this way for stages 3, …, k_{1} − 1, we finally obtain the equivalent system

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

A reduced system, involving only (P + 1)n variables and equations, can be formed by taking the last block row of each partition after the factorization above has been applied. As in [11], this system has the form

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

Since the form of (8) is identical to that of the original system, the elimination procedure just described can be applied recursively. Eventually, a 2n × 2n system of the form

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

is obtained. This can be solved by standard Gaussian elimination with partial pivoting. The “intermediate” unknowns can be recovered by back substitution, using the matrices U_{i} and E_{i} and the right-hand sides g_{i} which were obtained earlier in the factorization process.

Different variants of the overall algorithm can be distinguished by the number of levels used to obtain (9) from (3). A serial (one-level) algorithm would set P = 1 and hence go directly from (3) to (9). On a shared-memory machine, it may be appropriate to use a two-level algorithm, in which P is set equal to the number of available processors. The first stage of the reduction (which “compresses” (3) into (8)) can be carried out in parallel, while the compression of (8) to (9) can be performed on a single processor. On a vector processor, a “cyclic reduction” version, which uses log_{2} k levels, is most appropriate. At the first level, we set P = k/2, so that each partition consists of just two block rows. A reduced system of approximately half the size of the original system is produced. The code can be written so that the innermost (vectorized) loops have length k/2. This procedure is applied recursively; at each level the size of the remaining system is halved. For more details, see [11].

For future reference, note that the elimination step (6) can be recast as follows: suppose that l of the pivot rows (0 ≤ l ≤ n) are drawn from the first n rows of the original matrix (i.e., the C_{1} rows), so that n − l of the pivots come from the A_{2} rows. Then we can write

[ \begin{matrix}U_{1} \\ 0 \end{matrix} ] = [ \begin{matrix}L_{1}^{1} & 0 \\ L_{1}^{2} & I \end{matrix} ] [ \begin{matrix}R̄_{1} & 0 \\ 0 & 0 \end{matrix} ] [ \begin{matrix}P_{a} & P_{b} \\ P_{b} & P_{a} \end{matrix} ] [ \begin{matrix}Ō_{0} & 0 \\ 0 & Ō_{1} \end{matrix} ] [ \begin{matrix}C_{1} \\ A_{2} \end{matrix} ] , (10)

where R̄_{1}, Ō_{0} and Ō_{1} are n × n permutation matrices, L_{1}^{1} is n × n and unit lower triangular, L_{1}^{2} is n × n, and

P_{a} = [ \begin{matrix}0 & 0 \\ 0 & I_{l} \end{matrix} ] , P_{b} = [ \begin{matrix}I_{n−l} & 0 \\ 0 & 0 \end{matrix} ] . (11)

Table 1 Operation counts and storage requirements for four algorithms, assuming separated end conditions (k = number of mesh

points, 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)] 3kn^{2}
DECOMP/SOLVE  k[\frac{2}{3}n^{3} + (4R + 5p)n^{2} − 2np^{2}]  2kn^{2}
Structured QR  k[\frac{46}{3}n^{3} + (15R + 30)n^{2}]  4kn^{2}
Structured LU  k[\frac{23}{3}n^{3} + (8R − \frac{5}{2})n^{2}]  4kn^{2}

Table 2 Operation counts and storage requirements for four algorithms, assuming nonseparated end conditions (k = number of

mesh points, 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}
DECOMP/SOLVE k[\frac{14}{3}n^{3} + 4Rn^{2}] 3kn^{2}
Structured QR  k[\frac{46}{3}n^{3} + (15R + 30)n^{2}] 4kn^{2}
Structured LU  k[\frac{23}{3}n^{3} + (8R − \frac{5}{2})n^{2}]  4kn^{2}

Note that

[ \begin{matrix}L_{1}^{1} & 0 \\ L_{1}^{2} & I \end{matrix} ] = L̃_{1,n}L̃_{1,n−1} … L̃_{1,1}.

The notation here is similar to that of Mattheij [9]. This is no accident, since the style of stability analysis in Section 4 is similar to that used in [9].

3. Computational Results. Variants of the structured elimination 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. In this section, we compare these codes with implementations of other algorithms for solution of the linear system arising from finite difference and shooting methods. These are

• a plain row partial pivoting code. The ordering (4) is used if the end conditions are separated. Otherwise, left-hand and coupled conditions are listed first, and right-hand conditions are listed at the end. We denote these codes by ROWPP.

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

• the structured orthogonal factorization algorithm from [11]. These codes are generically denoted by

       SQR.

The structured elimination codes are referred to by the name SLU. The structured codes are subcategorized according to the number of levels. For example, the one-level, two-level, and cyclic reduction variants of the SQR code are named SQR-1, SQR-2, and SQR-CR, respectively.

Table 1 compares storage and operation counts for the four algorithms, assuming separated end conditions. In tabulating storage requirements, we assume that the multipliers and Householder vectors generated during the factorization are stored, for possible later use with different right-hand sides. For the two structured algorithms, the operation counts are independent of the number of recursive levels in the factorization process. Structured elimination requires about a third as many operations as structured orthogonal factorization, and they both have higher counts than the “serial” algorithms. However, operation counts are poor predictors of runtime, particularly on advanced computer architectures with their vectorization and pipelining capabilities. The structured algorithms are not as slow as Table 1 would suggest, even in serial mode. Table 2 contains the same information for the case in which the end conditions are completely coupled. The structured algorithms are indifferent to whether the end conditions are separated or nonseparated, whereas the operation counts for the serial algorithms increase markedly in the nonseparated case.

Table 3 Box method, ∞-norm error in the computed solution for Problem 1A

k = 32 k = 128 k = 1024

                  ROWPP, DECOMP, SQR-1, SLU-1 .28(-1) .17(-2) .26(-4)

Table 4 Multiple shooting, ∞-norm error in the computed solution for Problem 1A

k = 8 k = 16 k = 64

                    ROWPP, DECOMP, SQR-1, SLU-1 .17(-3) .14(-4) .35(-6)

The stability properties of the algorithms are tested by using them to solve two small problems based on the same ODE, but with different end conditions.

Problem 1A 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_{2}(π) = e^{π}
y_{1}(π) + 3y_{3}(π) = 4e^{π}.

The solution is y(t) = e^{t}(1, 1, 1)^{T}. A fundamental solution for this problem is

Y (t) = [ \begin{matrix}sin t & 0 & − cos t \\ 0 & 1 & 0 \\ cos t & 0 & sin t \end{matrix} ] diag(e^{20t}, e^{19t}, e^{−18t}).

Problem 1B (Mattheij [9]) Same as Problem 1A, but with end conditions

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

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

The code dverk is used as the IVP solver in the multiple shooting code, with the global error tolerance set to 10^{−10}. The k subintervals are equally spaced in both codes. A midpoint-rule discretization is used in the finite-difference code (yielding the “box method”). Results from single-precision implementation on the CRAY Y-MP are reported.

Tables 3 and 4 show maximum errors in the computed solutions for Problem 1A. The behavior of the four codes is identical. Since the end conditions are separated, the two “serial” codes use the ordering (4); the matrix is now banded, and any reasonable partial pivoting strategy can be expected to be stable. The errors are consistent with estimates that take into account discretization error (in the case of the finite difference method) and conditioning of the problem.

Tables 5 and 6 show errors in the computed solutions for Problem 1B, which has coupled end conditions. In this respect, the behavior of the codes ROWPP, SQR-1, and SLU-1 was indistinguishable. The code DECOMP does not properly decouple the modes, and so gives large errors.

Table 5 Box method, ∞-norm error in the computed solution for Problem 1B

k = 32 k = 128 k = 1024
ROWPP, SQR-1, SLU-1 .28(-1) .17(-2) .26(-4)
DECOMP  .24(+2) .24(+2) .24(+2)

Table 6 Multiple shooting, ∞-norm error in the computed solution for Problem 1B

k = 8 k = 16 k = 64
ROWPP, SQR-1, SLU-1 .17(-3) .14(-4) .35(-6)
DECOMP  .24(+2) .24(+2) .24(+2)

As expected, the compactification algorithm produces very large errors when applied to both 1A and 1B.

Four test problems are used to compare the relative execution speeds of the codes on the Alliant FX/8 and the CRAY Y-MP:

Problem 2 (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), y_{1}(0) = 1, y_{1}(1) = e,

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

Problem 3 (Brown and Lorenz [3]) 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{ϵ}}.

(We use ϵ = .001.)

Problem 4A a = 0, b = 1, n = 5,

M(t) = [ \begin{matrix}−λ_{1} cos 2ω_{1}t & 0 & ω_{1} + λ_{1} sin 2ω_{1}t & 0 & 0 \\ 0 & −λ_{2} cos 2ω_{2}t & 0 & ω_{2} + λ_{2} sin 2ω_{2}t & 0 \\ −ω_{1} + λ_{1} sin 2ω_{1}t & 0 & λ_{1} cos 2ω_{1}t & 0 & 0 \\ 0 & −ω_{2} + λ_{2} sin 2ω_{2}t & 0 & λ_{2} cos 2ω_{2}t & 0 \\ 0 & 0 & 0 & 0 & λ_{3} \end{matrix} ] ,

q(t) = [q_{1}(t), q_{2}(t), q_{3}(t), q_{4}(t), −(λ_{3} − 1)e^{t}]^{T} ,

with boundary conditions

\begin{aligned}y_{1}(0) & = & 1, \\ y_{2}(0) + 4y_{5}(0) & = & 5, \end{aligned}  \begin{aligned}y_{1}(1) & = & e, \\ −y_{3}(1) + y_{4}(1) & = & 0, \\ −4y_{2}(1) + 5y_{5}(1) & = & e. \end{aligned}

Table 7 Alliant FX/8, one-processor timings for linear system solvers, k = 1024 (times in seconds)

                  Problem ROWPP  DECOMP  SQR-1  SQR-CR  SLU-1  SLU-CR

2  .424  .667  .816  1.06  .576  .649
3  1.31  1.94  3.21  4.50  2.08  2.53
4A  1.94  2.71  5.42  7.79  3.48  4.40
4B  3.44  3.28  5.42  7.78  3.47  4.37

We use λ_{1} = 200, λ_{2} = 50, λ_{3} = 10, ω_{1} = 1, and ω_{2} = 25. The components of q(t) are chosen so that the solution is e^{t}(1, 1, 1, 1, 1)^{T}. This problem has a fundamental solution

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

and is obtained by replicating Problem 2.

Problem 4B Same ODE and solution as Problem 4A, except that three of the boundary conditions are coupled:

y_{1}(0) = 1
y_{2}(0) + 4y_{5}(0) + y_{3}(1) = 5 + e
−5y_{1}(0) + y_{1}(1) = −5 + e
3y_{2}(0) − y_{3}(1) + y_{4}(1) = 3
−4y_{2}(1) + 5y_{5}(1) = e.

Tables 7 and 8 show execution times on one processor of the Alliant FX/8. The -Og option was used during compilation, so the vectorization capabilities of the processor were not used. Comparing these tables with the operation counts in Tables 1 and 2, we see that the time required by the DECOMP code is greater than expected, while the time required by SQR-1 is less than expected. This result is at least partially due to the large amount of data manipulation performed by DECOMP (in the interests of economizing on storage) and the lack of such data manipulation in the SQR code, which avoids the need for pivoting altogether. The “cyclic reduction” codes take somewhat longer than their sequential counterparts, partly because their data access patterns are not localized. On problem 4B (the only one with some coupled end conditions), the time needed by ROWPP and DECOMP increases markedly over the time for Problem 4A. In fact, SLU becomes competitive here even without the benefit of vectorization or parallelization. In addition, we note that DECOMP always gives inaccurate answers for the coupled end-condition Problem 4B, because of its failure to decouple the increasing and decreasing fundamental modes.

Tables 9 and 10 give execution times for the codes when they have been compiled to run on eight processors of the Alliant without vectorization. The structured codes speed up as predicted, although the bottleneck of solving the reduced system on a single processor is particularly noticeable when k = 128. Defining speedup as the ratio of time for fastest serial algorithm to time for best parallel algorithm, we see from Table 11 that the SLU-2 code demonstrates excellent efficiency for the larger problems, especially when the boundary conditions are not separated.

Results for vectorized implementations on the CRAY Y-MP are shown in Tables 12 and 13. The “cyclic reduction” versions of the structured codes vectorize very efficiently, as can be seen by comparing their execution times with the corresponding sequential versions. Speedups over the best serial codes fall off

Table 8

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

Problem ROWPP  DECOMP  SQR-1  SQR-CR  SLU-1  SLU-CR

2  .065  .118  .119  .157  .085  .098
3  .172  .291  .416  .573  .273  .322
4A  .247  .399  .692  .989  .449  .543
4B  .445  .468  .697  .990  .446  .538

Table 9

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

Problem ROWPP  DECOMP  SQR-2  SLU-2

2  .359  .695  .136  .104
3  .749  1.65  .466  .311
4A  .963  2.14  .765  .506
4B  1.37  2.79  .766  .505

Table 10

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

Problem ROWPP  DECOMP  SQR-2  SLU-2

2  .051  .114  .037  .035
3  .096  .242  .091  .069
4A  .124  .301  .143  .100
4B  .181  .387  .142  .100

Table 11

Alliant FX/8. Speedups, defined as time for fastest one-processor algorithm divided by time for fastest eight-processor algorithm

Problem k = 1024 k = 128
2  4.08  1.86
3  4.21  2.49
4A  3.83  2.47
4B  6.50  4.45

Table 12 CRAY Y-MP, vectorized code on one processor, k = 1024 (times in milliseconds)

                  Problem ROWPP  DECOMP  SQR-1  SQR-CR  SLU-1  SLU-CR

2  34.3  48.8  117.  10.8  83.0  7.3
3  84.7  134.  300.  51.6  194.  40.0
4A  111.  171.  415.  89.4  263.  71.4
4B  233.  245.  413.  89.5  261.  71.8

Table 13 CRAY Y-MP, vectorized code on one processor, k = 128 (times in milliseconds)

                  Problem ROWPP  DECOMP  SQR-1  SQR-CR  SLU-1  SLU-CR

2  4.3  6.1  14.6  2.2  10.3  1.6
3  10.6  16.8  37.4  10.1  24.2  7.3
4A  13.8  21.4  51.8  17.3  32.8  12.5
4B  29.2  30.7  51.5  17.3  32.6  12.5

somewhat when the end conditions are separated and n = 4 or n = 5, since in these cases the codes ROWPP and DECOMP achieve a small amount of vectorization in the factorization of each of the blocks in the coefficient matrix in (4).

4. Stability of Structured Elimination. The stability of the SQR algorithms has been examined in [11] by using standard tools for the analysis of orthogonal factorizations from the numerical linear algebra literature (see, for example, Golub and Van Loan [5]). Similar tools can be used to analyze the SLU algorithms, since each variant of SLU is simply Gaussian elimination with row partial pivoting, applied to a row- and column-reordered version of the coefficient matrix in (3). However, it is well known that, even with row pivoting, element growth can occur in the upper triangular factor which is exponential in the size of the matrix (in our case (k + 1)n). Although such growth is unlikely in practice, it is nevertheless desirable to use the structure of the matrix in (3) and the properties of the underlying problem (1), (2) to obtain a more realistic bound on the error growth in the method. We do this by using a style of analysis that was developed by Mattheij in a series of papers (see, for example [8], [9]). Our approach is first to describe a block factorization algorithm, which is similar to “stable compactification” (Section 4.2). This algorithm performs the required decoupling of the growing and decaying modes. Then, in Section 4.3, we show that the “serial” variant of SLU is essentially equivalent to this block algorithm under certain reasonable assumptions.

4.1. Dichotomy and Conditioning. We start by briefly reviewing the salient points of the theory of dichotomy and conditioning, from Ascher, Mattheij, and Russell [2] and Mattheij [9].

The ODE (1) is said to have an exponential dichotomy if there exist real positive constants K, λ, and μ and an integer l with 0 ≤ l ≤ n, such that, for some fundamental solution Y (t) of (1),

∥Y (x)P_{a}Y^{−1}(t)∥ ≤ Ke^{−λ(x−t)},  x ≥ t,
(12)  ∥Y (x)P_{b}Y^{−1}(t)∥ ≤ Ke^{−μ(t−x)},  t ≥ x,

where P_{a} and P_{b} are as in (11). Any other fundamental solution Z(t) of (1) is related to Y (t) by

Z(t) = Y (t)H = Y (t) [ \begin{matrix}H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} ] ,

where H ∈ R^{n×n} is constant and nonsingular, H_{11} ∈ R^{(n−l)×(n−l)}, etc.

Similarly, we can consider the linear recurrence given by the discrete system (3), namely,

A_{i}s_{i} + C_{i}s_{i+1} = f_{i}, i = 1, …, k. (13)

This recurrence is said to have an exponential dichotomy if there exist positive constants K, σ > 1, ρ > 1, and an integer l as above, such that for some fundamental solution {Y_{i}} (with A_{i}Y_{i}+C_{i}Y_{i+1} = 0, i = 1, …, k) we have

∥Y_{i}P_{a}Y_{j+1}^{−1}∥ ≤ Kσ^{j−i},  i ≥ j + 1,
(14)  ∥Y_{i}P_{b}Y_{j}^{−1}∥ ≤ Kρ^{i−j},  i ≤ j.

As in the continuous case, any other fundamental solution {Z_{i}} of (13) can be related to {Y_{i}} by

Z_{i} = Y_{i}H = Y_{i} [ \begin{matrix}H_{11} & H_{12} \\ H_{21} & H_{22} \end{matrix} ] , i = 1, …, k. (15)

{Z_{i}} is said to be consistent with {Y_{i}} if H_{11} in (15) is nonsingular. Since the exact dichotomic solution {Y_{i}} is not obtainable in practice, any consistent solution {Z_{i}} can in theory be used to decouple the increasing and decreasing fundamental modes. It is sufficient to define consistency constant L by

L =^{def} \frac{∥Y_{1}^{2}H_{21}∥}{min_{∥u∥=1} ∥Y_{1}^{1}H_{11}u∥}

(where Y_{1}^{1} and Y_{1}^{2} are the first (n − l) and last l columns of Y_{1}, respectively) and assume that L is of “moderate size.”

It is known [2, p. 289] that if (3) is obtained by applying finite differencing to (1) with uniform mesh spacing h ≡ t_{i+1} −t_{i}, and if the ODE (1) has the dichotomy (12), then the recurrence (13) has the dichotomy (14) with σ = e^{λh} and ρ = e^{μh}.

From here until the end of the paper, we work solely with (3) and (13) and assume that the dichotomy (14) applies. We further assume that C_{i} = −I, i = 1, …, k. This assumption is true for multiple shooting and the condensed form of the collocation equations; it can be obtained for finite differencing by multiplying (13) through by −C_{i}^{−1}. We can therefore rewrite (13) as

s_{i+1} = A_{i}s_{i} − f_{i}. (16)

Suppose that we choose some orthogonal matrix Q_{0} ∈ R^{n×n} and generate a sequence {Q_{i}}, i = 0, …, k of orthogonal matrices by using the formula

A_{i}Q_{i−1} = Q_{i}Δ_{i}. (17)

We require that Δ_{i} have the form

Δ_{i} = [ \begin{matrix}Δ_{i}^{11} & Δ_{i}^{12} \\ 0 & Δ_{i}^{22} \end{matrix} ] . (18)

This could be obtained, for example, by defining Q_{i}, Δ_{i} as QR factors of A_{i}Q_{i−1}. If {Z_{i}} is a fundamental solution of (16), we can use (17) to write

Z_{i+1} = A_{i}Z_{i} ⇒ (Q_{i}^{T} Z_{i+1}) = Δ_{i}(Q_{i−1}^{T}Z_{i}),

so that W_{i} = Q_{i−1}^{T}Z_{i}, W_{1} = I, defines a fundamental solution of the homogeneous recurrence

w_{i+1} = Δ_{i}w_{i}, i = 1, …, k.

It follows from (18) that each W_{i} has the form

W_{i} = [ \begin{matrix}W_{i}^{11} & W_{i}^{12} \\ 0 & W_{i}^{22} \end{matrix} ] , i = 1, …, k + 1. (19)

Note, in particular, that

W_{i}^{11}(W_{j+1}^{11})^{−1} = ∏_{r=i}^{j}(Δ_{r}^{11})^{−1},  i ≤ j,
W_{i}^{22}(W_{j+1}^{22})^{−1} = Δ_{i−1}^{22}Δ_{i−2}^{22} … Δ_{j+1}^{22},  i ≥ j + 1.

As a special case of [2, Theorem 6.30], we have the following:

Theorem 4.1. Suppose the dichotomy (14) holds, with fundamental solution {Y_{i}}, and that the fundamental solution {Z_{i}} with initial condition Z_{1} = Q_{0} is consistent with {Y_{i}}, with L of moderate size. Then there are moderately sized constants K_{1} > 1 and K_{2} > 1 such that

∥ ∏_{r=i}^{j}(Δ_{r}^{11})^{−1}∥_{2} ≤ K_{1}ρ^{i−j},  i ≤ j,
∥Δ_{i−1}^{22}Δ_{i−2}^{22} … Δ_{j+1}^{22}∥_{2} ≤ K_{2}σ^{j−i},  i ≥ j + 1.

4.2. A Structured Block Factorization. Given the assumption C_{i} ≡ −I, we can negate all but the boundary condition rows of (3) and perform an obvious row and column reordering to obtain the matrix

A = [ \begin{matrix}I & & & & & 0 & −A_{1} \\ −A_{2} & I & & & & 0 & 0 \\ & −A_{3} & I & & & ⋮ & ⋮ \\ & & ... & ... \\ & & & & I & 0 & 0 \\ & & & & −A_{k} & I & 0 \\ & & & & & B_{b} & B_{a} \end{matrix} ] . (20)

In this section, we describe a block factorization algorithm for A and show that it gives a stable way of solving the system (3).

Using the orthogonal matrices Q_{0}, Q_{1}, …, Q_{k} of the previous section, define

R̃ = diag(Q_{1}^{T} , Q_{2}^{T} , …, Q_{k}^{T} , I)
D̃ = diag(Q_{1}, Q_{2}, …, Q_{k}, Q_{0}).

From (17) and (18) it follows that R̃AD̃ has the form

R̃AD̃ = [ \begin{matrix}I & 0 & −Δ_{1}^{11} & −Δ_{1}^{12} \\ 0 & I & 0 & −Δ_{1}^{22} \\ −Δ_{2}^{11} & −Δ_{2}^{12} & I & 0 \\ 0 & −Δ_{2}^{22} & 0 & I \\ & & −Δ_{3}^{11} & −Δ_{3}^{12} & I & 0 \\ & & 0 & −Δ_{3}^{22} & 0 & I \\ & & & ... & ... & ... & ... \\ & & & & & & I & 0 \\ & & & & & & 0 & I \\ −Δ_{k}^{11} & −Δ_{k}^{12} & I & 0 & 0 & 0 \\ 0 & −Δ_{k}^{22} & 0 & I & 0 & 0 \\ & & & & & & & & (B_{b}Q_{k}) & (B_{a}Q_{0}) \end{matrix} ] .

Using the n×n permutation matrices P_{a} and P_{b} from (11) as building blocks, we define (k +1)n-dimensional permutation matrices P_{i}, i = 1, …, (k − 1), by

P_{i} = [ \begin{matrix}I \\ & . & . & . \\ & & & & I \\ & & & & & P_{a} & P_{b} (block row i) \\ & & & & & P_{b} & P_{a} (block row i + 1) \\ & & & & & & I \\ . & . & . \\ & & & & & & I \end{matrix} ] .

Let

P = P_{k−1}P_{k−2} … P_{1}.

Applying these permutations to the rows of R̃AD̃, we obtain

P R̃AD̃ =

[ \begin{matrix}−Δ_{2}^{11} & −Δ_{2}^{12} & I & 0 & 0 & 0 \\ 0 & I & 0 & 0 & 0 & −Δ_{1}^{22} \\ & 0 & −Δ_{3}^{11} & −Δ_{3}^{12} & I \\ & −Δ_{2}^{22} & 0 & I & 0 \\ & & & 0 & −Δ_{4}^{11} & −Δ_{4}^{12} & I \\ & & & −Δ_{3}^{22} & 0 & I & 0 \\ & & & & ... & ... & ... & ... \\ −Δ_{k}^{11} & −Δ_{k}^{12} & I & 0 \\ 0 & I & 0 & 0 \\ I & 0 & 0 & 0 & 0 & 0 & −Δ_{1}^{11} & −Δ_{1}^{12} \\ 0 & 0 & 0 & −Δ_{k}^{22} & 0 & I & 0 & 0 \\ & & & & & & & & (B_{b}Q_{k}) & (B_{a}Q_{0}) \end{matrix} ] .

Now, without further pivoting, we can perform block elimination in an obvious way, stopping when the final 2n × 2n submatrix is reached. We thus obtain

P R̃AD̃ = L̃′Ũ′,

where

L̃′ = [ \begin{matrix}I \\ L̃_{1} & I \\ & L̃_{2} & I \\ & & ... & ... \\ & & & L̃_{k−2} & I \\ H̃_{1} & & & H̃_{k−2} & Ĥ_{k−1} & I \\ & & & & & & I \end{matrix} ] ,

Ũ′ = [ \begin{matrix}Ũ_{1} & Ẽ_{1} & & & & G̃_{1} \\ & Ũ_{2} & Ẽ_{2} & & & G̃_{2} \\ & & ... & ... & & ⋮ \\ & & & Ũ_{k−1} & Ẽ_{k−1} & G̃_{k−1} \\ & & & & R̂_{1} & R̂_{2} \\ & & & & B_{b}Q_{k} & B_{a}Q_{0} \end{matrix} ] ,

L̃_{i} = [ \begin{matrix}0 & 0 \\ 0 & −Δ_{i+1}^{22} \end{matrix} ] ,  Ũ_{i} = [ \begin{matrix}−Δ_{i+1}^{11} & −Δ_{i+1}^{12} \\ 0 & I \end{matrix} ] ,  Ẽ_{i} = [ \begin{matrix}I & 0 \\ 0 & 0 \end{matrix} ] ,

H̃_{i} = [ \begin{matrix}−Δ_{1}^{11}S_{i+1}^{−1} & −Δ_{1}^{11}S_{i+1}^{−1}Δ_{i+1}^{12} \\ 0 & 0 \end{matrix} ] , G̃_{i} = [ \begin{matrix}0 & 0 \\ 0 & −T_{i} \end{matrix} ] i = 1, …, k − 1,

Ĥ_{k−1} = H̃_{k−1} + L̃_{k−1},
R̂_{1} = [ \begin{matrix}−Δ_{1}^{11}S_{k}^{−1} & 0 \\ 0 & I \end{matrix} ] ,  R̂_{2} = [ \begin{matrix}−Δ_{1}^{11} & −Δ_{1}^{11}X \\ 0 & −T_{k} \end{matrix} ] ,
S_{i} = Δ_{i}^{11}Δ_{i−1}^{11} … Δ_{1}^{11}, T_{i} = Δ_{i}^{22}Δ_{i−1}^{22} … Δ_{1}^{22},

X = S_{1}^{−1}Δ_{1}^{12} − ∑_{j=1}^{k−1} S_{j+1}^{−1}Δ_{j+1}^{12}T_{j}.

Finally, we can find a 2n × 2n permutation matrix Π̃, unit lower triangular ˜ L̃ and upper triangular ˜ Ũ such that

˜L̃ ˜Ũ = Π̃ [ \begin{matrix}R̂_{1} & R̂_{2} \\ B_{b}Q_{k} & B_{a}Q_{0} \end{matrix} ] .

Incorporating this into the overall factorization, we obtain

[ \begin{matrix}I & 0 \\ 0 & Π̃ \end{matrix} ] P R̃AD̃ = L̃Ũ, (21)

where

L̃ = [ \begin{matrix}I \\ L̃_{1} & I \\ & L̃_{2} & I \\ & & ... & ... \\ & & & L̃_{k−2} & I \\ H̃_{1}^{1} & H̃_{2}^{1} & … & H̃_{k−2}^{1} & Ĥ_{k−1}^{1} \\ H̃_{1}^{2} & H̃_{2}^{2} & … & H̃_{k−2}^{2} & Ĥ_{k−1}^{2} & ˜L̃ \end{matrix} ]

Ũ = [ Ũ_{1} Ẽ_{1}  G̃_{1} ; Ũ_{2} Ẽ_{2}  G̃_{2} ; ... ...  ⋮ ; Ũ_{k−1} Ẽ_{k−1} G̃_{k−1} ; ˜Ũ ] ,

and

[ H̃_{i}^{1} ; H̃_{i}^{2} ] = Π̃ [ H̃_{i} ; 0 ] , i = 1, …, k − 1, [ Ĥ_{k−1}^{1} ; Ĥ_{k−1}^{2} ] = Π̃ [ Ĥ_{k−1} ; 0 ] .

The entire process is almost identical to stabilized compactification (as described in [2, pp. 157–161]), the only difference being in the handling of the final 2n × 2n block. Stability of the factorization (21) follows immediately from this relationship. Alternatively, from the standard viewpoint of numerical linear algebra, stability follows from the fact that element growth in the L̃ and Ũ factors is bounded by reasonable quantities. To see this, define

τ =^{def} max(1, ∥A_{1}∥_{2}, …, ∥A_{k}∥_{2}, ∥B_{a}∥_{2}, ∥B_{b}∥_{2}),

and note that

τ ≥ ∥A_{i}∥_{2} = ∥Δ_{i}∥_{2} ≥ ∥Δ_{i}^{12}∥_{2},
τ ≥ ∥Δ_{i}^{11}∥_{2}, τ ≥ ∥Δ_{i}^{22}∥_{2}, i = 1, …, k.

We have from Theorem 4.1 that

∥S_{i}^{−1}∥_{2} ≤ (ρK_{1})ρ^{−i},
∥T_{i}∥_{2} ≤ (K_{2}/σ)σ^{−i}.

Hence

∑_{i=1}^{k−2} ∥H̃_{i}∥_{2} + ∥Ĥ_{k−1}∥_{2} ≤ ∑_{i=1}^{k−1} [∥Δ_{1}^{11}S_{i+1}^{−1}∥_{2} + ∥Δ_{1}^{11}S_{i+1}^{−1}Δ_{i+1}^{12}∥_{2}] + ∥Δ_{k}^{22}∥_{2}

≤ (τ + τ^{2}) ∑_{i=1}^{k−1} ∥S_{i+1}^{−1}∥_{2} + τ
≤ (τ + τ^{2})K_{1} \frac{ρ}{ρ − 1} + τ.

In addition,

∥X∥_{2} ≤ (K_{1}ρ)(K_{2}/σ)τ [ρ^{−1} + ρ^{−2}σ^{−1} + … + ρ^{−k}σ^{−k+1}] ≤ \frac{ρK_{1}K_{2}τ}{ρσ − 1} ,
∥L̃_{i}∥_{2} ≤ τ,
∥Ũ_{i}∥_{2} ≤ 3τ,
∥Ẽ_{i}∥_{2} = 1,
∥G̃_{i}∥_{2} ≤ (K_{2}/σ)σ^{−i}.

Note that if σ = e^{λh} and ρ = e^{μh}, where h is a finite difference mesh width, and if λh ≪ 1 and μh ≪ 1, then terms such as 1/(ρ − 1) and 1/(ρσ − 1) are O(h^{−1}).

For future reference, we state the result as a theorem:

Theorem 4.2. Provided that the assumptions of Theorem 4.1 are satisfied, the factorization (21) is stable.

4.3. Relating the Block Factorization to SLU. We now examine the “serial” version of structured elimination and show how it is related to a block factorization for a particular choice of Q_{0}. This relationship can be used to prove the stability of SLU.

Recall that (10) and (11) depict the first of k + 1 stages of the factorization process. Note that in this stage, a total of n − l rows were exchanged between the first and second block rows. Throughout the remainder of the paper, we make the following assumption:

Assumption 1. The ODE (1) has an exponential dichotomy, with l decreasing fundamental modes (1 ≤ l ≤ n). At each stage i of the structured elimination algorithm (for i = 1, …, k), exactly n − l rows are exchanged, by row pivoting, between the i-th block row and the (i+1)-st block row. Moreover, the collection of n − l rows that is swapped out of the first block row is passed down, intact, through block rows i = 2, 3, …, k.

This assumption was always observed to hold for the examples reported in Section 3. It means that if, for example, row 1 of block 1 is swapped into block 2 during the first stage, then this same row is eventually pivoted to the bottom of the matrix, that is, to somewhere in the k-th or (k + 1)-st block rows. We stress that this does not mean that “reimbedding,” as described in Dieci, Osborne, and Russell [4] and Keller and Lentini [6], does not occur. On the contrary, in all the test examples, we observe reimbedding in the form of occasional changes in the internal pivot sequence from one stage to the next. The last statement of Assumption 1 refers to the collection of rows that are exchanged between successive stages, rather than the pivot sequence within a stage.

Using Assumption 1, we can view stage i of the factorization process (for i = 1, 2, …, k − 1) as a premultiplication of the remaining submatrix by

L_{i}R_{i}P_{i}O_{i},

where, for i = 1, 2, …, k − 1, P_{i} is as defined in Section 4.2, and

L_{i} = [ \begin{matrix}I \\ & . & . & . \\ & & & & I \\ & & & & & L_{i}^{1} & & (block row i) \\ & & & & & L_{i}^{2} & I & (block row i + 1) \\ . & . & . \\ . & . & . \\ & & & & & & & I \end{matrix} ] ,

R_{i} = [ \begin{matrix}I \\ & . & . & . \\ & & & & I \\ & & & & & R̄_{i} & (block row i) \\ & & & & & & I \\ . & . & . \\ & & & & & & I \end{matrix} ] ,

O_{i} = [ \begin{matrix}I \\ & . & . & . \\ & & & & I \\ & & & & & Ō_{i} & (block row i + 1) \\ & & & & & & I \\ . & . & . \\ & & & & & & I \end{matrix} ] .

Here, L_{i}^{1} is unit lower triangular, and R̄_{i} and Ō_{i} are permutation matrices. Before the very first stage, the matrix A is multiplied by some initial permutation O_{0} = diag (Ō_{0}, I, I, …, I). For the final stage, in which the remaining 2n × 2n block is factored, we define a 2n × 2n permutation matrix R̄_{k}, and premultiply by L_{k}R_{k}, where

L_{k} = [ \begin{matrix}I \\ & . & . & . \\ & & & & I \\ & & & & & L_{k}^{1} \end{matrix} ] ,

R_{k} = [ I ; ... ; I ; R̄_{k} ] ,

and L_{k}^{1} is 2n × 2n unit lower triangular. At the end of the process, a block upper triangular matrix U is obtained, with

(L_{k}R_{k})(L_{k−1}R_{k−1}P_{k−1}O_{k−1}) … (L_{1}R_{1}P_{1}O_{1})O_{0}A = U, (22)

where A is as in (20).

We now rework (22) so that it is more directly comparable with the block factorization (21). To start with, simple manipulation involving the permutation blocks in R_{i} and O_{i} yields

\begin{aligned}& (L_{k}R_{k})(L_{k−1}R_{k−1}P_{k−1}O_{k−1}) … (L_{1}R_{1}P_{1}O_{1})O_{0} \\ = & R(L̄_{k} ¯L̄_{k−1}P_{k−1} ¯L̄_{k−2}P_{k−2} … P_{2} ¯L̄_{1}P_{1})O, \end{aligned} (23)

where

R = diag (R̄_{1}, R̄_{2}, …, R̄_{k−1}, R̄_{k})
O = diag (Ō_{0}, Ō_{1}, Ō_{2}, …, Ō_{k−1}, I),

and, for i = 1, …, k − 1, ¯ L̄_{i} has the form

¯L̄_{i} = [ I ; ... ; I ; L̄_{i}^{1} ; L̄_{i}^{2} I ; I ; ... ; I ] ,

where

L̄_{i}^{1} = R̄_{i}^{T} L_{i}^{1}R̄_{i}
L̄_{i}^{2} = L_{i}^{2}R̄_{i}.

For the final elimination matrix, we have used the notation

L̄_{k} = R_{k}^{T} L_{k}R_{k} = diag (I, I, …, I, L̄_{k}^{1}),  where  L̄_{k}^{1} = R̄_{k}^{T} L_{k}^{1}R̄_{k}.

We can then use the following result, which is similar to property 8.56 of Mattheij [9], to isolate the “row exchange” matrices P_{1}, …, P_{k}. Its proof is simple and hence is not included here.

Lemma 4.3. For i = 1, …, k − 2,

P_{k−1}P_{k−2} … P_{i+1} ¯L̄_{i} = L̄_{i}P_{k−1}P_{k−2} … P_{i+1},

where

L̄_{i} = [ I ; ... ; I ; L̄_{i}^{1} ; P_{a}L̄_{i}^{2} I ; I ; ... ; I ; P_{b}L̄_{i}^{2} 0 0 … I 0 ; I ] .

Lemma 4.3 can be used, together with (22) and (23) and the definition of P, to yield

RL^{−1}P OA = U, (24)

where, defining L̄_{k} = ¯ L̄_{k−1} for consistency, we have

L^{−1} = L̄_{k}L̄_{k−1}L̄_{k−2} … L̄_{1}.

Apart from some permutations internal to the blocks, the back-substitution part of the solution process can be viewed as premultiplication by L^{−1}P. Suppose the vector z ∈ R^{(k+1)n} is partitioned as

z = (z_{b}^{1}, z_{a}^{1}, z_{b}^{2}, z_{a}^{2}, …, z_{b}^{k}, z_{a}^{k}, z_{b}^{k+1}, z_{a}^{k+1})^{T} ,

where z_{b}^{i} ∈ R^{n−l} and z_{a}^{i} ∈ R^{l}. Then

P z = (z_{b}^{2}, z_{a}^{1}, z_{b}^{3}, z_{a}^{2}, …, z_{b}^{k}, z_{a}^{k−1}, z_{b}^{1}, z_{a}^{k}, z_{b}^{k+1}, z_{a}^{k+1})^{T} .

For i = 2, …, k−1, the subvector (z_{b}^{i+1}, z_{a}^{i}) is affected only by premultiplication by L̄_{i−1} and L̄_{i}. Meanwhile, the next-to-last subvector (z_{b}^{1}, z_{a}^{k}) is affected by every L̄_{i}. We use the following notation:

• (z_{b}^{i+1}, z_{a}^{i}) is transformed to (z̄_{b}^{i+1}, z̄_{a}^{i}) when premultiplication by L̄_{i−1} is performed;

• (^{i}z_{b}^{1},^{i} z_{a}^{k}) denotes the transformed version of (z_{b}^{1}, z_{a}^{k}), after L̄_{1}, …, L̄_{i} have been applied.

For later reference, we also set

z̃^{i} = L̄_{i}^{1} [ \begin{matrix} & z̄_{b}^{i+1} \\ & z̄_{a}^{i} \end{matrix} ] , i = 1, …, k − 1, [ \begin{matrix} & z̃^{k} \\ & z̃^{k+1} \end{matrix} ] = L̄_{k}^{1} [ \begin{matrix}^{k−1}z_{b}^{1} \\ ^{k−1}z_{a}^{k} \\ z_{b}^{k+1} \\ z_{a}^{k+1} \end{matrix} ] , (25)

so that z̃ = (z̃^{1}, …, z̃^{k+1})^{T} is the final outcome of the application of L^{−1}P. Note that since the L̄_{i}^{1} (i = 1, …, k) are row and column permutations of unit lower triangular matrices, they are certainly nonsingular.

Using the notation of Lemma 4.3, we have, upon premultiplication by L̄_{i−1}, that

(26)  [ \begin{matrix} & z̄_{b}^{i+1} \\  & z̄_{a}^{i} \end{matrix} ] = P_{a}L̄_{i−1}^{2} [ \begin{matrix} & z̄_{b}^{i} \\  & z̄_{a}^{i−1} \end{matrix} ] + [ \begin{matrix}z_{b}^{i+1} \\ z_{a}^{i} \end{matrix} ]
= [ \begin{matrix}0 & 0 \\ (L̄_{i−1}^{2})_{21} & (L̄_{i−1}^{2})_{22} \end{matrix} ] [ \begin{matrix} & z̄_{b}^{i} \\  & z̄_{a}^{i−1} \end{matrix} ] + [ \begin{matrix}z_{b}^{i+1} \\ z_{a}^{i} \end{matrix} ] .

(Here, as earlier, we have used (H)_{11} to denote the leading (n − l) × (n − l) submatrix of a generic n × n matrix H, with the remaining blocks H_{12}, H_{21}, and H_{22} defined accordingly.) Examining the last block row, we see that

(27)  [ \begin{matrix}^{(i−1)}z_{b}^{1} \\ ^{(i−1)}z_{a}^{k+1} \end{matrix} ] = P_{b}L̄_{i−1}^{2} [ \begin{matrix} & z̄_{b}^{i} \\  & z̄_{a}^{i−1} \end{matrix} ] + [ \begin{matrix}^{(i−2)}z_{b}^{1} \\ ^{(i−2)}z_{a}^{k} \end{matrix} ]
= [ \begin{matrix}(L̄_{i−1}^{2})_{11} & (L̄_{i−1}^{2})_{12} \\ 0 & 0 \end{matrix} ] [ \begin{matrix} & z̄_{b}^{i} \\  & z̄_{a}^{i−1} \end{matrix} ] + [ \begin{matrix}^{(i−2)}z_{b}^{1} \\ ^{(i−2)}z_{a}^{k} \end{matrix} ] .

By combining and assembling the nontrivial parts of (26) and (27), we find that

[ \begin{matrix}^{(i−1)}z_{b}^{1} \\ & z̄_{a}^{i} \end{matrix} ] = [ \begin{matrix}I & (L̄_{i−1}^{2})_{12} \\ 0 & (L̄_{i−1}^{2})_{22} \end{matrix} ] [ \begin{matrix}^{(i−2)}z_{b}^{1} \\ & z̄_{a}^{i−1} \end{matrix} ] + [ \begin{matrix}(L̄_{i−1}^{2})_{11} & 0 \\ (L̄_{i−1}^{2})_{21} & I \end{matrix} ] [ \begin{matrix} & z̄_{b}^{i} \\ z_{a}^{i} \end{matrix} ] . (28)

We now show that (28) is a stable recurrence by relating it to the block factorization described in Section 4.2. The block factorization can be initialized in the following way: Define Q_{1} = Ō_{0}^{T}, and then choose Q_{0} in such a way that Q_{1}^{T} A_{1}Q_{0} is upper triangular. The remaining matrices Q_{2}, …, Q_{k} can then be obtained as in (17) and (18), leading to the factorization (21). We need the following assumption:

Assumption 2. The matrix Q_{0} obtained by the procedure just described induces a consistent fundamental solution {Z_{i}} (with Z_{1} = Q_{0}) of the recurrence (16), with a consistency constant of moderate size.

Intuitively, it seems almost certain that this assumption will be satisfied. The first n − l rows of the reordered matrix Q_{1}^{T} A_{1} = Ō_{0}A_{1} tend to correspond to the “larger” rows that are exchanged with the second stage. Postmultiplication by Q_{0} yields Δ_{1}^{11} and Δ_{1}^{22}, which can be expected to satisfy the inequalities of Theorem 4.1, with K_{1} and K_{2} not much larger than 1.

Comparing (21) with (24), we find that

RL^{−1} (P OR̃^{T} P^{T} [ \begin{matrix}I & 0 \\ 0 & Π̃^{T} \end{matrix} ]) L̃ = UD̃Ũ^{−1}. (29)

Clearly, by inspection of the structures of D̃, U and Ũ, we see that both sides of this equation are block upper triangular matrices. We examine the structure of the left-hand side more closely. Using the notation

V_{i} = Ō_{i−1}Q_{i}, i = 1, 3, …, k, (V_{1} = Ō_{0}Q_{1} = I), V_{k+1} = I,

we have OR̃^{T} = diag (V_{1}, V_{2}, …, V_{k}, V_{k+1}). Hence

P OR̃^{T} P^{T} = (30)

[ \begin{matrix}(V_{2})_{11} & 0 & 0 & (V_{2})_{12} \\ 0 & (V_{1})_{22} & 0 & 0 \\ 0 & 0 & (V_{3})_{11} & 0 & 0 & (V_{3})_{12} \\ (V_{2})_{21} & 0 & 0 & (V_{2})_{22} & 0 & 0 \\ & & 0 & 0 \\ & & (V_{3})_{21} & 0 & & ... & ... \\ & & & & ... & ... & ... \\ & & & & & & & (V_{k})_{11} & 0 & 0 & (V_{k})_{12} \\ & & & & & & & 0 & (V_{k−1})_{22} & 0 & 0 \\ & & & & & & & 0 & 0 & (V_{1})_{11} & 0 \\ & & & & & & & (V_{k})_{21} & 0 & 0 & (V_{k})_{22} \\ & & & & & & & & & & & I \end{matrix} ] ,

where the lower right identity block is n × n. Further, from the definition of L̃ in Section 4.2, it is easy to see that

[ \begin{matrix}I & 0 \\ 0 & Π̃^{T} \end{matrix} ] L̃ = [ \begin{matrix}I \\ L̃_{1} & I \\ & L̃_{2} & I \\ & & ... & ... \\ & & & & I \\ H̃_{1} ; 0 & H̃_{2} ; 0 & \begin{matrix}… \\ … \end{matrix} & H̃_{k−2} ; 0 & Ĥ_{k−1} ; 0 & Π̃^{T} ˜ L̃ \end{matrix} ] . (31)

When the matrices (30) and (31) are multiplied together, we obtain the following block lower Hessenberg matrix, which we denote by J^{0}:

J^{0} = (P OR̃^{T} P^{T}) [ \begin{matrix}I & 0 \\ 0 & Π̃^{T} \end{matrix} ] L̃ =

[ \begin{matrix}J_{11}^{0} & J_{12}^{0} \\ J_{21}^{0} & J_{22}^{0} & J_{23}^{0} \\ & ... & ... & ... \\ & & & & J_{k−2,k−1}^{0} \\ & & & J_{k−1,k−2}^{0} & J_{k−1,k−1}^{0} & J_{k−1,k}^{0} \\ J_{k,1}^{0} & J_{k,2}^{0} & … & … & J_{k,k−1}^{0} & J_{k,k}^{0} & J_{k,k+1}^{0} \\ 0 & 0 & … & … & 0 & J_{k+1,k}^{0} & J_{k+1,k+1}^{0} \end{matrix} ] .

Each J_{r,s}^{0} (1 ≤ r, s ≤ k + 1) denotes an n × n block, and they are defined as follows:

J_{i,i}^{0} = [ \begin{matrix}(V_{i+1})_{11} & −(V_{i+1})_{12}Δ_{i+1}^{22} \\ 0 & (V_{i})_{22} \end{matrix} ] , i = 1, …, k − 1,
J_{i,i+1}^{0} = [ \begin{matrix}0 & (V_{i+1})_{12} \\ 0 & 0 \end{matrix} ] , i = 1, …, k − 2,
J_{i+1,i}^{0} = [ \begin{matrix}0 & 0 \\ (V_{i+1})_{21} & −(V_{i+1})_{22}Δ_{i+1}^{22} \end{matrix} ] , i = 1, …, k − 2,
J_{k,i}^{0} = [ \begin{matrix}−(V_{1})_{11}Δ_{1}^{11}S_{i+1}^{−1} & −(V_{1})_{11}Δ_{1}^{11}S_{i+1}^{−1}Δ_{i+1}^{12} \\ 0 & 0 \end{matrix} ] , i = 1, …, k − 2,

J_{k,k−1}^{0} = [ \begin{matrix}−(V_{1})_{11}Δ_{1}^{11}S_{k}^{−1} & −(V_{1})_{11}Δ_{1}^{11}S_{k}^{−1}Δ_{k}^{12} \\ (V_{k})_{21} & −(V_{k})_{22}Δ_{k}^{22} \end{matrix} ] ,
J_{k−1,k}^{0} = [ \begin{matrix}0 & (V_{k})_{12} \\ 0 & 0 \end{matrix} ] Π̃^{T} ˜ L̃,

[ \begin{matrix}J_{k,k}^{0} & J_{k,k+1}^{0} \\ J_{k+1,k}^{0} & J_{k+1,k+1}^{0} \end{matrix} ] = [ \begin{matrix}(V_{1})_{11} \\ & (V_{k})_{22} \\ & & I \end{matrix} ] Π̃^{T} ˜ L̃.

We now define a sequence of matrices J^{i}, i = 0, 1, …, k by

J^{i} = L̄_{i} … L̄_{1} J^{0}. (32)

The following lemma summarizes some observations about each J^{i}:

Lemma 4.4.

(i) J_{i−2,i}^{j} = 0, j = 0, 1, …, k + 1, i = 3, 4, …, k + 1; (ii) J_{i+1,i}^{i} = 0, i = 1, 2, …, k; (iii) J_{j+2,i}^{j} = J_{j+2,i}^{0}, j = 0, 1, …, k − 3, i = 0, 1, …, k + 1. (iv) J_{j,i}^{k} = J_{j,i}^{j} for j = 1, 2, …, k − 1 and i = 0, 1, …, k + 1. Proof. (i) follows from the fact that each of the “block elementary” matrices L̄_{i} is lower triangular, and

so each J^{i} is block lower Hessenberg. For (ii), note from (29) that J^{k} is block upper triangular, and so J_{i+1,i}^{k} = 0. Premultiplication by L̄_{i} affects only block rows i, i + 1 and k of J^{i−1}, and no structural fill-in is created in block row i. It follows that J_{i+1,i}^{k} = J_{i+1,i}^{k−1} = … = J_{i+1,i}^{i} = 0. Note for (iii) that transformations L̄_{1}, …, L̄_{j} affect only block rows 1, 2, …, j + 1 and k of the target matrix J^{0}. In particular, block row j + 2 is unaltered. Assertion (iv) follows from the fact that rows 1, …, j are unaffected by the transformations L̄_{j+1}, …, L̄_{k}.

From Lemma 4.4 (ii), (iii), we have for i = 1, …, k − 2 that

0 = J_{i+1,i}^{i} = P_{a}L̄_{i}^{2}J_{i,i}^{i−1} + J_{i+1,i}^{i−1}
(33)  = [ \begin{matrix}0 & 0 \\ (L̄_{i}^{2})_{21} & (L̄_{i}^{2})_{22} \end{matrix} ] J_{i,i}^{i−1} + [ \begin{matrix}0 & 0 \\ (V_{i+1})_{21} & −(V_{i+1})_{22}Δ_{i+1}^{22} \end{matrix} ] .

For i = 2, …, k − 2,

J_{i,i}^{i−1} = P_{a}L̄_{i−1}^{2}J_{i−1,i}^{i−2} + J_{i,i}^{i−2}
(34)  = [ \begin{matrix}0 & 0 \\ (L̄_{i−1}^{2})_{21} & (L̄_{i−1}^{2})_{22} \end{matrix} ] J_{i−1,i}^{i−2} + [ \begin{matrix}(V_{i+1})_{11} & −(V_{i+1})_{12}Δ_{i+1}^{22} \\ 0 & (V_{i})_{22} \end{matrix} ] .

For i = 3, …, k − 2,

J_{i−1,i}^{i−2} = P_{a}L̄_{i−2}J_{i−2,i}^{i−3} + J_{i−1,i}^{i−3}.

By Lemma 4.4 (i), J_{i−2,i}^{i−3} = 0, and so

J_{i−1,i}^{i−2} = J_{i−1,i}^{i−3} = [ \begin{matrix}0 & (V_{i})_{12} \\ 0 & 0 \end{matrix} ] . (35)

By combining (35) and (34), we obtain

J_{i,i}^{i−1} = [ \begin{matrix}(V_{i+1})_{11} & −(V_{i+1})_{12}Δ_{i+1}^{22} \\ 0 & (L̄_{i−1}^{2})_{21}(V_{i})_{12} + (V_{i})_{22} \end{matrix} ] .

Substituting this expression in (33), and looking at the (2, 1) and (2, 2) blocks of the product, we find that for i = 3, …, k − 1,

(36)  0 = (L̄_{i}^{2})_{21}(V_{i+1})_{11} + (V_{i+1})_{21}
(37)  0 = − [(V_{i+1})_{22} + (L̄_{i}^{2})_{21}(V_{i+1})_{12}] Δ_{i+1}^{22} + (L̄_{i}^{2})_{22} [(V_{i})_{22} + (L̄_{i−1}^{2})_{21}(V_{i})_{12}] .

Consideration of the special cases i = 1 and i = 2 shows that (with the convention L̄_{0} = I) the formulae (36) and (37) hold in these cases as well. From (36),

[ (L̄_{i}^{2})_{21} I ] [ \begin{matrix}(V_{i+1})_{11} & (V_{i+1})_{12} \\ (V_{i+1})_{21} & (V_{i+1})_{22} \end{matrix} ] = [ 0 (L̄_{i}^{2})_{21}(V_{i+1})_{12} + (V_{i+1})_{22} ] .

Since [ (L̄_{i}^{2})_{21} I ] has full rank and V_{i+1} is orthogonal, (L̄_{i}^{2})_{21}(V_{i+1})_{12} + (V_{i+1})_{22} is nonsingular. Hence we can define the transformation

y_{a}^{i+1} = [(L̄_{i}^{2})_{21}(V_{i+1})_{12} + (V_{i+1})_{22}]^{−1}z̄_{a}^{i+1}, (38)

for i = 1, …, k − 2, with y_{a}^{1} = z̄_{a}^{1}. Substituting in (28) from (38), we have the following recurrence for i = 2, …, k:

y_{a}^{i} = Δ_{i}^{22}y_{a}^{i−1} + [(L̄_{i−1}^{2})_{21}(V_{i})_{12} + (V_{i})_{22}]^{−1}[(L̄_{i−1}^{2})_{21}z_{b}^{i} + z_{a}^{i}]. (39)

Theorem 4.1 implies that this recurrence is stable. Examination of the other term in (28) leads to

(i−1)z_{b}^{1} = (L̄_{i−1}^{2})_{12}z̄_{a}^{i−1} +^{(i−2)} z_{b}^{1} + (L̄_{i−1}^{2})_{11}z_{b}^{i+1}
= ∑_{j=1}^{i−1}(L̄_{j}^{2})_{12}z̄_{a}^{j} + z_{b}^{1} + ∑_{j=1}^{i−1}(L̄_{j}^{2})_{11}z_{b}^{j+2}.

Since the components of L̄_{j}^{2} are bounded in magnitude by 1, and since the recurrence underlying z̄_{a}^{j} is stable, there is no possibility of ∥^{i}z_{b}^{1}∥ being exponentially larger than ∥z∥.

Having shown that forward substitution involving the factor L is stable, we turn to the back substitution, which involves U.

From (29) and (32), we have that

J =^{def} J^{k} = UD̃Ũ^{−1} ⇔ JŨ = UD̃. (40)

Since J is both block lower Hessenberg and block upper triangular, it must be block upper bidiagonal. We wish to verify that the solution x of the system

Ux = z̃ (41)

does not have ∥x∥/∥z̃∥ exponential in k. By an orthogonal change of variables D̃^{T} x = x̄, (41) becomes

UD̃x̄ = z̃ ⇔ JŨx̄ = z̃. (42)

We can, therefore, prove the desired result by showing that back substitution with both J and Ũ is stable. For the case of Ũ, this is best seen by partitioning the (k+1)n−vectors x̄ and y in the usual way, namely,

x̄ = (x̄_{b}^{1}, x̄_{a}^{1}, x̄_{b}^{2}, x̄_{a}^{2}, …, x̄_{b}^{k+1}, x̄_{a}^{k+1})
y = (y_{b}^{1}, y_{a}^{1}, y_{b}^{2}, y_{a}^{2}, …, y_{b}^{k+1}, y_{a}^{k+1}).

To solve Ũx̄ = y, we first back substitute to find the solution of

˜Ũ [ \begin{matrix} & x̄_{b}^{k} \\  & x̄_{a}^{k} \\  & x̄_{b}^{k+1} \\  & x̄_{a}^{k+1} \end{matrix} ]

and then perform the following recurrence for i = k − 1, k − 2, …, 1:

Ũ_{i} [ \begin{matrix} & x̄_{b}^{i} \\ & x̄_{a}^{i} \end{matrix} ] + Ẽ_{i} [ \begin{matrix} & x̄_{b}^{i+1} \\ & x̄_{a}^{i+1} \end{matrix} ] + G̃_{i} [ \begin{matrix} & x̄_{b}^{k+1} \\ & x̄_{a}^{k+1} \end{matrix} ] = [ \begin{matrix}y_{b}^{i} \\ y_{a}^{i} \end{matrix} ] , i = k − 1, k − 2, …, 1. (43)

Componentwise, (43) is

−Δ_{i+1}^{11}x̄_{b}^{i} − Δ_{i+1}^{12}x̄_{a}^{i} + x̄_{b}^{i+1} = y_{b}^{i}
x̄_{a}^{i} − T_{i}x̄_{a}^{k+1} = y_{a}^{i}.

From Theorem 4.1, T_{i} decays exponentially with i, so the x̄_{a}^{i}, i = 1, 2, …, k − 1, remain bounded. The homogeneous part of the first recurrence is

x̄_{b}^{i} = (Δ_{i+1}^{11})^{−1}x̄_{b}^{i+1}.

Again, we appeal to Theorem 4.1 to deduce that this is stable.

For the back substitution with J, we again refer to the left-hand side of (29) to find expressions for J_{i,i} and J_{i,i+1}. First, from Lemma 4.4 and the fact that each J^{i} is lower Hessenberg, we have for i = 1, …, k −1, that

J_{i,i+1} = J_{i,i+1}^{i} = L̄_{i}^{1}J_{i,i+1}^{i−1} = L̄_{i}^{1}[P_{a}L̄_{i−1}^{2}J_{i−1,i+1}^{i−2} + J_{i,i+1}^{i−2}] = L̄_{i}^{1} [ \begin{matrix}0 & (V_{i+1})_{12} \\ 0 & 0 \end{matrix} ] .

Second,

J_{i,i} = J_{i,i}^{i} = L̄_{i}^{1}J_{i,i}^{i−1} = L̄_{i}^{1}[P_{a}L̄_{i−1}^{2}J_{i−1,i}^{i−2} + J_{i,i}^{i−2}],

and since

J_{i−1,i}^{i−2} = P_{a}L̄_{i−2}^{2}J_{i−2,i}^{i−3} + J_{i−1,i}^{i−3} = [ \begin{matrix}0 & (V_{i})_{12} \\ 0 & 0 \end{matrix} ] ,

we obtain by substitution that

J_{i,i} = L̄_{i}^{1} \{[ \begin{matrix}0 & 0 \\ 0 & (L̄_{i−1}^{2})_{21}(V_{i})_{12} \end{matrix} ] + [ \begin{matrix}(V_{i+1})_{11} & −(V_{i+1})_{12}Δ_{i+1}^{22} \\ 0 & (V_{i})_{22} \end{matrix} ]\} .

(This formula also holds in the special cases i = 1 and i = 2, with a simplified derivation.) Given the system Jy = z̃, partitioning in the usual way, and using (25), we obtain for i = 1, …, k − 1, that

J_{i,i}y^{i} + J_{i,i+1}y^{i+1} = (L̄_{i}^{1})^{−1}z̃^{i}
⇒ [ \begin{matrix}(V_{i+1})_{11} & −(V_{i+1})_{12}Δ_{i+1}^{22} \\ 0 & (L̄_{i−1}^{2})_{21}(V_{i})_{12} + (V_{i})_{22} \end{matrix} ] [ \begin{matrix}y_{b}^{i} \\ y_{a}^{i} \end{matrix} ] + [ \begin{matrix}0 & (V_{i+1})_{12} \\ 0 & 0 \end{matrix} ] [ \begin{matrix}y_{b}^{i+1} \\ y_{a}^{i+1} \end{matrix} ] = [ \begin{matrix} & z̄_{b}^{i} \\  & z̄_{a}^{i} \end{matrix} ] .

It follows immediately that

y_{a}^{i} = [(L̄_{i−1}^{2})_{21}(V_{i})_{12} + (V_{i})_{22}]^{−1}z̄_{a}^{i},

exactly as in (38), and from (39) we know that the sequence of y_{a}^{i}, i = 1, …, k − 1 is governed by a stable recurrence. For the remaining components,

y_{b}^{i} = (V_{i+1})_{11}^{−1}[z̄_{b}^{i} + (V_{i+1})_{12}Δ_{i+1}^{22}y_{a}^{i} − (V_{i+1})_{12}y_{a}^{i+1}].

(V_{i+1})_{11} is clearly nonsingular, since it is a diagonal block in the nonsingular upper triangular matrix J. Recall that for i = 2, 3, …, k, V_{i} = Ō_{i−1}Q_{i}, where Ō_{i−1} is a permutation matrix and Q_{i} is generated as in (17) and (18). It is reasonable to expect that the behavior of the sequence {Q_{i}} is governed by the behavior of some solution to the underlying continuous Lyapunov equation (see, for example, [2, 9]). Therefore, it tends to depend on i/k rather than k. Since, by observation, the permutation matrices O_{i−1}^{2} are also independent of the discretization (they tend to change at indices that correspond to fixed points in the interval [a, b]), we can safely make the following assumption:

Assumption 3. ∥(V_{i+1})_{11}^{−1}∥_{2} ≤ K_{3}, i = 1, …, k, where K_{3} is independent of k.

The end results of the discussion above can be stated as a theorem:

Theorem 4.5. If Assumptions 1, 2 and 3 hold, then the linear system (3) arising from (1)–(2) can be solved in a stable way by using the “serial” version of Algorithm SLU.

Stability of the multipartition versions of SLU can be proved in a similar way. As an aside, we note that the “standard” error analysis for Gaussian elimination factorizations (which would take account of the structure of the coefficient matrix in (3), but not the dichotomy properties of the underlying ODE) actually gives tighter bounds as the number of partitions increases, because the upper bound on element growth during the factorization decreases. For the serial version, element growth of O(2^{(k+1)n}) could in principle occur in the last block column during elimination with row partial pivoting. For the “cyclic reduction” version, this bound is only (2^{n log k}) (i.e., O(k^{n})).

5. Extension to Problems with Multipoint Conditions or Algebraic Parameters. Finally, we describe the extensions of SLU to problems with parameters, or multipoint/integral side conditions. We have already observed that all variants of SLU are simply Gaussian elimination with row partial pivoting, applied to a reordered coefficient matrix; for example, “serial” SLU is obtained by applying row partial pivoting to (20). In this section, we show serial SLU orderings analogous to (20) for these more general problems.

One form of the algebraic-parameter problem [2, p. 322] is

y′ = M(t, λ)y + q(t, λ), t ∈ [a, b], y ∈ R l ^{n}, λ ∈ R l ^{r}
B_{a}y(a, λ) + B_{b}y(b, λ) = d(λ), d ∈ R l ^{(n+r)×n}.

If this problem is linearized with respect to λ, and a finite-difference or multiple shooting method is applied, the following extended form of (3) is obtained:

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

Rearranging rows and columns, we obtain

A = [ \begin{matrix}C_{1} & & & & & & A_{1} & Z_{1} \\ A_{2} & C_{2} & & & & & 0 & Z_{2} \\ & A_{3} & C_{3} & & & & 0 & Z_{3} \\ & & ... & ... & & & ⋮ & ⋮ \\ & & & & C_{k−1} & & 0 & Z_{k−1} \\ & & & & A_{k} & C_{k} & 0 & Z_{k} \\ & & & & & B_{b} & B_{a} & Z_{0} \end{matrix} ] .

We can now proceed with SLU exactly as described in Section 2. The major difference is that the G_{i} blocks will now have n + r columns, rather than n.

The multipoint side condition problem (see [2, pp. 6,323]) is

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

∑_{j=1}^{J} B_{j}y(ξ_{j}) = d, d ∈ R l ^{n},
a = ξ_{1} < ξ_{2} < … < ξ_{J} = b.

When the points ξ_{1}, …, ξ_{J} are included in the finite difference or multiple shooting mesh, the generalization of (3), with some trivial row reordering, is

[ \begin{matrix}A_{1} & C_{1} \\ & A_{2} & C_{2} \\ & & ... & ... \\ & & & & A_{j} & C_{j} \\ & & & & & ... & ... \\ & & & & & & A_{k} & C_{k} \\ B_{1} & 0 & … & 0 & B_{j} & 0 & … & B_{J} \end{matrix} ] [ s_{1} ; s_{2} ; s_{3} ; ⋮ ; s_{k+1} ] = [ f_{1} ; f_{2} ; ⋮ ; f_{k} ; d ] . (45)

Taking the transpose of the coefficient matrix, we obtain

A^{T} = [ \begin{matrix}A_{1}^{T} & & & & & & B_{1}^{T} \\ C_{1}^{T} & A_{2}^{T} & & & & & 0 \\ & C_{2}^{T} & ... & & & & ⋮ \\ & & ... & A_{j}^{T} & & & B_{j}^{T} \\ & & & C_{j}^{T} & ... & & 0 \\ & & & & ... & A_{k}^{T} & ⋮ \\ & & & & & C_{k}^{T} & B_{J}^{T} \end{matrix} ] . (46)

Clearly, (46) differs from (20) only in that the last block column contains some extra nonzeros. Since SLU fills out this last column in any case, this is of little concern. Hence, (45) can be solved by applying SLU to A^{T} rather than A. SLU computes a factorization of the form

P A^{T} Q = LU,

where P and Q are permutation matrices. Since

As = f ⇔ U^{T} L^{T} (P s) = Q^{T} f,

the solution of (45) is obtained by performing a forward substitution with U^{T} and a back substitution with L^{T}.

If the multipoint condition in (45) is replaced by an integral condition of the form

∫_{a}^{b} B(τ)y(τ) dτ = d, (47)

then appropriate discretization of (47) will produce a problem of the form (45), with J = k+1. The resulting technique provides an alternative to the popular approach of introducing artificial dependent variables to obtain a new problem with boundary conditions only.

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, R. M. M. Mattheij, and R. D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, 1988.

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

[4] L. Dieci, M. R. Osborne, and R. D. Russell, A Riccati transformation method for solving linear BVPs. I: Theoretical aspects, SIAM Journal on Numerical Analysis, 25 (1988), pp. 1055–1073.

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

[6] H. B. Keller and M. Lentini, Invariant imbedding, the box scheme and an equivalence between them, SIAM Journal on Numerical Analysis, 19 (1982), pp. 942–962.

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

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

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

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

[11] S. J. Wright, Stable parallel algorithms for two-point boundary value problems, Tech. Rep. MCS–P178–0990, Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, Illinois, September 1990.