Yu. A. Kuznetsov

Fictitious domain and domain decomposition methods

The fictitious domains and domain decomposition methods are treated as iterative methods in subspaces, developed for solving the systems of linear algebraic equations. The first part of the paper is devoted to the algebraic theory of these two methods, and the second part to their application to solving large systems of grid equations. Detailed analysis is given of the application of the generalized conjugate gradient technique to increasing the rate of convergence of iterative procedures, and in particular, to its implementation as a computational process in a subspace. Estimates are given of the complexity of algorithm realizations in model problems. The imbedding of an original domain into a domain of a geometrically simpler form and its decomposition into subdomains are useful for the construction of methods for the solution of partial differential equations and of approximating systems of grid equations. These ideas are fairly popular in computational mathematics. The methods selected for this paper are those elaborated with active participation of the authors. They are the methods for the solution of systems of linear algebraic equations used as approximations of elliptic differential problems. The main feature of the proposed approach is the analysis and realization of iterative methods as computational processes in certain subspaces of the original vector space. This approach proved to be quite efficient and allowed us to design a number of effective (in the sense of the number of arithmetic operations) iterative methods [22, 23]. Such an approach was formerly used for the construction of efficient iterative methods for the solution of the neutron transport equation [26, 31]. The first section of this paper deals with the theory of the generalized conjugate gradient method in a subspace (GCGS) and also with some general aspects of its implementation in subspaces of special structure. This method, first formulated in [28, 30] and developed further in [20], generalizes the classical conjugate gradient method [13,14]. At the end of this section we propose a new approach to the implementation of the method; namely, by using three subspaces of, generally speaking, different structure, successively contained in each other. In Section 2 three versions of the fictitious component method are considered, using different extensions of the original system. Here, the presentation follows the original papers [24, 29, 32]. Note that similar versions of the method were proposed independently in [1, 2, 16] for studying specific grid equation systems. The methods described in Section 6 are based on the results of this section. Section 3 discusses three versions of the stationary block relaxation method which are then used as a basis for the construction of three corresponding versions of the generalized conjugate gradient method in a subspace. Two of these methods are well known: namely, the block Jacobi and the Gauss-Seidel methods [37]. The latter was 4 G. /. Marchuk et al. described, for example, in [17] as a method in subspace. The third version was constructed and analysed by a number of authors for different (grid, differential, operator) types of problems. The algebraic formulation of the method (in which we follow the papers [32, 33] is apparently proposed for the first time. Section 4 is concerned with the construction of systems of grid equations for a specific class of model problems. In the next section this system is used as a study case for the described iterative methods. In Section 5 we briefly discuss the concept of a partial solution of a system of linear algebraic equations. This notion was introduced in Section 1. Fast direct methods for the solution of specific partial problems are described in connection with the solution of problems involving matrices of a special type, appearing in the grid approximation of differential problems with separable variables. The formulation of these methods is based on the results of [3, 18, 25]. One of the most important computational features of these methods is the possibility of implementation with a significantly smaller amount of memory than the dimension of the original system. The application of the method of partial solutions is illustrated at the end of the section with a specific system of grid equations from Section 4. The last two sections deal with applications of the fictitious domain method (for the fictitious component method, see [27]) and of the domain decomposition method for the solution of the systems of grid equations from Section 4. Estimates of the convergence rate are constructed for these methods; the algorithms are formulated for their implementation in subspaces of small dimension, using partial solutions; and the estimates of computational cost (of computational complexity) and the amount of memory required for the implementation of these methods in a given subspace are obtained. It is worth noting that the described versions of the fictitious component method are closely connected with the capacitance matrix method [34]. 1. GENERALIZED CONJUGATE GRADIENT METHOD IN SUBSPACE Let R be the space of all N-dimensional real vectors with the usual scalar product N (U,V)= £ UfVi i= 1 and norm u = (u, w), w, veR. Let K be some subspace of R. An Ν χ N-matrix D is said to be self-adjoint in K if (Ώξ, η) = (ξ, Ώή) for any ξ,ηεΥ; it is said to be positive definite in Κ if (Ώξ, ξ)>0 for any non-zero ξεΥ. Likewise, an Ν χ JV-matrix D is said to be positive semidefinite in V if ( D & Q ^ O f o r a n y f e K . If V = R, matrix D is simply called self-adjoint (symmetric), positive definite, or positive semidefinite. If matrix D is self-adjoint and positive definite in a subspace K then we can define in V a new scalar product (ξ, η)0 = (Ώξ, η) and a norm Let Κ be invariant with respect to a matrix S, i.e. SV £ V. Then S is said to be D-selfadjoint in K if (

New algorithms for approximate realization of implicit difference schemes

ξ,η)Ώ = (ξ,8η)0, Υξ,τ/θΚ; it is said to be D-positive definite in V { (8ξ, ξ)Ό > 0 for any non-zero ξ e V. Obviously, S is D-self-adjoint (D-positive definite) in K if, and only if, the matrix DS is self-adjoint (positive definite) in the same subspace V. Fictitious domain and domain decomposition methods 5 Proposition 1 If matrix S is Z)-self-adjoint in V then V is the span of a system of Z)-orthonormal eigenvectors, corresponding to a set of eigenvalues of S. If matrix S is D-self-adjoint in its invariant subspace V then a set Λ(5, V) = {λΐ9...9λΜ} corresponds to S and K; namely, it is a set containing m = dim V real eigenvalues (possibly multiple) of 5 and such that the system of eigenvectors ψ !,..., ψη of 5, corresponding to these eigenvalues, forms a D-orthonormal basis in the subspace V. Proposition 2 Let us add to Proposition 1 an assumption of matrix S being D-positive definite in V. Then all eigenvalues of S in A(S, V) are positive. Hereafter we assume the eigenvalues of Λ(5, V) to be enumerated in non-decreasing order, i.e. λ1 ̂ λ2 *ζ ··· =ζ λη. Now let us consider a system of linear algebraic equations

On partial solution of systems of linear algebraic equations

The present paper treats the construction of estimates for Greens mesh function of finite difference and finite element operators involved in solving non-stationary heat conduction equations by implicit methods. Using these estimates, a new approach to approximate realization of implicit difference schemes is suggested on the basis of partitioning the spatial mesh domain into small mesh subdomains. Estimates are constructed for the computational cost of realizing implicit schemes with a prescribed accuracy. The present paper is devoted to approximate methods for solving systems of mesh equations involved in realization of implicit and explicit-implicit difference schemes for parabolic equations, namely, for problems of non-stationary heat conduction and diffusion. The particular feature of such mesh systems is that there is a small parameter, τ, a time step, at the mesh operator being a counterpart of the differential operator in spatial variables. If we denote by h a spatial mesh step, we can assume that τ takes a value in the interval Ic0h9clii], where c0 and cl are some positive constants. As a rule, we assume that τ is a function of ft, for example, τ = c h for implicit schemes and τ = c1h for explicit-implicit schemes. This fact makes it possible to obtain estimates for Greens mesh functions of the mesh systems considered, which indicate a sufficiently high rate of decrease of these functions depending on an increase in the distance from the mesh node at which the mesh source function is given. This paper shows that at the distance ~ σ^/τ Info from the source the value of Greens mesh function can be estimated by the quantity ~ h which, naturally, decreases as the distance from the source increases. The estimates obtained in the paper have enabled us to propose for approximate solution of the above-mentioned mesh systems new algorithms based on partitioning the mesh domain into mesh subdomains of a small size. In order to solve mesh systems in small subdomains, we suggested then different iterative methods which take into account the specific properties of the operators (matrices) of these systems. We have shown in the paper that asymptotically (for h «1 and, consequently, for τ «1) the cost of the approximate solution of the mesh systems to an accuracy of ε = h° amounts to Ο (σ h~ nlnh~) arithmetic operations. This result is, certainly, theoretical so far, the more so as it has been obtained using algebraic arguments. The impact of such approximate realization on the time stability of the difference scheme is not analysed, although it is obvious that the stability of the scheme thus perturbed exists at least for σ = 4. 1 Originally published in Russian as Preprint No. 142 of the Department of Numerical Mathematics of the USSR Academy of Sciences, Moscow, 1987. 100 Yu. A. Kuznetsov In Section 1 of this paper, we illustrate the approach suggested by considering the case of an one-dimensional heat conduction equation with constant coefficients. Estimates for Greens mesh function are constructed here on the basis of an explicit representation of the solution of a homogeneous mesh system. In Section 2, using the properties of monotone matrices (here, M-matrices) and the results of Section 1, we construct estimates for Greens mesh function and an approximate method of the solution of the mesh systems considered for the case of a two-dimensional heat conduction equation with constant coefficients subject to the Dirichlet boundary condition. In Section 3 which is the main part of the paper, we consider a sufficiently general two-dimensional heat conduction equation (variable coefficients, mixed boundary conditions). The well-known iterative technique, i.e. the conjugate gradient method, enables us to establish estimates for Greens mesh function which coincide with those in Section 2. Then making use of the idea of partitioning the mesh domain into small subdomains, we construct an algorithm for approximate solution of the original mesh system and estimate its arithmetic complexity. The paper ends with the discussion of some aspects of the generalization of the results obtained to domains with curvilinear boundaries and also their extension to the case of three-dimensional problems. 1. ONE-DIMENSIONAL HEAT CONDUCTION EQUATION Let us consider a very simple heat conduction equation where / =/(x, i) is a sufficiently smooth function with homogeneous boundary and initial conditions 11(0, i) = I<(1, t) = 0

