Linearly scaling direct method for accurately inverting sparse banded matrices

In many problems in Computational Physics and Chemistry, one finds a special kind of sparse matrices, termed "banded matrices". These matrices, which are defined as having non-zero entries only within a given distance from the main diagonal, need often to be inverted in order to solve the associated linear system of equations. In this work, we introduce a new O(n) algorithm for solving such a system, being n X n the size of the matrix. We produce the analytical recursive expressions that allow to directly obtain the solution, as well as the pseudocode for its computer implementation. Moreover, we review the different options for possibly parallelizing the method, we describe the extension to deal with matrices that are banded plus a small number of non-zero entries outside the band, and we use the same ideas to produce a method for obtaining the full inverse matrix. Finally, we show that the New Algorithm is competitive, both in accuracy and in numerical efficiency, when compared to a standard method based in Gaussian elimination. We do this using sets of large random banded matrices, as well as the ones that appear when one tries to solve the 1D Poisson equation by finite differences.Comment: 24 pages, 5 figures, submitted to J. Comp. Phy

Exact and efficient calculation of Lagrange multipliers in constrained biological polymers: Proteins and nucleic acids as example cases

In order to accelerate molecular dynamics simulations, it is very common to impose holonomic constraints on their hardest degrees of freedom. In this way, the time step used to integrate the equations of motion can be increased, thus allowing, in principle, to reach longer total simulation times. The imposition of such constraints results in an aditional set of Nc equations (the equations of constraint) and unknowns (their associated Lagrange multipliers), that must be solved in one way or another at each time step of the dynamics. In this work it is shown that, due to the essentially linear structure of typical biological polymers, such as nucleic acids or proteins, the algebraic equations that need to be solved involve a matrix which is banded if the constraints are indexed in a clever way. This allows to obtain the Lagrange multipliers through a non-iterative procedure, which can be considered exact up to machine precision, and which takes O(Nc) operations, instead of the usual O(Nc3) for generic molecular systems. We develop the formalism, and describe the appropriate indexing for a number of model molecules and also for alkanes, proteins and DNA. Finally, we provide a numerical example of the technique in a series of polyalanine peptides of different lengths using the AMBER molecular dynamics package.Comment: 29 pages, 10 figures, 1 tabl

Performance of Parallel Approximate Ideal Restriction Multigrid for Transport Applications

Algebraic multigrid (AMG) methods have been widely used to solve systems arising from the discretization of elliptic partial differential equations. In serial, AMG algorithms scale linearly with problem size. In parallel, communication costs scale logarithmically with the number of processors. Recently, a classical AMG method based on approximate ideal restriction (AIR) was developed for nonsymmetric matrices. AIR has already been shown to be effective for solving the linear systems arising from upwind discontinuous Galerkin (DG) finite element discretization of advection-diffusion problems, including the hyperbolic limit of pure advection. A new parallel version of AIR, pAIR, has been implemented in the hypre library. In this thesis, pAIR is tested for use solving the source iteration equations of the SN approximations to the transport equation. The performance is investigated with various meshes in two and three dimensions. Detailed profiling of parallel performance is also conducted to identify the most important areas for algorithm improvements. An improvement to the Local Ideal Approximate Restriction algorithm is introduced and discussed. Weak scaling results to 4,096 processors are presented. These results show total solve growing logarithmically with the number of processors used. Importantly, this result is shown on both uniform grids and unstructured grids in three dimensions. The unstructured mesh did not include reentrant cells

Incorporating Krylov Subspace Methods in the ETDRK4 Scheme

A modification of the (2,2)-Pade algorithm developed by Wade et al. for implementing the exponential time differencing fourth order Runge-Kutta (ETDRK4) method is introduced. The main computational difficulty in implementing the ETDRK4 method is the required approximation to the matrix exponential. Wade et al. use the fourth order (2,2)-Pade approximant in their algorithm and in this thesis we incorporate Krylov subspace methods in an attempt to improve efficiency. A background of Krylov subspace methods is provided and we describe how they are used in approximating the matrix exponential and how to implement them into the ETDRK4 method. The (2,2)-Pade and Krylov subspace algorithms are compared in solving the one and two dimensional Allen-Cahn equation with the ETDRK4 scheme. We find that in two dimensions, the Krylov subspace algorithm is faster, provided we have a spatial discretization that produces a symmetric matrix