Search

A widely used form of test matrix is the randsvd matrix constructed as the product $A = U \Sigma V^*$, where $U$ and $V$ are random orthogonal or unitary matrices from the Haar distribution and $\S...

A preconditioner for the mass matrix for high order finite element discretization on tetrahedra is presented and shown to give a condition number that is independent of both the mesh size and the p...

Methods are presented for automatically constructing coarse spaces of low dimension for domain decomposition algorithms. These constructions use equivalence classes of nodes on the interface between the subdomains into which the domain of a given elliptic problem has been subdivided, e.g., by a mesh partitioner such as METIS; these equivalence classes already play a central role in the design, analysis, and programming of many domain decomposition algorithms. The coarse space elements are well defined even for irregular subdomains, are continuous, and are well suited for use in two-level or multilevel preconditioners such as overlapping Schwarz algorithms. An analysis for scalar elliptic and linear elasticity problems reveals that significant reductions in the coarse space dimension can be achieved while not sacrificing the favorable condition number estimates for larger coarse spaces previously developed. These estimates depend primarily on the Lipschitz parameters of the subdomains. Numerical examples f...

This paper presents a new approach to constructing factorized approximate inverses for a symmetric and positive definite matrix $A$. The proposed strategy is based on adaptive dropping that reflects the quality of preserving the relation $UZ = I$ between the direct factor $U$ and the inverse factor $Z$ satisfying $A = U^TU$ and $A^{-1}=ZZ^T$. An important part of the approach is column pivoting, used to minimize the growth of the condition number of leading principal submatrices of $U$ that occurs explicitly in the dropping criterion. Numerical experiments demonstrate that the resulting approximate inverse factorization is robust as a preconditioner for solving large and sparse systems of linear equations.

An adaptive coarse space approach including a condition number bound for dual primal finite element tearing and interconnecting (FETI-DP) methods applied to three dimensional problems with coefficient jumps inside subdomains and across subdomain boundaries is presented. The approach is based on a known adaptive coarse space approach enriched by a small number of additional local edge eigenvalue problems. These edge eigenvalue problems serve to make the method robust and permit a condition number bound which depends only on the tolerance of the local eigenvalue problems and some properties of the domain decomposition. The introduction of the edge eigenvalue problems thus turns a well-known condition number indicator for FETI-DP and balancing domain decomposition by constraints (BDDC) methods into a condition number estimate. Numerical results are presented for linear elasticity and heterogeneous materials supporting our theoretical findings. The problems considered include those with random coefficients an...

To orthonormalize the columns of a dense matrix, the Cholesky QR (CholQR) requires only one global reduction between the parallel processing units and performs most of its computation using BLAS-3 kernels. As a result, compared to other orthogonalization algorithms, CholQR obtains superior performance on many of the current computer architectures, where the communication is becoming increasingly expensive compared to the arithmetic operations. This is especially true when the input matrix is tall-skinny. Unfortunately, the orthogonality error of CholQR depends quadratically on the condition number of the input matrix, and it is numerically unstable when the matrix is ill-conditioned. To enhance the stability of CholQR, we recently used mixed-precision arithmetic; the input and output matrices are in the working precision, but some of its intermediate results are accumulated in the doubled precision. In this paper, we analyze the numerical properties of this mixed-precision CholQR. Our analysis shows that ...

This paper introduces a high-order time stepping technique for solving the incompressible Navier--Stokes equations which, unlike coupled techniques, does not require solving a saddle point problem at each time step and, unlike projection methods, does not produce splitting errors and spurious boundary layers. The technique is a generalization of the artificial compressibility method; it is unconditionally stable (for the unsteady Stokes equations), can reach any order in time, and uncouples the velocity and the pressure. The condition number of the linear systems associated with the fully discrete vector-valued problems to be solved at each time step scales like $O(\tau h^{-2})$, where $\tau$ is the time step and $h$ is the spatial grid size. No Poisson problem or other second-order elliptic problem has to be solved for the pressure corrections. Unlike projection methods, optimal convergence is observed numerically with Dirichlet and mixed Dirichlet/Neumann boundary conditions.

