text:
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.


The Study and Application of Hadoop across Multiple Clusters

Hadoop is a widely applied tool for large-scale data-intensive computing in big data, but it can only be implemented on single cluster environment. In this paper, we focus on the application of Hadoop across multiple clusters and dedicate to solve the key problems of data sharing and task scheduling among clusters. A hierarchical distributed computing architecture of Hadoop across multiple clusters is designed. The virtual HDFS and job adapter are proposed to provide global data view and task allocation across multiple data centers. The job submitted by user to this platform is decomposed automatically into several sub-jobs and then allocated to corresponding cluster by location-aware manner. A prototype based on this architecture is presented and currently applied in the distributed spatial information processing across spatial data centers. Keywords-multiple clusters; Apache Hadoop; hierarchical distributed computing; virtual HDFS; job adapter

INTRODUCTION
The MapReduce paradigm has emerged as a highly successful programming model for large-scale dataintensive computing applications [1]. However, current Hadoop implementations cannot be applied to distributed data processing across data centers. Nowadays, dataintensive computing typically uses modern data center architectures and massive data processing paradigms [2]. The requirements for data-intensive analysis of scientific data across distributed clusters or data centers have grown observably in recent years. This research is devoted to a study on the distributed data processing model across multiple data centers.
In this paper, we try to improve the capability and flexibility of Hadoop. A new hierarchical distributed computing architecture of Hadoop across multiple data centers is proposed. Java Socket is used to transfer these tasks to the corresponding cluster and data center, where the data processing in supposed to be run without data moving. Virtual HDFS (Hadoop Distributed File System) is presented to provide the catalog service of these data centers, which can record the meta-data and access path of all data. The design and application of this architecture are illustrated in this paper and we hope this study may enlighten the attention of Hadoop global implementation.
The rest of this paper is organized as follows: Section 2 discusses related works of our research, Section 3 presents the design of a hierarchical distributed computing architecture of Hadoop, and a prototype system based on this architecture is shown in Section 3. Finally, Section 4 concludes the paper and points out the future work.
II. RELATED WORKS Cloud computing is a set of network enabled services to allow the centralized data storage and online access to computer services or resources [3]. It provides scalable, quality guaranteed, normally personalized, convenient computing infrastructure on demand. The MapReduce paradigm and its open sourced implementation-Hadoop have been recognized as a representative enabling technique for Cloud computing [4].
In distributed computing frameworks based on Hadoop, one job committed to master node (name node) may be divided into several same tasks and then run on several slave nodes (data nodes). The data of this job is stored in HDFS (Hadoop Distributed File System) which is distributed to data nodes of the same cluster. In the situation of distributed computing crossing clusters, the data copying from multiple sites to the computing center is the traditional method.
But the copy process is tedious and inefficient between global distributed data centers through wide-area network, especially in era of big data. Moving the computation instead of moving the data is the proper way to tackle this problem [5]. By using data parallel processing paradigms on multiple clusters, simulations can be run on multiple computing centers concurrently without the need of copying the data. In G-Hadoop [6], Gfarm and Torque were adopted to manage data files and resources for clusters. In P2P-MapReduce [7], the adaptive MapReduce International Conference on Mechatronics, Electronic, Industrial and Control Engineering (MEIC 2014) framework provided a more reliable MapReduce middleware in dynamic Cloud infrastructures. A Federated MapReduce [8] provided a transparent way to run original MapReduce jobs across multiple clusters without any extra programming burden.
The Apache Hadoop MapReduce implementation may be upgraded for a multi-cluster environment with a decision algorithm that would prefer local computers to the remote [9]. The research on Hadoop across data centers or clusters is rare, but we believe this issue is worthwhile. In our works, we try to employ Hadoop in the application of distributed spatial data processing with high performance.

