Shift-splitting iteration methods for a class of large sparse linear matrix equations

Finite-difference time-domain (FDTD) methods have long been the most popular solution for transient electromagnetic (TEM) modeling. Compared with FDTD methods, finiteelement time-domain (FETD) methods have the potential to reduce the number of unknowns and the number of time steps. However, with a requirement to solve an irregular and large sparse matrix equation at every time step, it is generally difficult for FETD methods to take advantage of massively parallel computers. In our work, we propose a parallel FETD solution for TEM diffusion phenomena. We design and implement a parallel iterative solver, which includes a preconditioner that is customized for the FETD method to achieve both fast convergence and good scalability. Using a seabed model, our experiments show that our parallel FETD solution provides the same accuracy as the FDTD methods. On the performance side, using 8 CPU cores, the iterative solver part runs 5 times faster than the serial verison, while the preconditioner part runs 3 times faster. In total, our parallel FETD method runs 3.5 times faster than the serial verison when simulating a seabed model with around 100,000 unknowns.

A transient electromagnetic vector finite-element forward-modeling method based on 3D anisotropic resistivity constrained by well-logging and seismic data is developed. This method integrates detailed stratigraphic division provided by well-logging data and structural information from seismic interpretation to establish a structural model that reflects the actual geologic conditions. Based on the structural model, a kriging interpolation method is applied to the logging resistivity data to obtain the anisotropic resistivity distribution information of each layer in the reservoir, thereby establishing a 3D anisotropic resistivity model. To ensure the smoothness of the model, we use a variogram function based on a spherical model to fit the discrete interpolation points. After the resistivity model is spatially discretized by locally refined unstructured tetrahedral meshes and temporally discretized by a nonuniform time step, the large sparse matrix equation is solved using the PARDISO direct solver to obtain the forward response. We apply this method to layered models of isotropic and vertical transversely isotropic (VTI) media and verify the accuracy of the forward simulation method by comparative analysis with quasi-analytical solutions. Doubling the nonuniform time step in our forward simulations does not affect the accuracy of the solution. Applying our method to China’s D oilfield, the electromagnetic response at observation points generally matches traditional methods well. Comparative analysis of the forward results between our method and the model without geologic data or homogeneous resistivity shows that we achieve better fitting with the field measurement data and provide more detailed electrical properties.

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.

We derive a new upper bound on the minimum rank of matrices belonging to an affine slice of the positive semidefinite cone, when the affine slice is defined according to a system of sparse linear matrix equations. It is shown that a feasible matrix whose rank is no greater than said bound can be computed in polynomial time. The bound depends on both the number of linear matrix equations and their underlying sparsity pattern. For certain problem families, this bound is shown to improve upon well known bounds in the literature. Several examples are provided to illustrate the efficacy of this bound.

In this letter, an efficient vertex-based domain decomposition technique is introduced into the unconditionally stable finite-difference time-domain (FDTD) method based on weighted Laguerre polynomials (WLPs). In the two dimensional (2-D) transversal electric (TE) case, the vertex-based partitioning scheme is implemented and the magnetic field is updated implicitly. After the whole computational domain is decomposed into multiple subdomains with the vertex-based partitioning scheme, the large sparse matrix equation generated by the 2-D WLP-FDTD method is transformed into some independent small ones. A numerical example is presented to validate the accuracy and efficiency of the proposed method.

In this paper, an efficient domain decomposition technique is introduced into the unconditionally stable finite-difference time-domain (FDTD) method based on weighted Laguerre polynomials to solve two-dimensional (2-D) electromagnetic scattering problems. The whole computational space is decomposed into multiple subdomains where there is no direct field coupling between any two different subdomains. For the large sparse matrix equation generated by the implicit scheme, the domain decomposition technique transforms this large scale equation into some independent smaller equations. With the total-field/scattered-field boundary and Mur's second-order absorbing boundary condition, the radar cross sections of two 2-D structures are calculated. The numerical examples verify the accuracy and efficiency of the proposed method.

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, an efficient algorithm for implementing Crank-Nicolson scheme in the finite-element time-domain (FETD) method is presented. Based on a direct discretization of the first-order coupled Maxwell curl equations, this algorithm employs edge elements (Whitney 1-form) to expand the electric field and face elements (Whitney 2-form) for the magnetic field. Since the curl of an edge-element is the linear combination of those face elements whose faces contain the given edge, only the Maxwell-Ampere equation composes a sparse linear matrix equation for the electric field update; the Maxwell-Faraday equation is explicit. The Crank-Nicolson scheme is implemented leading to an unconditionally stable vector FETD method and the matrix inverse is not required to be computed explicitly. Therefore, only one matrix equation is required to be solved at each time step. Numerical results demonstrate that the proposed method is efficient when compared with the conventional leap-frog mixed FETD method and the Crank-Nicolson FDTD method.

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.