The 2-Lagrange multiplier method is a domain decomposition method which can be used to parallelize the solution of linear problems arising from partial differential equations. In order to scale to large numbers of subdomains and processors, domain decomposition methods require a coarse grid correction to transport low frequency information more rapidly between subdomains that are far apart. We introduce two new 2-level methods by adding a coarse grid correction to 2-Lagrange multiplier methods. We prove that if we shrink $h$ (the grid parameter) while maintaining bounded the ratio $\frac{H}{h}$ (where $H$ is the size of the subdomains), the condition number of the method remains bounded. We confirm our analysis with experiments on the HECToR (High-End Computing Terascale Resource) supercomputer. This proves that the new methods scale weakly, opening the door to massively parallel implementations.

We study the solution of the linear least-squares problem $\min_x \| b - A x \|^2_2$ where the matrix $A \in {\Bbb R}^{m \times n}$ ($ m \ge n$) has rank $n$ and is large and sparse. We assume that $A$ is available as a matrix, not an operator. The preconditioning of this problem is difficult because the matrix $A$ does not have the properties of differential problems that make standard preconditioners effective. Incomplete Cholesky techniques applied to the normal equations do not produce a well-conditioned problem. We attempt to bypass the ill-conditioning by finding an $n \times n$ nonsingular submatrix $B$ of $A$ that reduces the Euclidean norm of $AB^{-1}$. We use $B$ to precondition a symmetric quasi-definite linear system whose condition number is then independent of the condition number of $A$ and has the same solution as the original least-squares problem. We illustrate the performance of our approach on some standard test problems and show it is competitive with other approaches.

Elliptic variational inequalities with multiple bodies in two dimensions are considered. It is assumed that an active set method is used to handle the nonlinearity of the inequality constraint, which results in auxiliary linear problems. For solving such linear problems we study two domain decomposition methods called the finite element tearing and interconnecting (FETI-FETI) and hybrid methods in this paper. Bodies are decomposed into several subdomains in both methods. The FETI-FETI method combines the one-level FETI and the dual-primal FETI (FETI-DP) methods. We present a proof that this combined method has a condition number that depends linearly on the number of subdomains across each body and polylogarithmically on the number of elements across each subdomain. Our numerical results, and those of others, suggest that this is the best possible bound. The hybrid method combines the one-level FETI and the balanced domain decomposition by constraints (BDDC) methods; we prove that the condition number of ...

We analyze the performance of the minimal residual (MINRES) method applied to linear Karush--Kuhn--Tucker systems arising in connection with inverse problems. Such optimality systems typically have a saddle point structure and have unique solutions for all $\alpha>0$, where $\alpha$ is the parameter employed in the Tikhonov regularization. Unfortunately, the associated spectral condition number is very large for small values of $\alpha$, which strongly indicates that their numerical treatment is difficult. Our main result shows that a broad range of linear ill posed optimality systems can be solved efficiently with the MINRES method. This result is obtained by carefully analyzing the spectrum of the associated saddle point operator: Except for a few isolated eigenvalues, the spectrum consists of three bounded intervals. Krylov subspace methods handle such problems very well. For severely ill posed cases, techniques based on Chebyshev polynomials are applied to prove that the number of iterations needed by...

We present a robust method for choosing multivariate polynomial interpolation nodes. Our algorithm is an optimization method to greedily minimize a measure of interpolant sensitivity, a variant of a weighted Lebesgue function. Nodes are therefore chosen that tend to control oscillations in the resulting interpolant. This method can produce an arbitrary number of nodes and is not constrained by the dimension of a complete polynomial space. Our method is therefore flexible: nested nodal sets are produced in spaces of arbitrary dimensions, and the number of nodes added at each stage can be arbitrary. The algorithm produces a nodal set given a probability measure on the input space, thus parameterizing interpolants with respect to finite measures. We present examples to show that the method yields nodal sets that behave well with respect to standard interpolation diagnostics: the Lebesgue constant, the Vandermonde determinant, and the Vandermonde condition number. We also show that a nongreedy version of the ...