III. THE DESIGN OF HIERARCHICAL DISTRIBUTED
COMPUTING ARCHITECTURE OF HADOOP With the rapid development of manufacture technology and observation technique, human can obtain massive spatial data of variant types from multiple sources. Researchers have devoted to collecting variant spatial data and providing spatial information service. Many spatial data centers have been set up all over the world for different purposes and diverse services. The processing and managing of spatial data must face the status of multiple data centers.
In the design and implementation of distributed computing architecture of Hadoop across multiple spatial data centers, there are mainly three problems need to be discussed and tackled. (1) How to obtain the view of all data in multiple data centers, which may decide the allocation of computing task to proper data center, (2) How to decompose the job into tasks and dispense them to corresponding clusters, which can synergistically control each Hadoop execution engine in multiple clusters, (3) How to deliver tasks transparently across different administrative domains, which need efficient and general way to keep communication among different clusters. In this paper, we try to employ the hierarchical methodology [10] and propose a hierarchical distributed computing architecture based on Hadoop across multiple clusters (named HDC-Hadoop). The overview of this architecture is shown as Fig .1. The HDC-Hadoop architecture represents a master/slave communication model. The master node is the central entity in this architecture and it is responsible for accepting jobs submitted by user, splitting the jobs into smaller sub-jobs and distributing these sub-jobs to the certain slave nodes where required data are located. The master node is also in charge of managing the metadata of all files available in every data nodes (virtual HDFS view of multiple HDFS). The slave node is installed on each participating cluster and enables it to run sub-jobs allocated by Sub-Job Adapter in the master node (we can also call it global name node).
The virtual HDFS provides a global view of all data files in HDFSs of every cluster. In HDC-Hadoop data must be managed in a location-aware manner in order to provide the required location information for the sub-job adapter on master node. The virtual HDFS just records meta-data of all HDFS files obtained from name nodes, such as data size, file format, access path and so on. When new file is input or existing file is modified, related name node (records the meta-data of HDFS in one cluster) must send the updating meta-data to the virtual HDFS. The relationship of virtual HDFS and HDFSs is shown in Fig .2.  Data Node local file system executed step by step in different clusters (so as to say, the program is organized into several sequential sub-jobs of different clusters).

The whole workflow of a job execution in HDC-Hadoop is composed of five steps:
Job submission: user submits a MapReduce program and configuration (including parameters, input files and additional resources) to master node of HDC-Hadoop. The job scheduler may set a job ID for the new job and push it in the queue of jobs.
Program compiling and decomposing: job in queue may be split into small location-aware sub-jobs based on required data. Based on meta-data in virtual HDFS, if all files in the MapReduce program are located in one cluster (data center), the program may be transferred directly to corresponding cluster (slave node). Else if the files in the MapReduce program are located in different clusters, the program may be recompiled to Map-Map-Reduce or Map-Reduce-Map-Reduce program and organized into several sequential sub-jobs of different clusters. The file names in Map and Reduce functions based on virtual HDFS are translated into local file path of the HDFS in clusters.
Sub-job localization and assignment: based on virtual HDFS, sub-job may be assigned properly to the name node of related cluster in data-aware manner. This allocation adopts socket communication method, which can inter-transfer message, program and file based on TCP (Transmission Control Protocol). And then a job execution (the method runJob() of Hadoop is called) can be run on the cluster based on local file system.
Task submission and Task execution: after the job is localized on the cluster, the JobTracker in name node splits into smaller tasks. When the TaskTraker in data node receives a new task, it localizes the tasks executable and resources in local file system. Then the job is executed by spawning a new JVM (Java Virtual Machine) on the computed node and running the corresponding task with the configured parameters.
Result feedback: after all tasks of a sub-job are completed in one cluster, the name node may return the result files paths to the sub-job adapter in master node. The adapter checks whether this sub-job is a simple job (all required files located in one cluster, that is to say there is only one sub-job corresponding to a user job). If it is a simple job then register the new files in virtual HDFS and return the file name list to the job scheduler. Else if this sub-job has other subsequent or relevant sub-jobs, the output of this sub-job may be used as the input of other sub-jobs which may then be allocated and executed on corresponding cluster. Till all the sub-jobs of one user job are completed, the output files are registered in virtual HDFS and the file name list is returned to job scheduler. Finally, job scheduler returns the results to corresponding user by job ID.
Based on the HDC-Hadoop architecture, distributed computing job across multiple clusters or data centers can be implemented.

