A direct solver with O(N) complexity for variable coefficient elliptic PDEs discretized via a high-order composite spectral collocation method

A numerical method for solving elliptic PDEs with variable coefficients on two-dimensional domains is presented. The method is based on high-order composite spectral approximations and is designed for problems with smooth solutions. The resulting system of linear equations is solved using a direct (as opposed to iterative) solver that has optimal O(N) complexity for all stages of the computation when applied to problems with non-oscillatory solutions such as the Laplace and the Stokes equations. Numerical examples demonstrate that the scheme is capable of computing solutions with relative accuracy of $10^{-10}$ or better, even for challenging problems such as highly oscillatory Helmholtz problems and convection-dominated convection diffusion equations. In terms of speed, it is demonstrated that a problem with a non-oscillatory solution that was discretized using $10^{8}$ nodes was solved in 115 minutes on a personal work-station with two quad-core 3.3GHz CPUs. Since the solver is direct, and the"solution operator"fits in RAM, any solves beyond the first are very fast. In the example with $10^{8}$ unknowns, solves require only 30 seconds.

1. Introduction 1.1. Problem formulation. The paper describes a numerical method with optimal O(N ) complexity for solving boundary value problems of the form (1) Au(x) = 0 x ∈ Ω, where Ω is a rectangle in the plane with boundary Γ, and where A is a coercive elliptic partial differential operator (2) [Au](x) = −c 11 The methodology is based on a high order composite spectral discretization and can be modified to handle a range of different domains, including curved ones. For problems with smooth solutions, we demonstrate that the method can easily produce answers with ten or more correct digits.
The proposed method is based on a direct solver which in a single sweep constructs an approximation to the solution operator of (1). This gives the solver several advantages over established linear complexity methods based on iterative solvers (e.g. GMRES or multigrid), perhaps most importantly, the new method can solve problems for which iterative methods converge slowly or not at all. The direct solver has O(N ) complexity for all stages of the computation. A key feature is that once the solution operator has been built, solves can be executed extremely rapidly, making the scheme excel when solving a sequence of equations with the same operator but different boundary data.