The increasing computational load required by most applications and the limits in hardware performances affecting scientific computing contributed in the last decades to the development of parallel software and architectures. In fluid-structure interaction (FSI) for haemodynamic applications, parallelization and scalability are key issues (see [L. Formaggia, A. Quarteroni, and A. Veneziani, eds., Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System, Modeling, Simulation and Applications 1, Springer, Milan, 2009]). In this work we introduce a class of parallel preconditioners for the FSI problem obtained by exploiting the block-structure of the linear system. We stress the possibility of extending the approach to a general linear system with a block-structure, then we provide a bound in the condition number of the preconditioned system in terms of the conditioning of the preconditioned diagonal blocks, and finally we show that the construction and evaluation of the devised preconditioner is modular. The preconditioners are tested on a benchmark three-dimensional (3D) geometry discretized in both a coarse and a fine mesh, as well as on two physiological aorta geometries. The simulations that we have performed show an advantage in using the block preconditioners introduced and confirm our theoretical results.

A scheme for the numerical solution of singular integral equations on piecewise smooth curves is presented. It relies on several techniques: reduction, Nyström discretization, composite quadrature, recursive compressed inverse preconditioning, and multipole acceleration. The scheme is fast and stable. Its computational cost grows roughly logarithmically with the precision sought and linearly with overall system size. When the integral equation models a boundary value problem, the achievable accuracy may be close to the condition number of that problem times machine epsilon. This is illustrated by application to elastostatic problems involving zigzag-shaped cracks with up to twenty thousand corners and branched cracks with hundreds of triple junctions.

We provide formulas for the norm of the error when solving nonsymmetric linear systems with the full orthogonalization method (FOM) and the generalized minimum residual method (GMRES) as well as relations between the error norm and the residual norm. From these formulas we are able to compute estimates of the norm of the error during the iterations. Since stopping criteria based on the norm of the residual may sometimes be misleading, such estimates could lead to a more robust way to stop the iterations. Numerical experiments show that the proposed norm estimates work nicely on difficult linear systems.

CG, SYMMLQ, and MINRES are Krylov subspace methods for solving symmetric systems of linear equations. When these methods are applied to an incompatible system (that is, a singular symmetric least-squares problem), CG could break down and SYMMLQ's solution could explode, while MINRES would give a least-squares solution but not necessarily the minimum-length (pseudoinverse) solution. This understanding motivates us to design a MINRES-like algorithm to compute minimum-length solutions to singular symmetric systems. MINRES uses QR factors of the tridiagonal matrix from the Lanczos process (where R is upper-tridiagonal). MINRES-QLP uses a QLP decomposition (where rotations on the right reduce R to lower-tridiagonal form). On ill-conditioned systems (singular or not), MINRES-QLP can give more accurate solutions than MINRES. We derive preconditioned MINRES-QLP, new stopping rules, and better estimates of the solution and residual norms, the matrix norm, and the condition number., Comment: 26 pages, 6 figures

The simulation of drift instabilities in the plasma edge leads to cubic polynomial PDE eigenvalue problems with parameter dependent coefficients. The aim is to determine the wave number which leads to the maximum growth rate of the amplitude of the wave. This requires the solution of a large number of PDE eigenvalue problems. Since we are only interested in a smooth eigenfunction corresponding to the eigenvalue with largest imaginary part, the Jacobi-Davidson method can be applied. Unfortunately, a naive implementation of this method is much too expensive for the large number of problems that have to be solved. In this paper we will present a multilevel approach for the construction of an appropriate initial search space. We will also discuss the efficient solution of the correction equation, and we will show how optimal scaling helps to accelerate the convergence.

A FETI-DP (dual-primal finite element tearing and interconnecting) algorithm for the three-dimensional Stokes problem is developed and analyzed. This is an extension of the previous work for the two-dimensional problem in [H. H. Kim, C.-O. Lee, and E.-H. Park, SIAM J. Numer. Anal., 47 (2010), pp. 4142-4162]. Advantages of this approach are the coarse problem without primal pressure unknowns and the use of a computationally cheap lumped preconditioner. Especially in three dimensions, these advantages provide a more practical FETI-DP algorithm. In three dimensions, the velocity unknowns at subdomain corners and the averages of velocity unknowns over common faces are selected as the primal unknowns in the FETI-DP formulation. Its condition number bound is analyzed to be $CH/h$, where $C$ is a positive constant which is independent of any mesh parameters and $H/h$ is the number of elements across each subdomain. Numerical results are included.