IV. APPLICATION AND PROTOTYPE
This section discusses a prototype system of this HDC-Hadoop architecture. In this prototype, two clusters are built on two data centers separately based on Hadoop. The schematic diagram of the deployment is shown in Fig .3. In Figure 3, there are two distributed data centers (two data centers for remote sensing raster data, one for MODIS data service and another for ETM data service). On the data storage servers, we installed virtual machines as data nodes of Hadoop and set the disks sharing, which can make the data accessible in Hadoop.
Moreover, in each data center, we added a server as name node of Hadoop system, which takes charge of the sub-job receiving from master node and task assigning to data nodes. In the master node, Web portal provide the interface for users to submit job and obtain results, job center is responsible for the job scheduling and managing, and MDDB (meta-data database) stores the meta-data of virtual HDFS.
When user submits a new job, he firstly searches the required data based on the data view of virtual HDFS, which can make the details of multiple data centers transparent, shown as Fig .4. Master Node machine, which may be transferred to the cluster where the required data is stored. In the job center, the execution of each job can be monitored in Web portal. The workflow overview of job execution is shown as Fig .5. The final execution results are expressed in form of KML (Keyhole Markup Language) files, and the access paths may be returned to the user interface in Web portal server. Based on the HDC-Hadoop architecture, the simulation system can realize data distributed computing based on Hadoop across multiple clusters, but this prototype is still a test bed currently and there are many unresolved problems yet.

V. CONCLUSIONS
The requirements for data-intensive computing across distributed clusters and multiple data centers have grown significantly in recent years. However, original MapReduce paradigm and Hadoop platform are developed to operate on single cluster environments and cannot be implemented in large-scale distributed data processing across multiple cluster and data centers. The goal of this research is to apply Hadoop framework in the large-scale distributed computing across multiple clusters. In this paper, a hierarchical distributed computing architecture is designed and presented. This architecture is based on master/slave communication model and virtual HDFS is proposed to provide global view of data sets across distributed clusters. The HDC-Hadoop architecture, virtual HDFS structure and job execution flowchart are illuminated, and also a prototype based on this architecture is built up and presented. The implementation results show that this design is feasible and has application prospect. To make HDC-Hadoop fully functional and implementable, next step we plan to enhance the distributed file system across wide area networks and design the security mechanism across multiple administrative domains for this architecture.


extfieldsofstudy:
Computer Science
Mathematics

text:
Modeling European Options

Option contracts can be valued by using the Black-Scholes equation, a partial differential equation with initial conditions. An exact solution for European style options is known. The computation time and the error need to be minimized simultaneously. In this paper, the authors have solved the Black-Scholes equation by employing a reasonably accurate implicit method. Options with known analytic solutions have been evaluated. Furthermore, an overall second order accurate space and time discretization has been accomplished in this paper.

price in response to external effects. It is represented by a random sample drawn from a normal distribution with mean zero denoted as dX  , Putting the above two factors together we obtain a stochastic differential equation, which is the mathematical representation for generating asset prices [1,4].
The equation (1) gives interesting and important information concerning the behavior of S in a probabilistic sense. At time ' t t  , suppose the price is ' S , and then ' S will be distributed about 0 S with a probability density function. The future asset price ' S is thus most likely to be close to 0 S and less likely to be far away. Thus, the equation generates time serieseach time the series is restarted a different path results. Each path is called a realization of the random walk.
The Black-Scholes Partial Differential Equation, the derivation of which is well-known, is given below If E S  at expiry, then in financial sense, for call option, there will be profit E S  , handing over an amount E , to obtain an asset worth, while, if E S  , then there will be loss of S E  . Thus the value of the call option at expiry can be written as As the time tends to expiry date the value of call option approaches (6), it is known as pay-off function for European Call Option. This is known as the final condition of PDE (12).
which means pay-off is also zero. Thus, the call option is worthless on 0  S even if there is long time to expiry. Hence, we have i.e. the asset price increases without bound it becomes ever more likely that the option will be exercised and the magnitude of the exercise price becomes less and less important.
And the initial condition   t P , 0 is determined by calculating the present value of an amount E received at timeT . For time-dependent interest rate, the boundary condition at And, Thus, the Black Scholes equation and boundary conditions for a European Put Option are given by equations (2), (6)-(8).
The analytical solution to the BSDE is not the focal point of the paper but given below, nonetheless as [2] The corresponding calculation for a European Call Option follows similar lines.