1.2.
Outline of solution procedure. The method in this paper is comprised of three steps: (1) The domain is first tessellated into a hierarchical tree of rectangular patches. For each patch on the finest level, a local "solution operator" is built using a dense brute force calculation. The "solution operator" will be defined in Section 1.3; for now we simply note that it encodes all information about the patch that is required to evaluate interactions with other patches.
(2) The larger patches are processed in an upwards pass through the tree, where each parent can be processed once its children have been processed. The processing of a parent node consists of forming its solution operator by "gluing together" the solution operators of its children. (3) Once the solution operators for all patches have been computed, a solution to the PDE can be computed via a downwards pass through the tree. This step is typically very fast.
1.3. Local solution operators. The "local solution operators" introduced in Section 1.2 take the form of discrete approximations to the Dirichlet-to-Neumann, or "DtN," maps. To explain what these maps do, first observe that for a given boundary function f , the BVP (1) has a unique solution u (recall that we assume A to be coercive). For x ∈ Γ, let g(x) = u n (x) denote the normal derivative in the outwards direction of u at x. The process for constructing the function g from f is linear, we write it as g = T f.
From a mathematical perspective, the map T is a slightly unpleasant object; it is a hyper-singular integral operator whose kernel exhibits complicated behavior near the corners of Γ. A key observation is that in the present context, these difficulties can be ignored since we limit attention to functions that are smooth. In a sense, we only need to accurately represent the projection of the "true" operator T onto a space of smooth functions (that in particular do not have any corner singularities).
Concretely, given a square box Ω τ we represent a boundary potential u| Γ and a boundary flux u n | Γ via tabulation at a set of r tabulation points on each side. (For a leaf box, we use r Gaussian nodes.) The DtN operator T τ is then represented simply as a dense matrix T τ of size 4r × 4r that maps tabulated boundary potentials to the corresponding tabulated boundary fluxes.
1.4. Computational complexity. A straight-forward implementation of the direct solver outlined in Sections 1.2 and 1.3 in which all solution operators T τ are treated as general dense matrices has asymptotic complexity O(N 1.5 ) for the "build stage" where the solution operators are constructed, and O(N log N ) complexity for the "solve stage" where the solution operator is applied for a given set of boundary data [16]. This paper demonstrates that by exploiting internal structure in these operators, they can be stored and manipulated efficiently, resulting in optimal O(N ) overall complexity.
To be precise, the internal structure exploited is that the off-diagonal blocks of the dense solution operators can to high accuracy be approximated by matrices of low rank. This property is a result of the fact that for a patch Ω τ , the matrix T τ is a discrete approximation of the continuum DtN operator T τ , which is an integral operator whose kernel is smooth away from the diagonal.
Remark 1.1. The proposed solver can with slight modifications be applied to non-coercive problems such as the Helmholtz equation. If the equation is kept fixed while N is increased, O(N ) complexity is retained. However, in the context of elliptic problems with oscillatory solutions, it is common to scale N to the wave-length so that the number of discretization points per wave-length is fixed as N increases. Our accelerated technique will in this situation lead to a practical speed-up, but will have the same O(N 1.5 ) asymptotic scaling as the basic method that does not use fast operator algebra.
1.5. Prior work. The direct solver outlined in Section 1.2 is an evolution of a sequence of direct solvers for integral equations dating back to [17] and later [9,8,10,4,3]. The common idea is to build a global solution operator by splitting the domain into a hierarchical tree of patches, build a local solution operator for each "leaf" patch, and then build solution operators for larger patches via a hierarchical merge procedure in a sweep over the tree from smaller to larger patches. In the context of integral equations, the "solution operator" is a type of scattering matrix while in the present context, the solution operator is a DtN operator.
The direct solvers [17,9,8,10,4], designed for dense systems, are conceptually related to earlier work on direct solvers for sparse systems arising from finite difference and finite element discretizations of elliptic PDEs such as the classical nested dissection method of George [6,11] and the multifrontal methods by Duff and others [5]. These techniques typically require O(N 1.5 ) operations to construct the LU-factorization of a sparse coefficient matrix arising from the discretization of an elliptic PDE on a planar domain, with the dominant cost being the formation of Schur complements and LU-factorizations of dense matrices of size up to O(N 0.5 ) × O(N 0.5 ). It was in the last several years demonstrated that these dense matrices have internal structure that allows the direct solver to be accelerated to linear or close to linear complexity, see, e.g., [22,7,13,14,19]. These accelerated nested dissection methods are closely related to the fast direct solver presented in this manuscript. An important difference is that the method in the present paper allows high order discretizations to be used without increasing the cost of the direct solver. To be technical, the solvers in [22,7,13,14,19] are based on an underlying finite difference or finite element discretization. High order discretization in this context leads to large frontal matrices (since the "dividers" that partition the grid have to be wide), and consequently very high cost of the LU-factorization.
Our discretization scheme is related to earlier work on spectral collocation methods on composite ("multi-domain") grids, such as, e.g., [12,23], and in particular Pfeiffer et al [18]. For a detailed review of the similarities and differences, see [16].
An O(N 1.5 ) complexity version of the direct solver described in this paper was presented in [16] which in turn is based on [15]. In addition to the improvement in complexity, this paper describes a new representation of the local solution operators that leads to cleaner implementation of the direct solvers and allows greater flexibility in executing the leaf computation, see Remark 3.1.
1.6. Outline of paper. Section 2 introduces the mesh of Gaussian nodes that forms our basic computational grid. Sections 3, 4, and 5 describe a relatively simple direct solver with O(N 1.5 ) complexity. Sections 6, 7, and 8 describe how to improve the asymptotic complexity of the direct solver from O(N 1.5 ) to O(N ) by exploiting internal structure in certain dense matrices. Section 9 describes numerical examples and Section 10 summarizes the key findings.

Discretization
Partition the domain Ω into a collection of square (or possibly rectangular) boxes, called leaf boxes. On the edges of each leaf, place q Gaussian interpolation points. The size of the leaf boxes, and the parameter q should be chosen so that any potential solution u of (1), as well as its first and second derivatives, can be accurately interpolated from their values at these points (q = 21 is often a good choice). Let {x k } N k=1 denote the collection of interpolation points on all boundaries. Next construct a binary tree on the collection of leaf boxes by hierarchically merging them, making sure that all boxes on the same level are roughly of the same size, cf. Figure 1. The boxes should be ordered so that if τ is a parent of a box σ, then τ < σ. We also assume that the root of the tree (i.e. the full box Ω) has index τ = 1. We let Ω τ denote the domain associated with box τ .
With each box τ , we define two index vectors I τ i and I τ e as follows: I τ e A list of all exterior nodes of τ . In other words, k ∈ I τ e iff x k lies on the boundary of Ω τ . I τ i For a parent τ , I τ i is a list of all its interior nodes that are not interior nodes of its children. For a leaf τ , I τ i is empty.
Let u ∈ R N denote a vector holding approximations to the values of u of (1), in other words,  Figure 1. The square domain Ω is split into 4×4 leaf boxes. These are then gathered into a binary tree of successively larger boxes as described in Section 5.1. One possible enumeration of the boxes in the tree is shown, but note that the only restriction is that if box τ is the parent of box σ, then τ < σ.
Finally, let v ∈ R N denote a vector holding approximations to the boundary fluxes of the solution u of (1), in other words , when x j lies on a horizontal edge, , when x j lies on a vertical edge.
Note the v(k) represents an outgoing flux on certain boxes and an incoming flux on others. This is a deliberate choice to avoid problems with signs when matching fluxes of touching boxes.