We develop and compare multilevel algorithms for solving constrained nonlinear variational problems via interior point methods. Several equivalent formulations of the linear systems arising at each iteration of the interior point method are compared from the point of view of conditioning and iterative solution. Furthermore, we show how a multilevel continuation strategy can be used to obtain good initial guesses ("hot starts") for each nonlinear iteration. Some minimal surface and volume-constrained image registration problems are used to illustrate the various approaches. © 2009 Society for Industrial and Applied Mathematics.

Balancing domain decomposition by constraints (BDDC) algorithms are constructed and analyzed for the system of almost incompressible elasticity discretized with Gauss-Lobatto-Legendre spectral elements in three dimensions. Initially mixed spectral elements are employed to discretize the almost incompressible elasticity system, but a positive definite reformulation is obtained by eliminating all pressure degrees of freedom interior to each subdomain into which the spectral elements have been grouped. Appropriate sets of primal constraints can be associated with the subdomain vertices, edges, and faces so that the resulting BDDC methods have a fast convergence rate independent of the almost incompressibility of the material. In particular, the condition number of the BDDC preconditioned operator is shown to depend only weakly on the polynomial degree $n$, the ratio $H/h$ of subdomain and element diameters, and the inverse of the inf-sup constants of the subdomains and the underlying mixed formulation, while being scalable, i.e., independent of the number of subdomains and robust, i.e., independent of the Poisson ratio and Young's modulus of the material considered. These results also apply to the related dual-primal finite element tearing and interconnect (FETI-DP) algorithms defined by the same set of primal constraints. Numerical experiments, carried out on parallel computing systems, confirm these results.

We present two new algorithms FastAccSum and FastPrecSum, one to compute a faithful rounding of the sum of floating-point numbers and the other for a result “as if” computed in $K$-fold precision. Faithful rounding means the computed result either is one of the immediate floating-point neighbors of the exact result or is equal to the exact sum if this is a floating-point number. The algorithms are based on our previous algorithms AccSum and PrecSum and improve them by up to 25%. The first algorithm adapts to the condition number of the sum; i.e., the computing time is proportional to the difficulty of the problem. The second algorithm does not need extra memory, and the computing time depends only on the number of summands and $K$. Both algorithms are the fastest known in terms of flops. They allow good instruction-level parallelism so that they are also fast in terms of measured computing time. The algorithms require only standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision.

Running relative error for evaluating polynomials is analyzed and applied to different evaluation algorithms: Horner, VS, de Casteljau, and Clenshaw algorithms. Numerical experiments show that the derived running relative error bounds are sharp and that the de Casteljau algorithm presents better stability properties than the remaining algorithms.

In this paper, we describe and analyze a novel tensor approximation method for discretized multidimensional functions and operators in $\mathbb{R}^d$, based on the idea of multigrid acceleration. The approach stands on successive reiterations of the orthogonal Tucker tensor approximation on a sequence of nested refined grids. On the one hand, it provides a good initial guess for the nonlinear iterations to find the approximating subspaces on finer grids; on the other hand, it allows us to transfer from the coarse-to-fine grids the important data structure information on the location of the so-called most important fibers in directional unfolding matrices. The method indicates linear complexity with respect to the size of data representing the input tensor. In particular, if the target tensor is given by using the rank-$R$ canonical model, then our approximation method is proved to have linear scaling in the univariate grid size $n$ and in the input rank $R$. The method is tested by three-dimensional (3D) electronic structure calculations. For the multigrid accelerated low Tucker-rank approximation of the all electron densities having strong nuclear cusps, we obtain high resolution of their 3D convolution product with the Newton potential. The accuracy of order $10^{-6}$ in max-norm is achieved on large $n\times n\times n$ grids up to $n=1.6\cdot10^4$, with the time scale in several minutes.

