A nested Krylov subspace method to compute the sign function of large complex matrices

We present an acceleration of the well-established Krylov–Ritz methods to compute the sign function of large complex matrices, as needed in lattice QCD simulations involving the overlap Dirac operator at both zero and nonzero baryon density. Krylov–Ritz methods approximate the sign function using a projection on a Krylov subspace. To achieve a high accuracy this subspace must be taken quite large, which makes the method too costly. The new idea is to make a further projection on an even smaller, nested Krylov subspace. If additionally an intermediate preconditioning step is applied, this projection can be performed without affecting the accuracy of the approximation, and a substantial gain in efficiency is achieved for both Hermitian and non-Hermitian matrices. The numerical efficiency of the method is demonstrated on lattice configurations of sizes ranging from 4⁴ to 10⁴, and the new results are compared with those obtained with rational approximation methods.     

An iterative method to compute the sign function of a non-Hermitian matrix and its application to the overlap Dirac operator at nonzero chemical potential

The overlap Dirac operator in lattice QCD requires the computation of the sign function of a matrix. While this matrix is usually Hermitian, it becomes non-Hermitian in the presence of a quark chemical potential. We show how the action of the sign function of a non-Hermitian matrix on an arbitrary vector can be computed efficiently on large lattices by an iterative method. A Krylov subspace approximation based on the Arnoldi algorithm is described for the evaluation of a generic matrix function. The efficiency of the method is spoiled when the matrix has eigenvalues close to a function discontinuity. This is cured by adding a small number of critical eigenvectors to the Krylov subspace, for which we propose two different deflation schemes. The ensuing modified Arnoldi method is then applied to the sign function, which has a discontinuity along the imaginary axis. The numerical results clearly show the improved efficiency of the method. Our modification is particularly effective when the action of the sign function of the same matrix has to be computed many times on different vectors, e.g., if the overlap Dirac operator is inverted using an iterative method.     

A new approach for solving singular systems in topology optimization using Krylov subspace methods

In topology optimization, elements without any contribution to the improvement of the objective function vanish by decrease of density of the design parameter. This easily causes a singular stiffness matrix. To avoid the numerical breakdown caused by this singularity, conventional optimization techniques employ additional procedures. These additional procedures, however, raise some problems. On the other hand, convergence of Krylov subspace methods for singular systems have been studied recently. Through subsequent studies, it has been revealed that the conjugate gradient method (CGM) does not converge to the local optimal solution in some singular systems but in those satisfying certain condition, while the conjugate residual method (CRM) yields converged solutions in any singular systems. In this article, we show that a local optimal solution for topology optimization is obtained by using the CRM and the CGM as a solver of the equilibrium equation in the structural analysis, even if the stiffness matrix becomes singular. Moreover, we prove that the CGM, without any additional procedures, realizes convergence to a local optimal solution in that case. Computer simulation shows that the CGM gives almost the same solutions obtained by the CRM in the case of the two-bar truss problem.     

Numerical methods for integrating chemical kinetic rate equations

In this paper we present a numerical method which is suitable for the integration of chemical rate equations. These equations are normally extremely stiff due to large differences in the kinetic rate coefficients. The method takes advantage of the fact that the Jacobian matrix is readily obtainable for this type of problem. Stability analysis will also be discussed in a general framework.     

Optimized Schwarz method without overlap for the gravitational potential equation on cluster of graphics processing unit

Many engineering and scientific problems need to solve boundary value problems for partial differential equations or systems of them. For most cases, to obtain the solution with desired precision and in acceptable time, the only practical way is to harness the power of parallel processing. In this paper, we present some effective applications of parallel processing based on hybrid CPU/GPU domain decomposition method. Within the family of domain decomposition methods, the so-called optimized Schwarz methods have proven to have good convergence behaviour compared to classical Schwarz methods. The price for this feature is the need to transfer more physical information between subdomain interfaces. For solving large systems of linear algebraic equations resulting from the finite element discretization of the subproblem for each subdomain, Krylov method is often a good choice. Since the overall efficiency of such methods depends on effective calculation of sparse matrix–vector product, approaches that use graphics processing unit (GPU) instead of central processing unit (CPU) for such task look very promising. In this paper, we discuss effective implementation of algebraic operations for iterative Krylov methods on GPU. In order to ensure good performance for the non-overlapping Schwarz method, we propose to use optimized conditions obtained by a stochastic technique based on the covariance matrix adaptation evolution strategy. The performance, robustness, and accuracy of the proposed approach are demonstrated for the solution of the gravitational potential equation for the data acquired from the geological survey of Chicxulub crater.