Constructing the Dirichlet-to-Neumann map for a leaf
This section describes a spectral method for computing a discrete approximation to the DtN map T τ associated with a leaf box Ω τ . In other words, if u is a solution of (1), we seek a matrix T τ of size 4q × 4q such that (3) v(I τ e ) ≈ T τ u(I τ e ). Conceptually, we proceed as follows: Given a vector u(I τ e ) of potential values tabulated on the boundary of Ω τ , form for each side the unique polynomial of degree at most q − 1 that interpolates the q specified values of u. This yields Dirichlet boundary data on Ω τ in the form of four polynomials. Solve the restriction of (1) to Ω τ for the specified boundary data using a spectral method on a local tensor product grid of q × q Chebyshev nodes. The vector v(I τ e ) is obtained by spectral differentiation of the local solution, and then re-tabulating the boundary fluxes to the Gaussian nodes in {x k } k∈I τ e . We give details of the construction in Section 3.2, but as a preliminary step, we first review a classical spectral collocation method for the local solve in Section 3.1 Remark 3.1. Chebyshev nodes are ideal for the leaf computations, and it is in principle also possible to use Chebyshev nodes to represent all boundary-to-boundary "solution operators" such as, e.g., T τ (indeed, this was the approach taken in the first implementation of the proposed method [16]). However, there are at least two substantial benefits to using Gaussian nodes that justify the trouble to retabulate the operators. First, the procedure for merging boundary operators defined for neighboring boxes is much cleaner and involves less bookkeeping since the Gaussian nodes do not include the corner nodes. (Contrast Section 4 of [16] with Section 4.) Second, and more importantly, the use of the Gaussian nodes allows for interpolation between different discretizations. Thus the method can easily be extended to have local refinement when necessary, see Remark 5.2.

Spectral discretization.
Let Ω τ denote a rectangular subset of Ω with boundary Γ τ , and consider the local Dirichlet problem where the elliptic operator A is defined by (2). We will construct an approximate solution to (4) using a classical spectral collocation method described in, e.g., Trefethen [21]: First, pick a small integer p and let {z k } p 2 k=1 denote the nodes in a tensor product grid of p × p Chebyshev nodes on Ω τ . Let D (1) and D (2) denote spectral differentiation matrices corresponding to the operators ∂/∂x 1 and ∂/∂x 2 , respectively. The operator (2) is then locally approximated via the p 2 × p 2 matrix where C 11 is the diagonal matrix with diagonal entries {c 11 (z k )} p 2 k=1 , and the other matrices C ij , C i , C are defined analogously.
Let w ∈ R p 2 denote a vector holding the desired approximate solution of (4). We populate all entries corresponding to boundary nodes with the Dirichlet data from h, and then enforce a spectral collocation condition at the interior nodes. To formalize, let us partition the index set in such a way that J e contains the 4(p − 1) nodes on the boundary of Ω τ , and J i denotes the set of (p − 2) 2 interior nodes, see Figure 2(a). Then partition the vector w into two parts corresponding to internal and exterior nodes via

J Analogously, partition A into four parts via
The potential at the exterior nodes is now given directly from the boundary condition: For the internal nodes, we enforce the PDE (4) via direct collocation: Solving (7) for w i , we find

3.2.
Constructing the approximate DtN. Now that we know how to approximately solve the local Dirichlet problem (4) via a local spectral method, we can build a matrix T τ such that (3) holds to high accuracy. The starting point is a vector u(I τ ) ∈ R 4q of tabulated potential values on the boundary of Ω τ . We will construct the vector v(I τ ) ∈ R 4q via four linear maps. The combination of these maps is the matrix T τ . We henceforth assume that the spectral order of the local Chebyshev grid matches the order of the tabulation on the leaf boundaries so that p = q.
Step 1 -re-tabulation from Gaussian nodes to Chebyshev nodes: For each side of Ω τ , form the unique interpolating polynomial of degree at most q − 1 that interpolates the q potential values on that side specified by u(I τ e ). Now evaluate these polynomials at the boundary nodes of a q × q Chebyshev grid on Ω τ . Observe that for a corner node, we may in the general case get conflicts. For instance, the potential at the south-west corner may get one value from extrapolation of potential values on the south border, and one value from extrapolation of the potential values on the west border. We resolve such conflicts by assigning the corner node the average of the two possibly different values. (In practice, essentially no error occurs since we know that the vector u(I τ e ) tabulates an underlying function that is continuous at the corner.) Step 2 -spectral solve: Step 1 populates the boundary nodes of the q × q Chebyshev grid with Dirichlet data. Now determine the potential at all interior points on the Chebyshev grid by executing a local spectral solve, cf. equation (8).
Step 3 -spectral differentiation: After Step 2, the potential is known at all nodes on the local Chebyshev grid. Now perform spectral differentiation to evaluate approximations to ∂u/∂x 2 for the Chebyshev nodes on the two horizontal sides, and ∂u/∂x 1 for the Chebyshev nodes on the two vertical sides.
Step 4 -re-tabulation from the Chebyshev nodes back to Gaussian nodes: After Step 3, the boundary fluxes on ∂Ω τ are specified by four polynomials of degree q − 1 (specified via tabulation on the relevant Chebyshev nodes). Now simply evaluate these polynomials at the Gaussian nodes on each side to obtain the vector v(I τ e ).
Putting everything together, we find that the matrix T τ is given as a product of four matrices where L i is the linear transform corresponding to " Step i" above. Observe that many of these transforms are far from dense, for instance, L 1 and L 4 are 4 × 4 block matrices with all off-diagonal blocks equal to zero. Exploiting these structures substantially accelerates the computation.
Remark 3.2. The grid of Chebyshev nodes {z k } p 2 j=1 introduced in Section 3.1 is only used for the local computation. In the final solver, there is no need to store potential values at these grid points -they are used merely for constructing the matrix T τ .