We present two algorithms for computing correctly rounded sums of arrays of floating-point numbers. First, iFastSum improves upon our previous FastSum by requiring no additional space beyond the original array, which is destroyed. It runs about 20% faster than FastSum in the general case and two times faster when extremely ill-conditioned data are used. The second algorithm is HybridSum, which combines three summation ideas together: splitting the mantissa, radix sorting, and using iFastSum. The result is that when the number of summands is greater than about $10^4$, for a given $n$ its running time is almost a constant, independent of the condition number. It runs almost as fast as iFastSum in the general case and much faster than iFastSum when ill-conditioned data are used. HybridSum requires only one pass through the input array and uses constant storage, and it is thus suitable for exact summation as an “online” algorithm. Neither algorithm requires extra precision accumulators, and both work in any base. Their accuracy is guaranteed independent of the condition number and the number of summands.

We show that several methods for solving time of arrival (TOA) geolocation problems depend crucially on what we call the large equal radius (LER) conditions. For example, we analyze Bancroft's method and show that it works well assuming the LER conditions are satisfied. However, when these conditions are violated, we show that the method significantly deteriorates and does not achieve the minimum possible variance for a least squares solution. A direct method is also presented that finds all the critical points of the least squares functional. The method is found to perform well when the LER conditions are not satisfied but suffers from conditioning problems when they are satisfied. We also show that a simple but statistically inferior method that we call the reduced least squares (RLS) method deteriorates when the LER conditions are satisfied.

When radial basis functions (RBFs) are made increasingly flat, the interpolation error typically decreases steadily until some point when Runge-type oscillations either halt or reverse this trend. Because the most obvious method to calculate an RBF interpolant becomes a numerically unstable algorithm for a stable problem in the case of near-flat basis functions, there will typically also be a separate point at which disastrous ill-conditioning enters. We introduce here a new method, RBF-QR, which entirely eliminates such ill-conditioning, and we apply it in the special case when the data points are distributed over the surface of a sphere. This algorithm works even for thousands of node points, and it allows the RBF shape parameter to be optimized without the limitations imposed by stability concerns. Since interpolation in the flat RBF limit on a sphere is found to coincide with spherical harmonics interpolation, new insights are gained as to why the RBF approach (with nonflat basis functions) often is the more accurate of the two methods.

A postprocessing technique based on negative order norm estimates for the discontinuous Galerkin methods was previously introduced by Cockburn, Luskin, Shu, and Suli [Proceedings of the International Symposium on Discontinuous Galerkin Methods, Springer, New York, pp. 291-300; Math. Comput., 72 (2003), pp. 577-606]. The postprocessor allows improvement in accuracy of the discontinuous Galerkin method for time-dependent linear hyperbolic equations from order $k$+1 to order 2$k$+1 over a uniform mesh. Assumptions on the convolution kernel along with uniformity in mesh size give a local translation invariant postprocessor that allows for simple implementation using small matrix-vector multiplications. In this paper, we present two alternatives for extending this postprocessing technique to include smoothly varying meshes. The first method uses a simple local $L^2$-projection of the smoothly varying mesh to a locally uniform mesh and uses this projected solution to compute the postprocessed solution. By using this local $L^2$-projection, recalculating the convolution kernel for every element can be avoided, and 2$k$+1 order accuracy of the postprocessed solution can be achieved. The second method uses the idea of characteristic length based upon the largest element size for the scaling of the postprocessing kernel. These two methods, local projection and characteristic length, are also applied to approximations over a mesh with elements that vary in size randomly. We discuss the computational issues in using these two techniques and demonstrate numerically that we obtain the 2$k$+1 order of accuracy for the smoothly varying meshes, and that although the 2$k$+1 order of accuracy is not fully realized for random meshes, there is significant improvement in the $L^2$-errors.

In this paper we present a new incomplete factorization of a square matrix into triangular factors in which we get standard $LU$ or $LDL^T$ factors (direct factors) and their inverses (inverse factors) at the same time. Algorithmically, we derive this method from the approach based on the Sherman-Morrison formula [R. Bru, J. Cerdan, J. Marin, and J. Mas, SIAM J. Sci. Comput., 25 (2003), pp. 701-715]. In contrast to the robust incomplete decomposition (RIF) algorithm [M. Benzi and M. Tůma, Numer. Linear Algebra Appl., 10 (2003), pp. 385-400] the direct and inverse factors here directly influence each other throughout the computation. Consequently, the algorithm to compute the approximate factors may mutually balance dropping in the factors and control their conditioning in this way. For the symmetric positive definite case, we derive the theory and present an algorithm for computing the incomplete $LDL^T$ factorization, and we discuss experimental results. We call this new approximate $LDL^T$ factorization the balanced incomplete factorization (BIF). Our experimental results confirm that this factorization is very robust and may be useful in solving difficult ill conditioned problems by preconditioned iterative methods. Moreover, the internal coupling of the computation of direct and inverse factors results in much shorter setup times (times to compute approximate decomposition) than RIF, a method of a similar and very high level of robustness. We also derive and present the theory for the general nonsymmetric case, but do not discuss its implementation.