Algebraic multigrid domain decomposition methods

Abstract—The problem of partial solution of systems of linear algebraic equations is formulated and consists in finding groups of components of the vector of the solution in the case of the right-hand side of the system belonging to a subspace. Such a problem arises, for example, in realization of block relaxation methods in subspaces. The problem is investigated in detail in the case of Jacobi matrices. The algorithms proposed are applied to constructing efficient direct methods for solving systems of difference equations approximating boundary value problems for Poissons equation in the rectangle.

Two-stage fictitious components method for solving the Dirichlet boundary value problem

The paper treats a new approach to constructing multilevel preconditioned for symmetric, positive definite and semi-definite matrices which arise in approximating elliptic problems by finite difference and finite element methods. The approach is based on multilevel domain decomposition methods with partitionings into small substructures and with inner Chebyshev iterative procedures. The paper shows that in the case of model elliptic boundary value problems, the iterative methods with the constructed preconditioners can be referred at the same time to the class of multigrid methods. Due to this fact these methods are called multigrid domain decomposition methods (MGDD-methods). The estimates of the convergence rate for the methods and of the computational cost of their realization are constructed for these model elliptic problems.

Block relaxation methods in subspaces, their optimization and application

The Dirichlet boundary value problem is considered for an elliptic equation with piecewisesmooth coefficients in twoand three-dimensional domains with complex internal and external curvilinear boundaries. To approximate the boundary value problem, the standard finite element method is used on rectangular meshes locally adapted to the boundaries. Systems of mesh equations arised are solved by a two-stage iterative method. This method involves the use of spectrally equivalent operators with constant coefficients as an outer iterative procedure and of the nonsymmetric version of the fictitious components method as an inner iterative procedure. The paper contains convergence rate estimates for the method discussed, proposes algorithms of its realization as a computational process in a subspace, gives estimates for the arithmetic and communication complexity of the algorithms suggested. The paper ends with the results of a numerical experiment to solve a specific three-dimensional problem of electrostatics. The fictitious components method has already a 15-year long history. Symmetric versions of this method were independently proposed in [2] and [17]. These papers also contained some convergence rate estimates for the symmetric fictitious components method. Especially note that the technique proposed in [2] and then developed in [3] for estimating the convergence rate of the symmetric fictitious components method for the Neumann boundary value problems makes up the foundation of present-day investigations in convergence of fictitious components methods and of many versions of domain decomposition methods. This technique is based on the theorems on continuation and traces of mesh functions, which are mesh analogues of the corresponding theorems in the theory of functions and functional analysis. The proof of such theorems has been given recently, for example, in [7,24,30]. The convergence rate estimates for the symmetric version of the fictitious components method for elliptic problems with the Dirichlet boundary conditions were obtained in [20,22]. The results of these investigations are used in this paper. The nonsymmetric version of the fictitious components method was independently proposed and investigated for solving elliptic problems with the Dirichlet boundary conditions in [11] and [23]. In [23], the author also established its close relation (duality, in a certain sense) to the symmetric fictitious components method for the Neumann boundary conditions. It was shown that the convergence rate estimates for the method constructed were implied by [2,3] and were also independent of the mesh size. Note that the symmetric and nonsymmetric versions of the fictitious components method are very close to capacitance matrix methods [4,25-27] by their idea, technique of constructing convergence rate estimates and realization algorithms. The two approaches have considerably enriched each other recently both in theory and in realization algorithms and applications. This paper is devoted to further development of theory and to applications of the fictitious components method (the capacitance matrix method) for elliptic problems 302 S. A. Finogenov and Yu. A. Kuznetsov with the Dirichlet boundary conditions for domains with curvilinear boundaries. It is a continuation of the previous investigations [13,14,19] made by the authors. The distinctive feature of these investigations lies in considering the fictitious components method for equations with piecewise-constant coefficients as a computational process in a subspace of a considerably smaller dimension as compared with the dimension of the original algebraic problem. The practical foundation for the realization of such computational processes is made up by the partial solution algorithms [5,13,15]. In Section 1 of this paper, the nonsymmetric fictitious components method is considered as an inner iterative procedure of another iterative method. The matrix properties which will be later used as preconditioners for an outer iterative method are investigated, and the estimates are given for the arithmetic and communication complexity of the algorithms. The results obtained in Section 1 are used for the construction of a two-stage iterative method for two-dimensional elliptic problems in Section 2 and for three-dimensional problems in Section 3. Algorithms for realizing two-stage processes are also discussed and the estimates of their complexity are given therein. Section 4 contains the results of a numerical experiment to solve a threedimensional problem with complex internal and external boundaries. 1. FICTITIOUS COMPONENTS METHOD Let Ω be a two-dimensional polygonal domain with the boundary ΟΩ. According to [8], let us consider the following Dirichlet boundary value problem: for a given geL2(il) find νεΗ^Ω) such that (1.1) Vi?°Vwdu= Jn Jo Embed Ω into a rectangle Π with the sides parallel to the coordinate axes, denote by ΟΠ the boundary of Π, and define a domain G = ΠΩ with the boundary dG. Construct in Π a rectangular, in general, non-uniform mesh ΠΛ assuming that for each mesh cell the ratio of the radius of the circle circumscribed to that of the circle inscribed is bounded from above by a mesh-independent constant. Then, partitioning the rectangular cells of the mesh ΠΛ into triangles, as shown, for example, in Fig. 1, Figure 1. A case of the mesh ΠΛ. Two-stage fictitious components method 303 construct a triangular mesh ΠΛ, and also triangular meshes ΩΛ = ΠήπΩ and Let us define 7Λ(Ω) as a space of functions continuous in Ω, linear in each triangle ΩΛ and vanishing on 5Ω. Consider the following approximate problem [8]: find ι/ΈΚΛ(Ω) such that Vv°Vwd&= gwdu VweVh(U). (1.2) Jn Jn This problem leads to a system of linear algebraic equations