Merging two DtN maps
Let τ denote a box in the tree with children α and β. In this section, we demonstrate that if the DtN matrices T α and T β for the children are known, then the DtN matrix T τ can be constructed via a purely local computation which we refer to as a "merge" operation.
We start by introducing some notation: Let Ω τ denote a box with children Ω α and Ω β . For concreteness, let us assume that Ω α and Ω β share a vertical edge as shown in Figure 3, so that We partition the points on ∂Ω α and ∂Ω β into three sets: Figure 3. Notation for the merge operation described in Section 4. The rectangular domain Ω is formed by two squares Ω α and Ω β . The sets J 1 and J 2 form the exterior nodes (black), while J 3 consists of the interior nodes (white).

J 1
Boundary nodes of Ω α that are not boundary nodes of Ω β . J 2 Boundary nodes of Ω β that are not boundary nodes of Ω α . J 3 Boundary nodes of both Ω α and Ω β that are not boundary nodes of the union box Ω τ . Figure 3 illustrates the definitions of the J k 's. Let u denote a solution to (1), with tabulated potential values u and boundary fluxes v, as described in Section 2. Set Recall that T α and T β denote the operators that map values of the potential u on the boundary to values of ∂ n u on the boundaries of the boxes Ω α and Ω β , as described in Section 3. The operators can be partitioned according to the numbering of nodes in Figure 3, resulting in the equations Our objective is now to construct a solution operator S τ and a DtN matrix T τ such that To this end, we write (10) as a single equation: The last equation directly tells us that (11) holds with . By eliminating u 3 from (13) by forming a Schur complement, we also find that (12) holds with (15) T