Analytical Solution:
for some constants  and  , subsequent differentiation yields , so that the PDE assumes the standard parabolic form The argument follows from the fact that Next, evaluate the inner integral, The corresponding calculation for a European Put Option follows similar lines. Its transformed pay-off is and can be computed as above. However, having evaluated the Call, a simpler way is to use the put-call parity formula for the value P of a Put, given the value of the Call. This yields

Definitions/Terminologies:
Writer: The person selling the asset.
Holder: The person buying the asset.
Expiry Date: At a prescribed time in the future, the holder/writer of the option may purchase/sell a prescribed asset. Underlying Asset: The holder/writer may purchase/sell a prescribed asset at a prescribed time.
Exercise Amount/Strike Price: The holder/writer purchases a prescribed asset at a prescribed time at a prescribed amount.
Arbitrage: In financial terms, there are never any opportunities to make an instantaneous risk-free profit.
Risk: Risk is commonly of two typesspecific and non-specific, specific risk is the component of risk associated with a single asset, whereas nonspecific risk is associated with factors affecting the whole market.

Notation table:
t -Time

Software System Attributes
Dept. of CSE, PES Institute of technology, Bangalore South Campus  Reliability: This implementation is based on the Black-Scholes and Crank Nicolson's Equations; it is mathematically sound with no error included during coding. This does not slow down the system or does not cause any harm to the environment.
Availability: This software is built with a simple user interface and easily available resources.
Portability: Since this product uses platform independent HTML and C codes, it is portable. Assign to y maximum value between z-K and 0, which symbolizes the call option else Assign to y maximum value between K-z and 0, which symbolizes the put option Step3: stop


The Effect of Flashcache and Bcache on I/O Performance

Solid state drives (SSDs) provide significant improvements in random I/O performance over traditional rotating SATA and SAS drives. While the cost of SSDs has been steadily declining over the past few years, high density SSDs continue to remain prohibitively expensive when compared to traditional drives. Currently, 1 TB SSDs generally cost more than USD $1,000, while 1 TB SATA drives typically retail for under USD $100. With ever-increasing x86_64 server CPU core counts, and therefore job slot counts, local scratch space density and random I/O performance have become even more important for HEP/NP applications. Flashcache and Bcache are Linux kernel modules which implement caching of SATA/SAS hard drive data on SSDs, effectively allowing one to create hybrid SSD drives using software. In this paper, we discuss our experience with Flashcache and Bcache, and the effects of this software on local scratch storage performance.