In this paper we construct optimal preconditioners for the discrete mathematical models arising in modeling the elastic responses of lattice block materials. We present extensive numerical experiments to show that the preconditioned system has a uniformly bounded condition number with respect to the size of problem and with respect to the parameter relating the stretching and bending of the beams in a lattice. Using the limiting system of partial differential equations, we show theoretically that for square lattices the proposed preconditioners are efficient by proving a uniform bound on the condition number of the preconditioned system.

Given a vector of floating-point numbers with exact sum $s$, we present an algorithm for calculating a faithful rounding of $s$, i.e., the result is one of the immediate floating-point neighbors of $s$. If the sum $s$ is a floating-point number, we prove that this is the result of our algorithm. The algorithm adapts to the condition number of the sum, i.e., it is fast for mildly conditioned sums with slowly increasing computing time proportional to the logarithm of the condition number. All statements are also true in the presence of underflow. The algorithm does not depend on the exponent range. Our algorithm is fast in terms of measured computing time because it allows good instruction-level parallelism, it neither requires special operations such as access to mantissa or exponent, it contains no branch in the inner loop, nor does it require some extra precision: The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithm are proved to be optimal.

This paper investigates the performance of tensor methods for solving small- to large-scale systems of nonlinear equations where the Jacobian matrix at the root is ill-conditioned or singular. This condition occurs on many classes of problems, such as identifying or approaching turning points in path-following problems. The singular case has been studied more than the highly ill-conditioned case, for both Newton and tensor methods. It is known that Newton-based methods do not work well with singular problems because they converge linearly to the solution and, in some cases, with poor accuracy. On the other hand, direct tensor methods have performed well on singular problems and have superlinear convergence on such problems under certain conditions. This behavior originates from the use of a special, restricted form of the second-order term included in the local tensor model that provides information lacking in a (nearly) singular Jacobian. With several implementations available for large-scale problems, tensor methods now are capable of solving larger problems. We compare the performance of tensor methods and Newton-based methods for small- to large-scale problems over a range of conditionings, from well-conditioned to ill-conditioned to singular. Previous studies with tensor methods concerned only the ends of this spectrum. Our results show that tensor methods are increasingly superior to Newton-based methods as the problem grows more ill-conditioned.

The established expressions for the condition numbers of a root of a polynomial assume that the perturbations in its coefficients are small, such that only the lowest order terms in a Taylor series need be considered. This restriction to the lowest order terms is considered in this paper, and it is shown that there exists a well-defined class of polynomial for which it is not valid. In particular, if a root $x_{0}$ of a polynomial $p(x)$ is close to but distinct from its neighboring roots $x_{i},i=1,\dots,r$, it is shown that it is necessary to consider higher order terms in the Taylor series. These additional terms imply that a closed form expression for the condition number of $x_{0}$ cannot be obtained, and numerical methods must therefore be employed. In this paper, a high order Taylor series is used to compute the region in the complex plane in which the perturbed root $x_{0} + \delta x_{0}$ lies, and this enables the variation of the condition number of $x_{0}$ with the distances to its neighboring roots to be established. It is shown that this extra information, which cannot be obtained if only the lowest order terms are used, is significant.

Maximum weight matchings have become an important tool for solving highly indefinite unsymmetric linear systems, especially in direct solvers. In this study we investigate the benefit of reorderings and scalings based on symmetrized maximum weight matchings as a preprocessing step for incomplete $\mathrm{LDL^T}$ factorizations. The reorderings are constructed such that the matched entries form $1 \times 1$ or $2 \times 2$ diagonal blocks in order to increase the diagonal dominance of the system. During the incomplete factorization only tridiagonal pivoting is used. We report results for this approach and comparisons with other solution methods for a diverse set of symmetric indefinite matrices, ranging from nonlinear elasticity to interior point optimization.