Efficient fast direct method of solving Poisson’s equation on a parallelepiped and its implementation in an array processor

One more approach is proposed to numerical realization and optimization of block relaxation methods. This approach is based on the idea of subspaces naturally arising in realization of block relaxation methods. To optimize block relaxation methods, we make use of the generalized conjugate gradient method in subspace. The paper contains estimates of the convergence rate of the method and of the arithmetic cost of the solution of elliptic mesh systems by this method. The paper treats a new approach to algorithmic realization and optimization of block relaxation methods. The first section deals with general aspects of block relaxation methods and their optimization in subspaces (see [8, 9, 11]). The second section describes special algorithms of partial solution of special-form systems, which are closely related to the realization of the block relaxation method in subspace and gives the corresponding estimates of the number of operations and memory locations required. The third section shows the utilization of the method for solving systems of mesh equations arising in the approximation by the finite-difference method of the Dirichlet problem for the Laplace equation in a domain composed of rectangles. The forth section induces estimates of the convergence rate and the computational cost for a particular method of representation of the original domain in the form of a union of rectangles. It is shown here that in the case of an arbitrary non-uniform mesh in both variables, the estimate obtained coincides in order with the estimate of the successive overrelaxation method with the Chebyshev optimization (see [17,18]), and in a particular case providing for the utilization of the Fast Fourier Transform, it corresponds to the estimate of [1]. The final section outlines possibilities for generalization and development of the results of the paper. 1. GENERAL FORMULATION OF METHOD Let us consider a compatible system of linear algebraic equations Au = f (1.1) with the real symmetric and positive semi-definite matrix A of order n and the vector feimA = AEn. Here, by En we denote a space of η-dimensional real vectors with the scalar product (w, ι;) = Σ?= ιit d the norm || ν = (v, v), u, veEn. Write down the matrix A in the block form t Originally published in Russian in Variatsionno-Raznostnye Metody v Matematicheskoy Fizike. Comp. Cent. Sib. Branch, USSR Acad. Sei., Novosibirsk, 1978, pp. 178-212. 434 Yu. A. Kuznetsov