Introduction
Over the past few years, x86 64 server processor core counts have continued to increase. Sandy Bridge processors used in servers today already commonly provide 16 logical cores per physical CPU. In the not so distant future, Haswell-based server CPUs are expected to provide up to 40 logical cores in a single physical package [1]. In the traditional High Energy Physics and Nuclear Physics (HEP/NP) batch processing model, one job slot is allocated for each logical core in a system. This has lead to a constant increase in the number of jobs per host, and therefore in demand for local scratch and remote storage I/O performance and density. There has been some success in the development of multithreaded HEP/NP batch processing applications and frameworks, which may help alleviate increasing random I/O performance requirements, potentially through the use of scatter/gather I/O models [2]. However, as of the time of writing, multithreaded batch processing has not yet been widely adopted by the experiments utilizing our computing center, the RHIC/ATLAS Computing Facility (RACF) at Brookhaven National Laboratory (BNL).
The current method we've adopted for improving scratch I/O performance on our worker nodes has been to populate these systems with multiple commodity drives, and combine them into software RAID0 arrays. This can be somewhat expensive, however, as it requires that many drives be purchased for each host. Furthermore, the increase in core density is outpacing the increase in available spindle count in typical worker node server hardware, and therefore, this solution may not scale at some point in the near future. SSDs provide excellent random I/O performance characteristics, and we considered adopting them for scratch storage in our worker nodes. However, they are extremely expensive for the capacities we require. For instance, as of September 2013, 1 TB SSD drives typically cost near or above USD $1,000, with enterprise models exceeding USD $2,500. In contrast, a commodity 1 TB SATA drive retails for less than USD $100. Therefore, we decided to look into the possibility of using hybrid hard drives, which have a flash memory cache in front of traditional rotating media. While Seagate offers such products close to the price of traditional SATA drives, they have not yet been validated for use in servers offered by large scale manufacturers such as Dell. Therefore, this technology was not a viable option for us.
Hybrid devices can also be implemented in software, and there are currently a number of Linux kernel drivers available which implement this functionality, two of which are Flashcache and Bcache. We tested both of these driver implementations to determine if they could be used to eliminate the need for multi-spindle RAID0 array scratch storage on worker nodes, or increase the performance of these arrays.

Flashcache and Bcache
While Flashcache and Bcache provide similar functionality, there are a number of important differences between the two projects. Bcache has been integrated into the Linux kernel as stable software since the 3.10 release. Flashcache, on the other hand, is available only as separate software, which is compiled as a module outside of the kernel source. Whereas Flashcache is controlled via the kernel's sysctl interface, Bcache uses Sysfs (/sys) for configuration parameter modification. Finally, Bcache utilizes a btree cache structure [3], while Flashcache's cache is structured as a set associative hash [4]. The steps we used to setup and configure our Flashcache and Bcache devices are available in Table 1.  There are number of configuration options available for both drivers. Most importantly, all of our tests were performed with the Flashcache and Bcache device SSD caches in writeback mode. A listing of some additional relevant parameter settings for our tests are present in Table  2. The number of available options is quite large: testing variations on each was not possible given our time constraints. For simplicity, in general, configuration defaults were untouched, and only altered when it was believed that significant performance gains could be achieved, or when such modifications were recommended in documentation.

Hardware and Operating System Configuration
Details on the hardware configuration used during our benchmarking is available in Table 3. 64-bit Scientific Linux (SL) 6.4 was used as the operating system in all tests. However, the use of a custom configured/compiled vanilla upstream 3.11.1 kernel was necessary during Bcache testing. This was because it was not possible (at least not without heavy modification) to extract the Bcache source from the upstream 3.x kernel tree, and integrate it into the stock SL6 2.6.32 kernel source.
For the Flashcache and Bcache benchmarks, two hardware configurations were tested: one with a single SSD cache in front of a single SATA drive, and one with a single SSD cache fronting a 7-spindle software RAID0 array. For comparison, the same benchmarks were run on a single SSD, a single SATA drive, and an 8-spindle software RAID0 device. Due to the fact that SSDs are not subject to the I/O timings required for efficient access to rotating media, the "deadline" or "noop" kernel I/O schedulers should be used with these drives. During our tests, we set the kernel I/O scheduler to "deadline" for the SSD, and TRIM was not enabled at mount time.
Only one SSD model was tested in our benchmarks due to time and available hardware constraints. It's possible that better performance could have been obtained from Flashcache and Bcache through the use of a higher performing SSD, or multiples thereof. However, the performance of the SSD utilized in this study was in line with the results of other SSD makes/models we've benchmarked in the past.

