Recursion relations for the extended Krylov subspace method

The extended Krylov subspace method and orthogonal Laurent polynomials

Modern power systems have experienced a significant increase in its complexity, and extremely large-scale system models have to be addressed in the study of stability analysis and control. Model reduction is a technique for developing an approximate system model with lower dimensions that shares similar properties to the original system. This work proposes a computationally efficient approach for linear system model reduction in large-scale power systems that is based on the balanced truncation method using an extended Krylov subspace technique. Key algorithm improvements, including sparsity handling of the linearized power system models and efficient computational techniques for solving dual Lyapunov equations, are discussed in detail. In addition, application of the proposed model reduction is extended to unstable systems, and a power system stabilizer (PSS) parameter optimization method based on reduced system models is used to validate the effectiveness of the proposed method. Case studies show that the resulting reduced models preserve the unique characteristics of the original high-dimensional models, including time and frequency domain responses as well as eigenvalues. Numerical results based on test systems with up to 12 $\thinspace$ 685 buses demonstrate the computational efficiency and validity of the proposed approach.

<abstract> By utilizing an inner-outer iteration strategy, a shift-splitting (SS) iteration method to solve a class of large sparse linear matrix equation $ AXB = C $ is proposed in this work. Two convergence theorems for differential forms are studied in depth. Moreover, the quasi-optimal parameters which minimize the upper bound for the spectral radius of SS iteration matrix are given. Two numerical examples illustrate the high-efficiency of SS iteration method, especially when coefficient matrices are ill-conditioned. </abstract>

A projected preconditioned conjugate gradient algorithm for computing many extreme eigenpairs of a Hermitian matrix

The representation of an integral operator in wavelet bases is a large numerical sparse matrix. Thus, a suitable approach is needed that converts the matrix from full storage mode into sparse storage mode, which not only saves the storage, but also saves the execution time. In this paper, the storage of this kind of matrix is discussed in the symmetric case. First, we present a standard storage approach. Then, a modified one is proposed, which is based on a special decomposition of the symmetric matrix. When iterative solvers are used to solve the corresponding linear system, the main cost is the multiplication of the symmetric sparse matrix and a given vector. The modified sparse storage mode is optimized for this purpose and reduces nearly half of the number of multiplications.

We introduce an economical Gram--Schmidt orthogonalization on the extended Krylov subspace originated by actions of a symmetric matrix and its inverse. An error bound for a family of problems arising from the elliptic method of lines is derived. The bound shows that, for the same approximation quality, the diagonal variant of the extended subspaces requires about the square root of the dimension of the standard Krylov subspaces using only positive or negative matrix powers. An example of an application to the solution of a 2.5-D elliptic problem attests to the computational efficiency of the method for large-scale problems.

We consider the computation of a few eigenvectors and corresponding eigenvalues of a large sparse nonsymmetric matrix using shift-invert Arnoldi's method with and without implicit restarts. For the inner iterations we use preconditioned GMRES as the inexact iterative solver. The costs of the solves are measured by the number of inner iterations needed by the iterative solver at each outer step of the algorithm. We first extend the relaxation strategy developed by Simoncini [SIAM J. Numer. Anal., 43 (2005), pp. 1155–1174] to implicitly restarted Arnoldi's method, which yields an improvement in the overall costs of the method. Secondly, we apply a new preconditioning strategy to the inner solver. We show that small rank changes to the preconditioner can produce significant savings in the total number of iterations. The combination of the new preconditioner with the relaxation strategy in implicitly restarted Arnoldi produces enhancement in the overall costs of around 50 percent in the examples considered here. Numerical experiments illustrate the theory throughout the paper.

We describe a procedure for determining a few of the largest singular values of a large sparse matrix. The method by Golub and Kent which uses the method of modified moments for estimating the eigenvalues of operators used in iterative methods for the solution of linear systems of equations is appropriately modified in order to generate a sequence of bidiagonal matrices whose singular values approximate those of the original sparse matrix. A simple Lanczos recursion is proposed for determining the corresponding left and right singular vectors. The potential asynchronous computation of the bidiagonal matrices using modified moments with the iterations of an adapted Chebyshev semi-iterative (CSI) method is an attractive feature for parallel computers. Comparisons in efficiency and accuracy with an appropriate Lanczos algorithm (with selective re-orthogonalization) are presented on large sparse (rectangular) matrices arising from applications such as information retrieval and seismic reflection tomography. This procedure is essentially motivated by the theory of moments and Gauss quadrature.

In this paper, we present new error bounds for the Lanczos method and the shift-and-invert Lanczos method for computing $e^{-\tau A} v$ for a large sparse symmetric positive semidefinite matrix $A$. Compared with the existing error analysis for these methods, our bounds relate the convergence to the condition numbers of the matrix that generates the Krylov subspace. In particular, we show that the Lanczos method will converge rapidly if the matrix $A$ is well-conditioned, regardless of what the norm of $\tau A$ is. Numerical examples are given to demonstrate the theoretical bounds.

This paper discusses the design and development of a code to calculate the eigenvalues of a large sparse real unsymmetric matrix that are the rightmost, leftmost, or are of the largest modulus. A subspace iteration algorithm is used to compute a sequence of sets of vectors that converge to an orthonormal basis for the invariant subspace corresponding to the required eigenvalues. This algorithm is combined with Chebychev acceleration if the rightmost or leftmost eigenvalues are sought, or if the eigenvalues of largest modulus are known to be the rightmost or leftmost eigenvalues. An option exists for computing the corresponding eigenvectors. The code does not need the matrix explicitly since it only requires the user to multiply sets of vectors by the matrix. Sophisticated and novel iteration controls, stopping criteria, and restart facilities are provided. The code is shown to be efficient and competitive on a range of test problems.