On some versions of the domain decomposition method

The paper treats problems of adjusting the structure of the computational algorithms to the architecture of a computer system with an array processor. The authors suggest a method of simulating the computational algorithms. Simulation makes it possible to determine the parameters of efficient implementation of the algorithm for the a priori prescribed level of performance of a specific computer system. As an example of the implementation of the method suggested, algorithms are analysed for the realization of a new fast direct method of solving Poissons equation on a parallelepiped. The results of numerical experiments are presented. Computer systems with attached processors (AP) constitute one of the important directions of progress in designing high-performance parallel computers. The attraction of such systems stems, first of all, from their relatively low cost and, second, from actually achievable performance comparable to that of supercomputers. These factors resulted in an impressive number of AP users. Nevertheless, the specifics of the architecture of such computer systems and the broad class of problems they solve generated difficulties in selecting methods of computation and in designing efficient implementation techniques. The decisive indices of efficiency in solving large-scale scientific and technical problems, covering the stages of formulation of the mathematical model through the processing of the results obtained, are the consumed astronomical time of software development and the run time for the problem. These cost indices greatly depend on the computational algorithms used and on the approach to designing the software. In most cases, an analysis of the properties of the algorithm is the necessary condition for its efficacious implementation on a computer system with AP. This paper suggests a special technique of simulating the computational algorithm on an AP computer system; this technique makes it possible to achieve an a priori prescribed real performance of the AP computer system and to improve the productivity of programmers developing an application-oriented software package. Such simulation answers several important questions: Does there exist an implementation of a given computational algorithm which achieves the prescribed real performance of a specific AP computer system? Assuming that this implementation exists, what software tool must be selected? If the given level of performance cannot be achieved, what are the reasons? As an example illustrating the technique suggested, we consider the solution of a system of linear algebraic equations arising in the approximation of three-dimensional Originally published in Russian as Preprint No. 147 of the Department of Numerical Mathematics of the USSR Academy of Sciences, Moscow, 1987. 2 A. A. Abakumov, A. Yu. Yeremin and Yu. A. Kuznetsov Poissons equation on a non-uniform rectangular grid. A new fast direct method is used for solving this system [6]. The efficacious realization of this method using AP computer systems is made possible by an advantageous combination of the parallel organization at the global and local levels. The paper is organized as follows. Section 1 describes the structure of a specific computer system with an attached processor and the software, and discusses the problems encountered in designing efficacious software. The simulation technique suggested by the authors is treated in Section 2. Section 3 presents a fast direct method of solving three-dimensional Poissons equation. Simulation of a computational algorithm using this method and a technique for designing and implementing the application-oriented software package for an AP computer system are described in Section 4. This section also gives the results of numerical experiments and analyses them. 1. DESCRIPTION OF THE COMPUTER SYSTEM AND PROGRAMMING TOOLS The architecture of the computer system consisting of a host scalar computer (HOST) and a high-performance attached processor ES2706 (of AP-190L type) [4,8] is schematically shown in Fig. 1. The following notation is used in Fig. 1 and throughout the text: ch—channel of data transfer between HOST and AP; P1-P4—four pages of data memory (DM), 64k words per page; TM—memory for constants (table memory), of 4k words capacity, with reduced access time; DPX, DPY—data pad registers, 32 words each; FAU—floating arithmetic unit with peak performance of 12 MFLOPS, consisting of addition and multiplication pipelines; PS—program source memory, with capacity of 2k 64-bits words. AP is composed of independent functional units controlled by a horizontal macroinstruction in lock-step mode. Up to ten subcommands can be initialized simultaneously for execution at each clock cycle. The following operations can be performed simultaneously: read or write the data on a data page, exchange the data located on another page with the host computer, read from TM, read/write on DPX and DPY, and compute using the pipeline FAU. The synchronization of functional units is organized by the user through programming tools. The extensive memory hierarchy is another important feature of the AP architecture. The data to be processed can be stored on four pages (P1-P4), each consisting of

Multigrid method for the plane problem of elasticity theory

The paper suggests a new approach to construction of preconditioned and convergence rate estimates for iterative methods to solve simultaneous linear algebraic equations arising in finite element approximations of elliptic problems. The approach is based on the decomposition of the original mesh domain into superelements. It is proved that for many important practical problems the convergence rate estimates of the methods considered do not depend on the mesh, and also on coefficients and boundary conditions of the original differential problem. INTRODUCTION Let Ω be a bounded two-dimensional polygonal domain with the boundary ΟΩ, and let Γ0 be a closed subset of <3Ω consisting of a finite number of straight line segments. Following [19] define the Sobolev space H = W(Q) and its subspace v = Qon Γ0}. In particular, if Γ0 = δΩ, then V = Hj =

Domain decomposition method to realize an implicit difference scheme for the one-phase Stefan problem

Τ(Ω and if Γ0 = 0, we have V = H. If mes Γ0 > 0, the relation f Γ ΓΥ&Λ /3ιΛΊ 1 ioH (o-i + hr~ r UnlA*!/ *2/ J J defines in the space V a norm which is equivalent to the H-norm and in the case where f ο = 0, the relation