Direct and approximate algorithmic methods for solving self-adjoint elliptic PDEs in three-space dimensions

The numerical implementation of the extended to the limit sparse LDL^T factorization solution methods for three-dimensional self-adjoint elliptic partial differential equations [3] is given. Two FORTRAN routines for the approximate (or exact) factorization of the coefficient matrix and solution of the resulting finite difference equations are supplied. The amount of fill-in terms can be controlled by the user through parameters R1, R2 the limiting case being when the matrix is factorized exactly.     

Optimal damping of selected eigenfrequencies using dimension reduction

We consider linear vibrational systems described by a system of second‐order differential equations of the form , where M and K are positive definite matrices, representing mass and stiffness, respectively. The damping matrix D is assumed to be positive semidefinite. We are interested in finding an optimal damping matrix that will damp a certain (critical) part of the eigenfrequencies. For this, we use an optimization criterion based on the minimization of the average total energy of the system. This is equivalent to the minimization of the trace of the solution of the corresponding Lyapunov equation AX + XA^T = ?GG^T, where A is the matrix obtained from linearizing the second‐order differential equation, and G depends on the critical part of the eigenfrequencies to be damped. The main result is the efficient approximation and the corresponding error bound for the trace of the solution of the Lyapunov equation obtained through dimension reduction, which includes the influence of the right‐hand side GG^T and allows us to control the accuracy of the trace approximation. This trace approximation yields a very accelerated optimization algorithm for determining the optimal damping. Copyright © 2011 John Wiley & Sons, Ltd.     

Solving Differential Matrix Equations using Parareal