The full hierarchical scheme
At this point, we know how to construct the DtN operator for a leaf (Section 3), and how to merge two such operators of neighboring patches to form the DtN operator of their union (Section 4). We are ready to describe the full hierarchical scheme for solving the Dirichlet problem (1). This scheme takes the Dirichlet boundary data f , and constructs an approximation to the solution u. The output is a vector u that tabulates approximations to u at the Gaussian nodes {x k } N k=1 on all interior edges that were defined in Section 2. To find u at an arbitrary set of target points in Ω, a post-processing step described in Section 5.3 can be used.
5.1. The algorithm. Partition the domain into a hierarchical tree as described in Section 2. Then execute a "build stage" in which we construct for each box τ the following two matrices: S τ For a parent box τ , S τ is a solution operator that maps values of u on ∂Ω τ to values of u at the interior nodes. In other words, u( (Recall that the index vector I τ e and I τ i were defined in Section 2.) The build stage consists of a single sweep over all nodes in the tree. Any bottom-up ordering in which any parent box is processed after its children can be used. For each leaf box τ , an approximation to the local DtN map T τ is constructed using the procedure described in Section 3. For a parent box τ with children σ 1 and σ 2 , the matrices S τ and T τ are formed from the DtN operators T σ 1 and T σ 2 via the process described in Section 4. Algorithm 1 summarizes the build stage.
Once all the matrices {S τ } τ have been formed, a vector u holding approximations to the solution u of (1) can be constructed for all discretization points by starting at the root box Ω and moving down the tree toward the leaf boxes. The values of u for the points on boundary of Ω can be obtained by tabulating the boundary function f . When any box τ is processed, the value of u is known for all nodes on its boundary (i.e. those listed in I τ e ). The matrix S τ directly maps these values to the values of u on the nodes in the interior of τ (i.e. those listed in I τ i ). When all nodes have been processed, approximations to u have constructed for all tabulation nodes on interior edges. Algorithm 2 summarizes the solve stage.
Remark 5.1. The merge stage is exact when performed in exact arithmetic. The only approximation involved is the approximation of the solution u on a leaf by its interpolating polynomial.
Remark 5.2. To keep the presentation simple, we consider in this paper only the case of a uniform computational grid. Such grids are obviously not well suited to situations where the regularity of the solution changes across the domain. The method described can in principle be modified to handle locally refined grids quite easily. A complication is that the tabulation nodes for two touching boxes will typically not coincide, which requires the introduction of specialized interpolation operators. Efficient refinement strategies also require the development of error indicators that identify the regions where the grid need to be refined. This is work in progress, and will be reported at a later date. We observe that our introduction of Gaussian nodes on the internal boundaries (as opposed to the Chebyshev nodes used in [16]) makes re-interpolation much easier.

Asymptotic complexity.
In this section, we determine the asymptotic complexity of the direct solver. Let N leaf = 4q denote the number of Gaussian nodes on the boundary of a leaf box, and let q 2 denote the number of Chebychev nodes used in the leaf computation. Let L denote the number of levels in the binary tree. This means there are 4 L boxes. Thus the total number of discretization nodes N is approximately 4 L q = (2 L q) 2 q . (To be exact, N = 2 2L+1 q + 2 L+1 q.) Algorithm 1 (build solution operators) This algorithm builds the global Dirichlet-to-Neumann operator for (1). It also builds all matrices S τ required for constructing u at any interior point. It is assumed that if node τ is a parent of node σ, then τ < σ.
Let σ 1 and σ 2 be the children of τ .
Split I σ 1 e and I σ 2 e into vectors I 1 , I 2 , and I 3 as shown in Figure 3.
end if (11) end for Algorithm 2 (solve BVP once solution operator has been built) This program constructs an approximation u to the solution u of (1). It assumes that all matrices S τ have already been constructed in a pre-computation. It is assumed that if node τ is a parent of node σ, then τ < σ. The cost to process one leaf is approximately O(q 6 ). Since there are N q 2 leaf boxes, the total cost of pre-computing approximate DtN operators for all the bottom level is N q 2 × q 6 ∼ N q 4 . Next, consider the cost of constructing the DtN map on level ℓ via the merge operation described in Section 4. For each box on the level ℓ, the operators T τ and S τ are constructed via (14) and (14). These operations involve matrices of size roughly 2 −ℓ N 0.5 × 2 −ℓ N 0.5 . Since there are 4 ℓ boxes per level. The cost on level ℓ of the merge is The total cost for all the merge procedures has complexity L ℓ=1 2 −ℓ N 1.5 ∼ N 1.5 .
Finally, consider the cost of the downwards sweep which solves for the interior unknowns. For any non-leaf box τ on level ℓ, the size of S τ is 2 l q × 2 l (6q) which is approximately ∼ 2 −ℓ N 0.5 × 2 −ℓ N 0.5 . Thus the cost of applying S τ is roughly (2 −ℓ N 0.5 ) 2 = 2 −2ℓ N . So the total cost of the solve step has complexity L−1 l=0 2 2ℓ 2 −2ℓ N ∼ N log N.
In Section 8, we explain how to exploit structure in the matrices T and S to improve the computational cost of both the precomputation and the solve steps.

5.3.
Post-processing. The direct solver in Algorithm 1 constructs approximations to the solution u of (1) at tabulation nodes at all interior edges. Once these are available, it is easy to construct an approximation to u at an arbitrary point. To illustrate the process, suppose that we seek an approximation to u(y), where y is a point located in a leaf τ . We have values of u tabulated at Gaussian nodes on ∂Ω τ . These can easily be re-interpolated to the Chebyshev nodes on ∂Ω τ . Then u can be reconstructed at the interior Chebyshev nodes via the formula (8); observe that the local solution operator −A −1 i,i A i,e was built when the leaf was originally processed and can be simply retrieved from memory (assuming enough memory is available). Once u is tabulated at the Chebyshev grid on Ω τ , it is trivial to interpolate it to y or any other point.

Compressible matrices
The cost of the direct solver given as Algorithm 1 is dominated by the work done at the very top levels; the matrix operations on lines (7) and (8) where N is the total number of discretization nodes, resulting in O(N 1.5 ) overall cost. It turns out that these dense matrices have internal structure that can be exploited to greatly accelerate the matrix algebra. Specifically, the off-diagonal blocks of these matrices are to high precision rank deficient, and the matrices can be represented efficiently using a hierarchical "data-sparse" format known as Hierarchically Block Separable (HBS) (and sometimes Hierarchically Semi-Separable (HSS) matrices [20,1]). This section briefly describes the HBS property, for details see [8].
6.1. Block separable. Let H be an mp × mp matrix that is blocked into p × p blocks, each of size m × m. We say that H is "block separable" with "block-rank" k if for τ = 1, 2, . . . , p, there exist m × k matrices U τ and V τ such that each off-diagonal block H σ,τ of H admits the factorization Observe that the columns of U σ must form a basis for the columns of all off-diagonal blocks in row σ, and analogously, the columns of V τ must form a basis for the rows in all the off-diagonal blocks in column τ . When (16) holds, the matrix H admits a block factorization    6.3. A binary tree structure. The HBS representation of an M × M matrix H is based on a partition of the index vector I = [1, 2, . . . , M ] into a binary tree structure. We let I form the root of the tree, and give it the index 1, I 1 = I. We next split the root into two roughly equi-sized vectors I 2 and I 3 so that I 1 = I 2 ∪ I 3 . The full tree is then formed by continuing to subdivide any interval that holds more than some preset fixed number m of indices. We use the integers ℓ = 0, 1, . . . , L to label the different levels, with 0 denoting the coarsest level. A leaf is a node corresponding to a vector that never got split. For a non-leaf node τ , its children are the two boxes σ 1 and σ 2 such that I τ = I σ 1 ∪ I σ 2 , and τ is then the parent of σ 1 and σ 2 . Two boxes with the same parent are called siblings. These definitions are illustrated in Figure 4.

Definition of the HBS property.
We now define what it means for an M × M matrix H to be hierarchically block separable with respect to a given binary tree T that partitions the index vector J = [1, 2, . . . , M ]. For simplicity, we suppose that for every leaf node τ the index vector I τ holds precisely m points, so that M = m 2 L . Then H is HBS with block rank k if the following two conditions hold: (1) Assumption on ranks of off-diagonal blocks at the finest level: For any two distinct leaf nodes τ and τ ′ , define the m × m matrix (18) H τ,τ ′ = H(I τ , I τ ′ ).
Then there must exist matrices U τ , V τ ′ , andH τ,τ ′ such that (2) Assumption on ranks of off-diagonal blocks on level ℓ = L − 1, L − 2, . . . , 1: The rank assumption at level ℓ is defined in terms of the blocks constructed on the next finer level ℓ + 1: For any distinct nodes τ and τ ′ on level ℓ with children σ 1 , σ 2 and σ ′ 1 , σ ′ 2 , respectively, define  Then there must exist matrices U τ , V τ ′ , andH τ,τ ′ such that An HBS matrix is now fully described if the basis matrices U τ and V τ are provided for each node τ , and in addition, we are for each leaf τ given the m × m matrix (22) D τ = H(I τ , I τ ), and for each parent node τ with children σ 1 and σ 2 we are given the 2k × 2k matrix Observe in particular that the matricesH σ 1 ,σ 2 are only required when {σ 1 , σ 2 } forms a sibling pair. Figure 5 summarizes the required matrices.
Furthermore, we letH (ℓ) denote the block matrix whose diagonal blocks are zero, and whose offdiagonal blocks are the blocksH τ,τ ′ for all distinct τ, τ ′ on level ℓ. With these definitions, and finally

Fast arithmetic operations on HBS matrices
Arithmetic operations involving dense HBS matrices of size M ×M can often be executed in O(M ) operations. This fast matrix algebra is vital for achieving linear complexity in our direct solver. This section provides a brief introduction to the HBS matrix algebra. We describe the operations we need (inversion, addition, and low-rank update) in some detail for the single level "block separable" format. The generalization to the multi-level "hierarchically block separable" format is briefly described for the case of matrix inversion. A full description of all algorithms required is given in [7], which is related to the earlier work [2].
Before we start, we recall that a block separable matrix H consisting of p × p blocks, each of size m × m, and with "HBS-rank" k < m, admits the factorization Inversion of a block separable matrix. The decomposition (31) represents H as a sum of one term UHV * that is "low rank," and one term D that is easily invertible (since it is block diagonal). By modifying the classical Woodbury formula for inversion of a matrix perturbed by the addition of a low-rank term, it can be shown that (see Lemma 3.1 of [8]) assuming the inverses in formulas (32) -(36) all exist. Now observe that the matricesD, E, F, and G can all easily be computed since the formulas defining them involve only block-diagonal matrices. In consequence, (32) reduces the task of inverting the big (size mp × mp) matrix H to the task of inverting the small (size kp × kp) matrixH +D.
When H is not only "block separable", but "hierarchically block separable", the process can be repeated recursively by exploiting thatH +D is itself amenable to accelerated inversion, etc. The resulting process is somewhat tricky to analyze, but leads to very clean codes. To illustrate, we include Algorithm 3 which shows the multi-level O(M ) inversion algorithm for an HBS matrix H. The algorithm takes as input the factors {U τ , V τ , D τ , B τ } τ representing H (cf. Figure 5), and outputs an analogous set of factors {E τ , F τ , G τ } τ representing H −1 . With these factors, the matrix-vector multiplication y = H −1 x can be executed via the procedure described in Algorithm 4.

Addition of two block separable matrices. Let H A and H B be block separable matrices with factorizations
and Then H = H A + H B can be written in block separable form via To restore (37) to block separable form, permute the rows and columns of U A U B and V A V B to attain block diagonal form, then re-orthogonalize the diagonal blocks. This process in principle results in a matrix H whose HBS-rank is the sum of the HBS-ranks of H A and H B . In practice, this Algorithm 3 (inversion of an HBS matrix) Given factors {U τ , V τ , D τ , B τ } τ representing an HBS matrix H, this algorithm constructs factors loop over all levels, finer to coarser, ℓ = L, L − 1, . . . , 1 loop over all boxes τ on level ℓ, if τ is a leaf nodẽ D τ = D τ else Let σ 1 and σ 2 denote the children of τ .
Algorithm 4 (application of the inverse of an HBS matrix) Given x, compute y = H −1 x using the factors {E τ , F τ , G τ } τ resulting from Algorithm 3. loop over all leaf boxes τ x τ = F * τ x(I τ ). end loop loop over all levels, finer to coarser, ℓ = L, L − 1, . . . , 1 loop over all parent boxes τ on level ℓ, Let σ 1 and σ 2 denote the children of τ .
end loop end loop loop over all leaf boxes τ y(I τ ) = E τqτ + G τ x(I τ ). end loop rank increase can be combated by numerically recompressing the basis matrices, and updating the middle factor as needed. For details, as well as the extension to a multi-level scheme, see [2,7]. 7.3. Addition of a block separable matrix with a low rank matrix. Let H B = QR be a k-rank matrix where Q and R * are of size mp × k. We would like to add H B to the block separable matrix H A . Since we already know how to add two block separable matrices, we choose to rewrite H B in block separable form. Without loss of generality, assume Q is orthogonal. Partition Q into p blocks of size m × k. The blocks make up the matrix U B . Likewise partition R into p blocks of size k × m. The block matrix D B has entries D τ = Q τ R τ for τ = 1, . . . , p. To construct the matrices V B , for each τ = 1, . . . , p, the matrix R τ is factorized intoR τ V τ * where the matrix V τ is orthogonal. The matricesR τ make up the entries ofH B .

Accelerating the direct solver
This section describes how the fast matrix algebra described in Sections 6 and 7 can be used to accelerate the direct solver of Section 5 to attain O(N ) complexity. We recall that the O(N 1.5 ) cost of Algorithm 1 relates to the execution of lines (7) and (8) at the top levels, since these involve dense matrix algebra of matrices of size O(N 0.5 ) × O(N 0.5 ). The principal claims of this section are: • The matrices T σ 1 1,3 , T σ 1 3,1 , T σ 2 2,3 , T σ 2 3,2 have low numerical rank. • The matrices T σ 1 1,1 , T σ 2 2,2 , T σ 1 3,3 , T σ 2 3,3 are HBS matrices of low HBS rank. To be precise, the ranks that we claim are "low" scale as log(1/ν) × log(m) where m is the number of points along the boundary of Ω τ , and ν is the computational tolerance. In practice, we found that for problems with non-oscillatory solutions, the ranks are extremely modest: when ν = 10 −10 , the ranks range between 10 and 80, even for very large problems.
The cause of the rank deficiencies is that the matrix T τ is a highly accurate approximation to the Dirichlet-to-Neumann operator on Ω τ . This operator is known to have a smooth kernel that is nonoscillatory whenever the underlying PDE has non-oscillatory solutions. Since the domain boundary ∂Ω τ is one-dimensional, this makes the expectation that the off-diagonal blocks have low rank very natural, see [8]. It is backed up by extensive numerical experiments (see Section 9), but we do not at this point have rigorous proofs to support the claim.
Once it is observed that all matrices in lines (7) and (8) of Algorithm 1 are structured, it becomes obvious how to accelerate the algorithm. For instance, line (7) is executed in three steps: (i) Add the HBS matrices T σ 1 3,3 and −T σ 2 3,3 . (ii) Invert the sum of the HBS matrices. (iii) Apply the inverse (in HBS form) to one of the low rank factors of −T α 3,1 T β 3,2 . The result is an approximation to S τ , represented as a product of two thin matrices. Executing line (8) is analogous: First form the matrix products T σ 1 1,3 S τ and T σ 2 2,3 S τ , exploiting that all factors are of low rank. Then perform a low-rank update to a block-diagonal matrix whose blocks are provided in HBS-form to construct the new HBS matrix T τ .
Accelerating the solve stage in Algorithm 2 is trivial, simply exploit that the matrix S τ on line (3) has low numerical rank.
Remark 8.1. Some of the structured matrix operators (e.g. adding two HBS matrices, or the lowrank update) can algebraically lead to a large increase in the HBS ranks. We know for physical reasons that the output should have rank-structure very similar to the input, however, and we combat the rank-increase by frequently recompressing the output of each operation.
Remark 8.2. In practice, we implement the algorithm to use dense matrix algebra at the finer levels, where all the DtN matrices T τ are small. Once they get large enough that HBS algebra outperforms dense operations, we compress the dense matrices by brute force, and rely on HBS algebra in the remainder of the upwards sweep.

Numerical Examples
In this section, we illustrate the capabilities of the proposed method for the boundary value problem where Ω = [0, 1] 2 , Γ = ∂Ω, and c 1 (x), c 2 (x), and c(x) are smooth, cf. (2). The choice of the functions c 1 (x), c 2 (x), and c(x) will vary for each example. All experiments reported in this section were executed on a machine with two quad-core Intel Xeon E5-2643 processors and 128GB of RAM. The direct solver was implemented in Matlab, which means that the speeds reported can very likely be improved, although the asymptotic complexity should be unaffected.
In Section 9.1 we apply the direct solver to four problems with known analytic solutions. This allows us to very accurately investigate the errors incurred, but is arguably a particularly favorable environment. Section 9.2 presents results from more general situations where the exact solution is not known, and errors have to be estimated.
In all experiments, the number of Gaussian points per leaf edge q is fixed at 21, and 21 × 21 tensor product grids of Chebyshev points are used in the leaf computations. Per Remark 8.2, we switch from dense computations to HBS when a box has more than 2000 points on the boundary. In a first experiment, we prescribed the tolerance in the "fast" matrix algebra to be ǫ = 10 −7 .
where u app denotes the approximate solution constructed by the direct solver. Our expectation is for all problems except Helmholtz III to exhibit optimal linear scaling for both the build and the solve stages. Additionally, we expect the cost of applying the Dirichlet-to-Neumann operator T 1 to scale as N 0.5 for all problems except Helmholtz III. The numerical results clearly bear this out for Laplace and Helmholtz I. For Helmholtz II, it appears that linear scaling has not quite taken hold for the range of problem sizes our hardware could manage. The Helmholtz III problem clearly does not exhibit linear scaling, but has not yet reached its O(N 1.5 ) asymptotic at the largest problem considered, which was of size roughly 426 × 426 wave-lengths. We remark that the cost of the solve stage is tiny. For example, a problem involving 11 million unknowns (corresponding to approximately 100 million discretization points) takes 115 minutes for the build stage and then only 30 seconds for each additional solve. The cost for applying the Dirichlet-to-Neumann operator is even less at 0.36 seconds per boundary condition. Figure 6 illustrates the scaling via log-log plots.
In a second set of experiments, we investigated the accuracy of the computed solutions, and in particular how the accuracy depends on the tolerance ǫ in the fast matrix algebra. In addition to reporting E pot , Table 2 reports where u app denotes the approximate solution constructed by the direct solver for tolerances ε = 10 −7 , 10 −10 , and 10 −12 . The number of discretization points was fixed problem to be N = 693504 (N tot = 7252225).  Table 2. Errors for solving the PDEs in Section 9.1 with different user prescribed tolerances when the number of discretization points is fixed at N = 693504 (N tot = 7252225).
The solution obtains (or nearly obtains) the prescribed tolerance while the normal derivative suffers from a three digit loss in accuracy. This loss is likely attributable to the unboundedness of the Dirichlet-to-Neumann map. The compressed representation captures the high end of the spectrum to the desired accuracy while the low end of the spectrum is captured to three digits less than the desired accuracy.  Figure 6. The first three graphs give the times required for building the direct solver (T build ), solving a problem (T solve ) and applying the approximate Dirichlet-to-Neumann operator on ∂Ω (T apply ). The fourth graph gives the memory R in MB required to store the solution operator.
For the three experiments, the boundary data is given by f (x) = cos(2x 1 )(1 − 2x 2 ).  Table 4. Times in seconds for solving the PDE's in Section 9.2 with a user subscribed tolerance of ǫ = 10 −12 .
To check for convergence, we post-process the solution as described in Section 5.3 to get the solution on the Chebyshev grid. Let u N denote the solution on the Chebyshev grid. Likewise,let u N n denote the normal derivative on the boundary at the Chebyshev boundary points. We compare the solution and the normal derivative on the boundary pointwise at the locationŝ The tolerance for the compressed representations is set to ǫ = 10 −12 . Table 3 reports the pointwise errors. We see that high accuracy is obtained in all cases, with solutions that have ten correct digits for the potential and about seven correct digits for the boundary flux.
The computational costs of the computations are reported in Table 4. The memory R reported now includes the memory required to store all local solution operators described in Section 5.3.

Conclusions
We have described a direct solver for variable coefficient elliptic PDEs in the plane, under the assumption that the solution and all coefficient functions are smooth. For problems with nonoscillatory solutions such as the Laplace and Stokes equations, the asymptotic complexity of the solver is O(N ), with a small constant of proportionality. For problems with oscillatory solutions, high practical efficiency is retained for problems of size up to several hundred wave-lengths.
Our method is based on a composite spectral discretization. We use high order local meshes (typically of size 21 × 21) capable of solving even very large scale problems to ten correct digits or more. The direct solver is conceptually similar to the classical nested dissection method [6]. To improve the asymptotic complexity from the classical O(N 1.5 ) to O(N ), we exploit internal structure in the dense "frontal matrices" at the top levels in a manner similar to recent work such as, e.g., [2,7,13,14,19]. Compared to these techniques, our method has an advantage in that high order discretizations can be incorporated without compromising the speed of the linear algebra. The reason is that we use a formulation based on Dirichlet-to-Neumann operators. As a result, we need high order convergence only in the tangential direction on patch boundaries.
The direct solver presented requires more storage than classical iterative methods, but this is partially off-set by the use of high-order discretizations. More importantly, the solver is characterized by very low data-movement. This makes the method particularly well suited for implementation on parallel machines with distributed memory.