A FETI-DP (dual-primal finite element tearing and interconnecting) formulation for the two-dimensional Stokes problem with mortar methods is considered. Separate sets of unknowns are used for velocity on interfaces, and the mortar constraints are enforced on the velocity unknowns by Lagrange multipliers. Average constraints on edges are further introduced as primal constraints to solve the Stokes problem correctly and to obtain a scalable FETI-DP algorithm. A Neumann--Dirichlet preconditioner is shown to give a condition number bound, C maxi = 1,...,N {(1 + log( Hi / hi))2}, where Hi and hi are the subdomain size and the mesh size, respectively, and the constant C is independent of the mesh parameters Hi and hi.

In this paper we compare various preconditioners for the numerical solution of partial differential equations. We compare the well-known balancing preconditioner used in domain decomposition methods with a so-called deflation preconditioner. We prove that the effective condition number of the deflated preconditioned system is always, i.e., for all deflation vectors and all restrictions and prolongations, below the condition number of the system preconditioned by the balancing preconditioner. Even more, we establish that both preconditioners lead to almost the same spectra. The zero eigenvalues of the deflation preconditioned system are replaced by eigenvalues which are one if the balancing preconditioner is used. Moreover, we prove that the A-norm of the errors of the iterates built by the deflation preconditioner is always below the A-norm of the errors of the iterates built by the balancing preconditioner. Depending on the implementation of the balancing preconditioner the amount of work of one iteration of the deflation preconditioned system is less than or equal to the amount of work of one iteration of the balancing preconditioned system. If the amount of work is equal, both preconditioners are sensitive with respect to inexact computations. Finally, we establish that the deflation preconditioner and the balancing preconditioner produce the same iterates if one uses certain starting vectors. Numerical results for porous media flows emphasize the theoretical results.

In this paper we propose and analyze a wavelet-based numerical method for exact controllability of the one-dimensional wave equation, concentrating on the particular case of Dirichlet controls. We use the Hilbert uniqueness method (HUM) and prove this result in the context of the finite element space semidiscretization with spline wavelet basis functions. The algorithm for the numerical construction of the controls is based on a conjugate gradient method proposed by Glowinski [J. Comput. Phys., 103 (1992), pp. 189-221] and a two-grid filtering technique in the wavelet basis. We show that wavelet filtering yields numerically exact controllability in the optimal time of the continuous problem and the convergence of the numerical controls using a multilevel discretization filtering technique. Classical arguments then allow proving the uniform boundedness of the numerical controls and passing to the limit as the mesh size tends to zero. We prove that the condition number of the discretized HUM operator $\Lambda$ with wavelet filtering is uniformly bounded with respect to the mesh width and prove convergence of the numerical controls obtained with the two-grid filtering technique. Numerical experiments show that wavelet filtering yields numerical controls with the same control time $T$ as in the continuous case.

Engineering tests generally acquire a large set of response data at once while focusing on a small number of variables. Higher-order response models are better for capturing and displaying the characteristics such data can reveal, but response surface models traditionally are linear or quadratic (first or second order) only. If attempted for a higher order, the matrix X of the test inputs formed for fitting the test responses often immediately stops satisfying the full-rank assumption on X traditionally made for computing the regression coefficients. The resulting product matrix XT X is singular and hence cannot be inverted, directly or indirectly, for obtaining the desired least-squares regression coefficients from the system of normal equations. Ill-conditioning of X when near rank deficiency, not unusual even in the traditional second-order response modeling, can also render the computed coefficients unreliable. This paper incorporates singular value decomposition (SVD) to handle the rank deficiency and ill-conditioning problems directly. The proposed SVD-based fitting algorithm enables response modeling to the desirable high-order reliably and efficiently with no difficulty despite rank deficiency in X. Two proposed trimming algorithms enable reduction of redundant terms or insignificant coefficients of the fitted polynomials on a sound mathematical base.