Benchmarks
For performance testing, we utilized both Bonnie++ [5] and Iozone [6]. When these benchmarks were run, we were careful to ensure that filesize was set larger than total system memory size. Bonnie++ only tests sequential I/O performance, but it supports a synchronized multiprocess/parallel execution mode. In this way, the benchmark can be used to create a random workload, albeit not heavily randomized due to kernel I/O buffering and scheduling optimizations that attempt to sequentialize all I/O to disk. In particular, this benchmark, with multiple processes performing sequential I/O, is likely a good model for HEP/NP batch process scratch access. This is especially the case at RACF where some of the experiments we support stream input/output files to/from local scratch storage at the beginning and end of their batch jobs. For this study, we ran Bonnie++ as both a single process for sequential I/O testing, and as 24 synchronized processes in parallel for random I/O testing. Iozone can run in multithreaded (throughput mode), or single-threaded modes, and performs several I/O tests, including both sequential and random. Only single-threaded Iozone tests were run, and we were primarily interested in the random I/O benchmark results. Iozone was run in "automatic mode," where record size was varied multiple times during execution. Bonnie++ was run with only two record sizes: 8 KB, and single byte (character).  As the parallel Bonnie++ results were most relevant to our use case, the Bonnie++ benchmarks were run with both clean and dirty Flashcache and Bcache SSD caches. Dirty test results were obtained by running the benchmark a number of times in succession before the final run, where results were collected. Clean tests were first runs of the benchmark on a newly created device and filesystem. The Iozone results presented can be considered dirty, since the benchmark was run multiple times with varying filesizes before the final 64 GB file tests.

Results
While the single SSD provided excellent random I/O benchmark results, particularly for small record sizes, it did not provide the performance of the multi-spindle software RAID0 configuration for larger records. The software RAID0 configuration provided roughly double the random I/O performance, when compared to the SSD for large records, and for parallel workloads. However, it consisted of 8 times the number of drives. As expected, single SSD random I/O performance was significantly better than a single SATA drive.
Flashcache and Bcache with an SSD cache generally augmented the I/O performance of a single SATA disk for files that fit within the cache. Flashcache usually yielded better random write performance in the single disk configuration, while Bcache provided for better overall random read performance. Smaller gains, or performance losses were typically seen when the cache was preloaded with dirty data during Bonnie++ testing.
Fronting a 7-spindle software RAID0 array with a single SSD cache via Flashcache and Bcache generally reduced the performance of the array, when compared to an 8-spindle RAID0 device without an SSD cache. It's possible that Flashcache and Bcache may improve the performance of RAID5/6 devices. However, we did not test these RAID levels, as we were mainly interested in the performance of local scratch storage, where they are not commonly used. Bcache backed by a 7-spindle software RAID0 array came close to the performance of the 8-spindle cacheless software RAID0 device, and in some tests surpassed it. This is likely due to Bcache's monitoring of SSD congestion, and its default automatic cache bypass capability during periods of high latency SSD access. Figure 3. Bonnie++ parallel sequential read/write performance -generates a randomized workload.

Conclusions
We tested Bcache and Flashcache, two software hybrid device implementations for the Linux kernel, to determine if they could be used to address the continued increase in I/O demands from HEP/NP jobs running on worker nodes with progressively larger core/slot counts. We were primarily interested in increasing local worker node scratch storage performance, and/or reducing the cost of this storage. While it appears that Flashcache and Bcache can be used to improve the performance of a single SATA drive, these kernel drivers, with a single SSD cache, cannot provide the overall performance of multi-spindle software RAID0 devices, or be used to greatly improve the performance of these arrays. Therefore, while Flashcache and Bcache may be well suited for use in situations where sets of relatively small files are repeatedly and randomly accessed on a single drive or slow redundant storage, this software does not appear to be particularly well suited for use in HEP/NP worker node local scratch storage environments. Until the cost of high capacity SSDs significantly drops, or the I/O performance and/or density requirements for HEP/NP jobs is reduced (i.e. through the widespread adoption of multithreaded batch processing models), the use of software RAID0 arrays consisting of many traditional rotating drives appears to be the best option for local scratch storage.  Figure 4. Bonnie++ sequential read/write performance.