The Lanczos algorithm is among the most frequently used iterative techniques for computing a few dominant eigenvalues of a large sparse non-symmetric matrix. At the same time, it serves as a building block within biconjugate gradient (BiCG) and quasi-minimal residual (QMR) methods for solving large sparse non-symmetric systems of linear equations. It is well known that, when implemented on distributed-memory computers with a huge number of processes, the synchronization time spent on computing dot products increasingly limits the parallel scalability. Therefore, we propose synchronization-reducing variants of the Lanczos, as well as BiCG and QMR methods, in an attempt to mitigate these negative performance effects. These so-called s-step algorithms are based on grouping dot products for joint execution and replacing time-consuming matrix operations by efficient vector recurrences. The purpose of this paper is to provide a rigorous derivation of the recurrences for the s-step Lanczos algorithm, introduce s-step BiCG and QMR variants, and compare the parallel performance of these new s-step versions with previous algorithms.

In many scientific computing applications, sparse Cholesky factorization is used to solve large sparse linear equations in distributed environment. GPU computing is a new way to solve the problem. However, sparse Cholesky factorization on GPU is hardly to achieve excellent performance due to the structure irregularity of matrix and the low GPU resource utilization. A hybrid CPU-GPU implementation of sparse Cholesky factorization is proposed based on multifrontal method. A large sparse coefficient matrix is decomposed into a series of small dense matrices (frontal matrices) in the method, and then multiple GEMM (General Matrix-matrix Multiplication) operations are computed on them. GEMMs are the main operations in sparse Cholesky factorization, but they are hardly to perform better in parallel on GPU. In order to improve the performance, the scheme of multiple task queues is adopted to perform multiple GEMMs parallelized with multifrontal method; all GEMM tasks are scheduled dynamically on GPU and CPU based on computation scales for load balance and computing-time reduction. Experimental results show that the approach can outperform the implementations of cuBLAS, achieving up to 1.98× speedup on GTX460 (Fermi micro-architecture) and 3.06× speedup on K20m (Kepler micro-architecture), respectively.

Overhead transmission lines are subjected to long-term vibration loads, which can easily produce cracks and broken strands, affecting the safe operation of the power grid. Steel-core aluminum stranded conductor is commonly used in overhead transmission lines, conductor strand structure is complex, therefore, the need to "nominal" bending stiffness to characterize its bending capacity, that is, the equivalent bending stiffness of overhead transmission lines. When the equivalent bending stiffness of the conductor suspension point and near the exit of the fixture changes due to fatigue, it will amplify the vibration level of the place, which in turn affects the arc sag and air gap of the conductor. The paper for the occurrence of transmission line breeze vibration characteristics, analysis of the equivalent bending stiffness near the suspension point of the tension tower conductor, the establishment of a 3D model of the conductor, using a large sparse solution matrix as a calculation method, the calculation of the equivalent bending stiffness of the conductor under different constraint length,on this basis, further analysis of the number of broken strands andbroken strand displacement on the equivalent stiffness of the impact. The calculation results show that the equivalent bending stiffness of the conductor decreases with the increase of the constraint length, decreases with the increase of the number of broken strands, and shows a trend of first decreasing and then increasing with the increase of the displacement of broken strands.

ABSTRACTAn accurate and efficient 3D finite‐difference (FD) forward algorithm for DC resistivity modelling is developed. In general, the most time‐consuming part of FD calculation is to solve large sets of linear equations: Ax=b, where A is a large sparse band symmetric matrix. The direct method using complete Choleski decomposition is quite slow and requires much more computer storage. We have introduced a row‐indexed sparse storage mode to store the coefficient matrix A and an incomplete Choleski conjugate‐gradient (ICCG) method to solve the large linear systems. By taking advantage of the matrix symmetry and sparsity, the ICCG method converges much more quickly and requires much less computer storage. It takes approximately 15 s on a 533 MHz Pentium computer for a grid with 46 020 nodes, which is approximately 700 times faster than the direct method and 2.5 times faster than the symmetric successive over‐relaxation (SSOR) conjugate‐gradient method. Compared with 3D finite‐element resistivity modelling with the improved ICCG solver, our algorithm is more efficient in terms of number of iterations and computer time. In addition, we solve for the secondary potential in 3D DC resistivity modelling by a simple manipulation of the FD equations. Two numerical examples of a two‐layered model and a vertical contact show that the method can achieve much higher accuracy than solving for the total potential directly with the same grid nodes. In addition, a 3D cubic body is simulated, for which the dipole–dipole apparent resistivities agree well with the results obtained with the finite‐element and integral‐equation methods. In conclusion, the combination of several techniques provides a rapid and accurate 3D FD forward modelling method which is fundamental to 3D resistivity inversion.

An implicit non-linear lower-upper symmetric Gauss-Seidel (LU-SGS) solution algorithm has been developed for a high-order spectral difference Navier-Stokes solver on unstructured hexahedral grids. The non-linear LU-SGS solver is preconditioned by the block element matrix, and the system of equations is then solved with a LU decomposition approach. The large sparse Jacobian matrix is computed numerically, resulting in extremely simple operations for arbitrarily complex residual operators. Several viscous test cases were performed to evaluate the performance. The implicit solver has shown speedup of 1 to 2 orders of magnitude over the multi-stage Runge-Kutta time integration scheme.