Differential matrix equations appear in many applications like optimal control of partial differential equations, balanced truncation model order reduction of linear time varying systems and many more. Here, we will focus on differential Riccati equations (DRE). Solving such matrix-valued ordinary differential equations (ODE) is a highly time consuming process. We present a Parareal based algorithm applied to Rosenbrock methods for the solution of the matrix-valued differential Riccati equations. Considering problems of moderate size, direct matrix equation solvers for the solution of the algebraic Lyapunov equations arising inside the time intgration methods are used. (© 2016 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

On large-scale diagonalization techniques for the Anderson model of localization

We propose efficient preconditioning algorithms for an eigenvalue problem arising in quantum physics, namely the computation of a few interior eigenvalues and their associated eigenvectors for large-scale sparse real and symmetric indefinite matrices of the Anderson model of localization. Our preconditioning approaches for the shift-and-invert symmetric indefinite linear system are based on maximum weighted matchings and algebraic multi-level incomplete LDL^T factorizations. These techniques can be seen as a complement to the alternative idea of using more complete pivoting techniques for the highly ill-conditioned symmetric indefinite Anderson matrices. Our numerical examples reveal that recent algebraic multi-level preconditioning solvers can accelerate the computation of a large-scale eigenvalue problem corresponding to the Anderson model of localization by several orders of magnitude. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

ADI preconditioned Krylov methods for large Lyapunov matrix equations

In the present paper, we propose preconditioned Krylov methods for solving large Lyapunov matrix equations AX+XA^T+BB^T=0. Such problems appear in control theory, model reduction, circuit simulation and others. Using the Alternating Direction Implicit (ADI) iteration method, we transform the original Lyapunov equation to an equivalent symmetric Stein equation depending on some ADI parameters. We then define the Smith and the low rank ADI preconditioners. To solve the obtained Stein matrix equation, we apply the global Arnoldi method and get low rank approximate solutions. We give some theoretical results and report numerical tests to show the effectiveness of the proposed approaches.     

Boundary characteristic point regularity for semilinear reaction-diffusion equations: towards an ode criterion

The classical problem of regularity of boundary characteristic points for semilinear heat equations with homogeneous Dirichlet conditions is considered. The Petrovskii ( 2?{loglog} ) \left( {2\sqrt {{\log \log }} } \right) criterion (1934) of the boundary regularity for the heat equation can be adapted to classes of semilinear parabolic equations of reaction–diffusion type and takes the form of an ordinary differential equation (ODE) regularity criterion. Namely, after a special matching with a boundary layer, the regularity problem reduces to a onedimensional perturbed nonlinear dynamical system for the first Fourier-like coefficient of the solution in an inner region. A similar ODE criterion, with an analogous matching procedures, is shown formally to exist for semilinear fourth order biharmonic equations of reaction-diffusion type. Extensions to regularity problems of backward paraboloid vertices in \mathbbR^N {\mathbb{R}^N} are discussed. Bibliography: 54 titles. Illustrations: 1 figure.     

Detection of near-singularity in Cholesky and LDLT factorizations

In sparse matrix applications it is often important to implement the Cholesky and LDL^T factorization methods without pivoting in order to avoid excess fillin. We consider methods for detection of a nearly singular matrix by means of these factorizations without pivoting and demonstrate that a technique based on estimation of the smallest eigenvalue via inverse iteration will always reveal a nearly singular matrix.     

Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization

Interior-point methods are among the most efficient approaches for solving large-scale nonlinear programming problems. At the core of these methods, highly ill-conditioned symmetric saddle-point problems have to be solved. We present combinatorial methods to preprocess these matrices in order to establish more favorable numerical properties for the subsequent factorization. Our approach is based on symmetric weighted matchings and is used in a sparse direct LDL ^T factorization method where the pivoting is restricted to static supernode data structures. In addition, we will dynamically expand the supernode data structure in cases where additional fill-in helps to select better numerical pivot elements. This technique can be seen as an alternative to the more traditional threshold pivoting techniques. We demonstrate the competitiveness of this approach within an interior-point method on a large set of test problems from the CUTE and COPS sets, as well as large optimal control problems based on partial differential equations. The largest nonlinear optimization problem solved has more than 12 million variables and 6 million constraints.     

On splitting of linear nonstationary systems of differential equations

Using the matrix asymptotic method, we solve the problem of the splitting of n-order linear nonstationary systems of differential equations into n first-order independent equations for the case where the eigenvalues m_i( t) 1 0,i = 1,...,n,t ? [ 0,l ] {\mu_i}\left( \tau \right) \ne 0,i = 1,...,n,\tau \in \left[ {0,\ell } \right] , of the matrix U( t) U\left( \tau \right) determined by the initial system of differential equations are identically equal.     

Limited memory preconditioners for symmetric indefinite problems with application to structural mechanics

This paper presents a class of limited memory preconditioners (LMP) for solving linear systems of equations with symmetric indefinite matrices and multiple right‐hand sides. These preconditioners based on limited memory quasi‐Newton formulas require a small number k of linearly independent vectors and may be used to improve an existing first‐level preconditioner. The contributions of the paper are threefold. First, we derive a formula to characterize the spectrum of the preconditioned operator. A spectral analysis of the preconditioned matrix shows that the eigenvalues are all real and that the LMP class is able to cluster at least k eigenvalues at 1. Secondly, we show that the eigenvalues of the preconditioned matrix enjoy interlacing properties with respect to the eigenvalues of the original matrix provided that the k linearly independent vectors have been prior projected onto the invariant subspaces associated with the eigenvalues of the original matrix in the open right and left half‐plane, respectively. Third, we focus on theoretical properties of the Ritz‐LMP variant, where Ritz information is used to determine the k vectors. Finally, we illustrate the numerical behaviour of the Ritz limited memory preconditioners on realistic applications in structural mechanics that require the solution of sequences of large‐scale symmetric saddle‐point systems. Numerical experiments show the relevance of the proposed preconditioner leading to a significant decrease in terms of computational operations when solving such sequences of linear systems. A saving of up to 43% in terms of computational effort is obtained on one of these applications. Copyright © 2016 John Wiley & Sons, Ltd.     

Perturbation of Eigenpairs of Factored Symmetric Tridiagonal Matrices

Suppose that an indefinite symmetric tridiagonal matrix permits triangular factorization T = LDL ^t . We provide individual condition numbers for the eigenvalues and eigenvectors of T when the parameters in L and D suffer small relative perturbations. When there is element growth in the factorization, then some pairs may be robust while others are sensitive. A 4 × 4 example shows the limitations of standard multiplicative perturbation theory and the efficacy of our new condition numbers.     

Full-rank block LDL ? decomposition and the inverses of n×n block matrices

Full-rank block LDL ^? decomposition of a Hermitian n×n block matrix A is examined, where the iterative procedure evaluating the sub-matrices appearing in L and D is provided. This factorization is used to evaluate the inverse and Moore-Penrose inverse of a Hermitian n×n block matrix. The method for the calculation of the Moore-Penrose inverse of an arbitrary 2×2 block matrix is also provided. Therefore, matrix products A ^? A and AA ^? and the corresponding full-rank block LDL ^? factorizations are observed. Also, a simple explicit formulae calculating the solution vector components of the normal system of equations is stated, where the LDL ^? decomposition of the system matrix is done.     

How Competitive is the ADI for Tensor Structured Equations?

In [1] we presented an extension of the alternating direction implicit (ADI) method for the solution of Lyapunov equations (1) to higher dimensional problems. The vectorized form of the Lyapunov equation is We considered the generalization of this equation of the form (2) The tensor train structure is one possible generalization of the low rank factorization we find in the right hand side of (1). Therefor we assume B to be of tensor train structure. We showed that in analogy to the low rank ADI case the solution X can be generated in tensor train structure, too. Further we provided an algorithm that computes X using a generalization of the ADI method. Here we compare our new tensor ADI method with an density matrix renormalization group (DMRG) solver for tensor train matrix equations and with matrix equation solvers to investigate the competitiveness of our new solver. (© 2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Solution of sparse rectangular systems using LSQR and CRAIG

We examine two iterative methods for solving rectangular systems of linear equations: LSQR for over-determined systemsAx b, and Craig's method for under-determined systemsAx = b. By including regularization, we extend Craig's method to incompatible systems, and observe that it solves the same damped least-squares problems as LSQR. The methods may therefore be compared on rectangular systems of arbitrary shape.Various methods for symmetric and unsymmetric systems are reviewed to illustrate the parallels. We see that the extension of Craig's method closes a gap in existing theory. However, LSQR is more economical on regularized problems and appears to be more reliable if the residual is not small.In passing, we analyze a scaled augmented system associated with regularized problems. A bound on the condition number suggests a promising direct method for sparse equations and least-squares problems, based on indefiniteLDL ^T factors of the augmented matrix.Dedicated to Professor Åke Björck in honor of his 60th birthdayPresented at the 12th Householder Symposium on Numerical Algebra, Lake Arrowhead, California, June 1993.Partially supported by Department of Energy grant DE-FG03-92ER25117, National Science Foundation grant DMI-9204208, and Office of Naval Research grant N00014-90-J-1242.     

Nonconvex quadratic optimization on a parallelepiped

An approximate combinatorial method for solving optimization problems is used to search for a global maximum of a quadratic function on a parallelepiped. Approximating functions in this method are majorants of an object function which are defined on subsets of a parallelepiped of admissible solutions. The method is based on a diagonal or block-diagonal LDL ^T -factorization of a matrix of an object function.     

Equations for the Missing Boundary Values in the Hamiltonian Formulation of Optimal Control Problems

Partial differential equations for the unknown final state and initial costate arising in the Hamiltonian formulation of regular optimal control problems with a quadratic final penalty are found. It is shown that the missing boundary conditions for Hamilton’s canonical ordinary differential equations satisfy a system of first-order quasilinear vector partial differential equations (PDEs), when the functional dependence of the H-optimal control in phase-space variables is explicitly known. Their solutions are computed in the context of nonlinear systems with ℝ ⁿ -valued states. No special restrictions are imposed on the form of the Lagrangian cost term. Having calculated the initial values of the costates, the optimal control can then be constructed from on-line integration of the corresponding 2n-dimensional Hamilton ordinary differential equations (ODEs). The off-line procedure requires finding two auxiliary n×n matrices that generalize those appearing in the solution of the differential Riccati equation (DRE) associated with the linear-quadratic regulator (LQR) problem. In all equations, the independent variables are the finite time-horizon duration T and the final-penalty matrix coefficient S, so their solutions give information on a whole two-parameter family of control problems, which can be used for design purposes. The mathematical treatment takes advantage from the symplectic structure of the Hamiltonian formalism, which allows one to reformulate Bellman’s conjectures concerning the “invariant-embedding” methodology for two-point boundary-value problems. Results for LQR problems are tested against solutions of the associated differential Riccati equation, and the attributes of the two approaches are illustrated and discussed. Also, nonlinear problems are numerically solved and compared against those obtained by using shooting techniques.     

Extending the eigCG algorithm to nonsymmetric Lanczos for linear systems with multiple right‐hand sides

The technique that was used to build the eigCG algorithm for sparse symmetric linear systems is extended to the nonsymmetric case using the BiCG algorithm. We show that, similar to the symmetric case, we can build an algorithm that is capable of computing a few smallest magnitude eigenvalues and their corresponding left and right eigenvectors of a nonsymmetric matrix using only a small window of the BiCG residuals while simultaneously solving a linear system with that matrix. For a system with multiple right‐hand sides, we give an algorithm that computes incrementally more eigenvalues while solving the first few systems and then uses the computed eigenvectors to deflate BiCGStab for the remaining systems. Our experiments on various test problems, including Lattice QCD, show the remarkable ability of eigBiCG to compute spectral approximations with accuracy comparable with that of the unrestarted, nonsymmetric Lanczos. Furthermore, our incremental eigBiCG followed by appropriately restarted and deflated BiCGStab provides a competitive method for systems with multiple right‐hand sides. Copyright © 2013 John Wiley & Sons, Ltd.     

A remark on uniqueness for some ODE

Abstract We present a uniqueness result for the Cauchy problem associated to a particular type of ordinary differential equation (ODE), under the only assumption of continuity of the right hand side at the initial point. Keywords: Polar coordinates, Tangent vector, Inner product Mathematics Subject Classification (2000): 34A12     

On the Wl q-Regularity of Solutions to Systems of Differential Equations in the Case When the Equations Are Constructed from Discontinuous Functions

Some solution, final in a sense from the standpoint of the theory of Sobolev spaces, is obtained to the problem of regularity of solutions to a system of (generally) nonlinear partial differential equations in the case when the system is locally close to elliptic systems of linear equations with constant coefficients. The main consequences of this result are Theorems 5 and 8. According to the first of them, the higher derivatives of an elliptic C ^l -smooth solution to a system of lth-order nonlinear partial differential equations constructed from C ^l -smooth functions meet the local Hoelder condition with every exponent , 0<<1. Theorem 8 claims that if a system of linear partial differential equations of order l with measurable coefficients and right-hand sides is uniformly elliptic then, under the hypothesis of a (sufficiently) slow variation of its leading coefficients, the degree of local integrability of lth-order partial derivatives of every W ^l _q,loc-solution, q>1, to the system coincides with the degree of local integrability of lower coefficients and right-hand sides.     

Diagonally compensated reduction and related preconditioning methods

When solving linear algebraic equations with large and sparse coefficient matrices, arising, for instance, from the discretization of partial differential equations, it is quite common to use preconditioning to accelerate the convergence of a basic iterative scheme. Incomplete factorizations and sparse approximate inverses can provide efficient preconditioning methods but their existence and convergence theory is based mostly on M-matrices (H-matrices). In some application areas, however, the arising coefficient matrices are not H-matrices. This is the case, for instance, when higher-order finite element approximations are used, which is typical for structural mechanics problems. We show that modification of a symmetric, positive definite matrix by reduction of positive offdiagonal entries and diagonal compensation of them leads to an M-matrix. This diagonally compensated reduction can take place in the whole matrix or only at the current pivot block in a recursive incomplete factorization method. Applications for constructing preconditioning matrices for finite element matrices are described.