American Institute of Mining, Metallurgical, and Petroleum Engineers, Inc. Abstract An iterative method is proposed for the inversion of sparse band-structured matrices of the type that are common in numerical reservoir simulators. The method utilizes orthogonalizations and minimizations to achieve a fast convergence rate. Tests have shown that the method is highly competitive compared with other iterative techniques. The rate of convergence is insensitive to the use of iteration parameters, non-symmetry of the matrix, and ratios of off-diagonal bands in symmetrical matrices. 1. Introduction A time-consuming part of the calculation in numerical reservoir simulators is the inversion of a large set of simultaneous linear equations. (1) For example when effecting a solution by the Implicit Potential - Explicit Saturation method x would represent the pressure or potential field, whilst A would be a matrix of transmissibility coefficients describing the interflow of fluids between grid blocks. In the common finite difference approaches, A is a sparse banded matrix, as depicted in figure 1 for a five-point finite difference representation of a two dimensional system. Often, the diagonal element of A is approximately equal to the negative sum of the off-diagonal elements in a given row. Because of the sparseness of A, iterative methods appear at first sight attractive compared with direct methods. However, specialized direct routines with increased efficiency have become available in recent years. Of the iterative methods SIP (strongly implicit procedure) is possibly the most widely used. SIP has the possibly the most widely used. SIP has the disadvantage that its convergence rate depends strongly on a set of iteration parameters whose optimal values can be difficult to find. A new highly competitive iterative method for symmetric matrices, a minimization process conceptually based on the conjugate-gradient technique, has lately been developed. The numerical examples and proofs quoted in reference 4 indicate that minimization processes can be expected to offer the following advantages: Convergence more readily guaranteed. No need for iteration parameters. Insensitivity to number of equations. Insensitivity to transmissibility ratios (the ratios of the off-diagonal bands in figure 1).

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.

In this paper, a modified Hermitian and skew-Hermitian splitting (MHSS) iteration method for solving the complex linear matrix equation $AXB = C$ has been presented. As the theoretical analysis shows, the MHSS iteration method will converge under certain conditions. Each iteration in this method requires the solution of four linear matrix equations with real symmetric positive definite coefficient matrices, although the original coefficient matrices are complex and non-Hermitian. In addition, the optimal parameter of the new iteration method is proposed. Numerical results show that MHSS iteration method is efficient and robust.

Finite-Difference Frequency-Domain methods (FDFD) require solution of large linear systems of equations. These large systems are represented by matrix equations including highly sparse coefficient matrices, and they can often only be solved by using iterative methods. This paper presents an algorithm in which the matrix-equation solution approach in an iterative method is replaced by a multi-step solution process. Instead of using a coefficient matrix, the coefficients in the FDFD formulations are kept as three-dimensional arrays, and they are treated as operators. The algorithm is used together with the Bi-Conjugate Gradients Stabilized (BICGSTAB) method. This is applied to a three-dimensional FDFD method to solve for scattering from dielectric objects. It is also applied to two other FDFD methods (a single-grid and a double-grid FDFD) to solve for scattering from chiral objects. It has been shown that the presented algorithm effectively reduces the solution time and memory requirements.

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.

Multivariate regression is known as a multivariate extension of multiple regression, which explain/predict the variations in multiple dependent variables by multiple independent variables. Recently, various procedures for Sparse Multivariate Regression (SMR) have been proposed, in which a sparse regression coefficient matrix (having a number of zero elements) is obtained aiming to facilitate its interpretation. The procedures for SMR can be classified into the following two types; penalized and cardinality-constrained procedures. In them, the resulting number of zeros in the regression coefficient matrix is controlled/constrained by a prespecified penalty parameter or cardinality value. In this research, we propose another approach for SMR, referred to as Layered Multivariate Regression (LMR). In LMR, the regression coefficient matrix is assumed to be a sum of several sparse matrices, which is called layer. Therefore, the sparseness of the resulting coefficient matrix is controlled by how many layers are used. In LMR, k-th layer can be viewed as the coefficient matrix in the regression of a partial residual (i.e., the residual for all but k-th layer) on independent variables, and thus the variance explained by LMR gets closer to that for the unconstrained regression as the number of layers increases. We present an alternating least squares algorithm for LMR and a procedure for determining how many layers should be used. LMR is assessed in a simulation study and illustrated with a real data example. As an application of LMR, procedures for sparse estimation in some multivariate analysis techniques (e.g., principal component analysis) are also presented.