We perform a spectral analysis of the preconditioned Hermitian/skew-Hermitian splitting (PHSS) method applied to multilevel block Toeplitz linear systems in which the coefficient matrix *T*_{n}(*f*) is associated with a Lebesgue integrable matrix-valued function *f*. When the preconditioner is chosen as a Hermitian positive definite multilevel block Toeplitz matrix *T*_{n}(*g*), the resulting sequence of PHSS iteration matrices *M*_{n} belongs to the generalized locally Toeplitz class. In this case, we are able to compute the symbol *ϕ*(*f*,*g*) describing the asymptotic eigenvalue distribution of *M*_{n}when *n**∞* and the matrix size diverges. By minimizing the infinity norm of the spectral radius of the symbol *ϕ*(*f*,*g*), we are also able to identify effective PHSS preconditioners *T*_{n}(*g*) for the matrix *T*_{n}(*f*). A number of numerical experiments are presented and commented, showing that the theoretical results are confirmed and that the spectral analysis leads to efficient PHSS methods. Copyright © 2015 John Wiley & Sons, Ltd.

We present an algebraic structured preconditioner for the iterative solution of large sparse linear systems. The preconditioner is based on a multifrontal variant of sparse LU factorization used with nested dissection ordering. Multifrontal factorization amounts to a partial factorization of a sequence of logically dense frontal matrices, and the preconditioner is obtained if structured factorization is used instead. This latter exploits the presence of low numerical rank in some off-diagonal blocks of the frontal matrices. An algebraic procedure is presented that allows to identify the hierarchy of the off-diagonal blocks with low numerical rank based on the sparsity of the system matrix. This procedure is motivated by a model problem analysis, yet numerical experiments show that it is successful beyond the model problem scope. Further aspects relevant for the algebraic structured preconditioner are discussed and illustrated with numerical experiments. The preconditioner is also compared with other solvers, including the corresponding direct solver. Copyright © 2015 John Wiley & Sons, Ltd.

In various applications, for instance, in the detection of a Hopf bifurcation or in solving separable boundary value problems using the two-parameter eigenvalue problem, one has to solve a generalized eigenvalue problem with 2 × 2 operator determinants of the form

We present efficient methods that can be used to compute a small subset of the eigenvalues. For full matrices of moderate size, we propose either the standard implicitly restarted Arnoldi or Krylov–Schur iteration with shift-and-invert transformation, performed efficiently by solving a Sylvester equation. For large problems, it is more efficient to use subspace iteration based on low-rank approximations of the solution of the Sylvester equation combined with a Krylov–Schur method for the projected problems. Copyright © 2015 John Wiley & Sons, Ltd.

We construct, analyze, and implement SSOR-like preconditioners for non-Hermitian positive definite system of linear equations when its coefficient matrix possesses either a dominant Hermitian part or a dominant skew-Hermitian part. We derive tight bounds for eigenvalues of the preconditioned matrices and obtain convergence rates of the corresponding SSOR-like iteration methods as well as the corresponding preconditioned GMRES iteration methods. Numerical implementations show that Krylov subspace iteration methods such as GMRES, when accelerated by the SSOR-like preconditioners, are efficient solvers for these classes of non-Hermitian positive definite linear systems. Copyright © 2015 John Wiley & Sons, Ltd.

This article discusses a more general and numerically stable Rybicki Press algorithm, which enables inverting and computing determinants of covariance matrices, whose elements are sums of exponentials. The algorithm is true in exact arithmetic and relies on introducing new variables and corresponding equations, thereby converting the matrix into a banded matrix of larger size. Linear complexity banded algorithms for solving linear systems and computing determinants on the larger matrix enable linear complexity algorithms for the initial semi-separable matrix as well. Benchmarks provided illustrate the linear scaling of the algorithm. Copyright © 2015 John Wiley & Sons, Ltd.

In this work, we consider the numerical solution of large nonlinear eigenvalue problems that arise in thermoacoustic simulations involved in the stability analysis of large combustion devices. We briefly introduce the physical modeling that leads to a nonlinear eigenvalue problem that is solved using a nonlinear fixed point iteration scheme. Each step of this nonlinear method requires the solution of a complex non-Hermitian linear eigenvalue problem. We review a set of state-of-the-art eigensolvers and discuss strategies to recycle spectral information from one nonlinear step to the next. More precisely, we consider the Jacobi–Davidson algorithm, the Implicitly Restarted Arnoldi method, the Krylov–Schur solver and its block-variant, and the subspace iteration method with Chebyshev acceleration. On a small test example, we study the relevance of the different approaches and illustrate on a large industrial test case the performance of the parallel solvers best suited to recycle spectral information for large-scale thermoacoustic stability analysis. Copyright © 2015 John Wiley & Sons, Ltd.

Discrete representations of the Helmholtz operator generally give rise to extremely difficult linear systems from an iterative solver perspective. This is due in part to the large oscillatory near null space of the linear system. Typical iterative methods do not effectively reduce error components in the subspace associated with this near null space. Traditional coarse grids used within multilevel solvers also cannot capture these components because of their oscillatory nature. While the shifted Laplacian is a popular method allowing multigrid preconditioners to achieve improved convergence, it is still unsatisfactory for higher-frequency problems. This paper analyzes the shifted Laplacian through polynomial smoothers to show its effect on the spectrum of the discrete Helmholtz operator. This analysis reveals desirable features of the shifted Laplacian preconditioners as well as limitations. Shifted Laplacian techniques are first extended to smoothed aggregation algebraic multigrid. Motivated by analysis, we propose augmenting the algebraic multigrid shifted Laplacian with a two-grid error correction step. This additional step consists of the generalized minimal residual iterations applied to a projected version of the unshifted equations on a single auxiliary coarse grid. Results indicate that augmentation can improve the convergence significantly and can reduce the overall solve time on two-dimensional and three-dimensional problems. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we consider finding a sparse solution of nonnegative least squares problems with a linear equality constraint. We propose a projection-based gradient descent method for solving huge size constrained least squares problems. Traditional Newton-based methods require solving a linear system. However, when the matrix is huge, it is neither practical to store it nor possible to solve it in a reasonable time. We therefore propose a matrix-free iterative scheme for solving the minimizer of the captured problem. This iterative scheme can be explained as a projection-based gradient descent method. In each iteration, a projection operation is performed to ensure the solution is feasible. The projection operation is equivalent to a shrinkage operator, which can guarantee the sparseness of the solution obtained. We show that the objective function is decreasing. We then apply the proposed method to the inverse problem of constructing a probabilistic Boolean network. Numerical examples are then given to illustrate both the efficiency and effectiveness of our proposed method. Copyright © 2015 John Wiley & Sons, Ltd.

A cheap symmetric stiffness-based preconditioning of the Hessian of the dual problem arising from the application of the finite element tearing and interconnecting domain decomposition to the solution of variational inequalities with varying coefficients is proposed. The preconditioning preserves the structure of the inequality constraints and affects both the linear and nonlinear steps, so that it can improve the rate of convergence of the algorithms that exploit the conjugate gradient steps or the gradient projection steps. The bounds on the regular condition number of the Hessian of the preconditioned problem, which are independent of the coefficients, are given. The related stiffness scaling is also considered and analysed. The improvement is demonstrated by numerical experiments including the solution of a contact problem with variationally consistent discretization of the non-penetration conditions. The results are relevant also for linear problems. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents parallel preconditioners and multigrid solvers for solving linear systems of equations arising from stochastic polynomial chaos formulations of the diffusion equation with random coefficients. These preconditioners and solvers are extensions of the preconditioner developed in an earlier paper for strongly coupled systems of elliptic partial differential equations that are norm equivalent to systems that can be factored into an algebraic coupling component and a diagonal differential component. The first preconditioner, which is applied to the norm equivalent system, is obtained by sparsifying the inverse of the algebraic coupling component. This sparsification leads to an efficient method for solving these systems at the large scale, even for problems with large statistical variations in the random coefficients. An extension of this preconditioner leads to stand-alone multigrid methods that can be applied directly to the actual system rather than to the norm equivalent system. These multigrid methods exploit the algebraic/differential factorization of the norm equivalent systems to produce variable-decoupled systems on the coarse levels. Moreover, the structure of these methods allows easy software implementation through re-use of robust high-performance software such as the Hypre library package. Two-grid matrix bounds will be established, and numerical results will be given. Copyright © 2015 John Wiley & Sons, Ltd.

Using the Schur complement of a matrix, we propose a computational framework for performing constrained maximum likelihood estimation in which the unknown parameters can be partitioned into two sets. Under appropriate regularity conditions, the corresponding estimating equations form a non-linear system of equations with constraints. Solving this system is typically accomplished via methods which require computing or estimating a Hessian matrix. We present an alternative algorithm that solves the constrained non-linear system in block coordinate descent fashion. An explicit form for the solution is given. The overall algorithm is shown in numerical studies to be faster than standard methods that either compute or approximate the Hessian as well as the classical Nelder–Mead algorithm. We apply our approach to a motivating problem of evaluating the effectiveness of Road Safety Policies. This includes several numerical studies on simulated data. Copyright © 2015 John Wiley & Sons, Ltd.

Combining the modified matrix–vector equation approach with the technique of Lyapunov majorant function and the Banach fixed point theorem, we obtain improved rigorous perturbation bounds for the LU and QR factorizations with normwise perturbation in the given matrix. Each of the improved rigorous perturbation bounds is a rigorous version of the first-order perturbation bound derived by the matrix–vector equation approach in the literature, and we present their explicit expressions. These bounds are always tighter than those given by Chang and Stehlé in the paper entitled “Rigorous perturbation bounds of some matrix factorizations”. This fact is illustrated by numerical examples. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we study the convergence property of a block improvement method (BIM) for the bi-quadratic polynomial optimization problem over unit spheres. We establish the global convergence of the method generally and establish its linear convergence rate under the second-order sufficient condition. We also extend the BIM to inhomogeneous polynomial optimization problems over unit spheres. Numerical results reported in this paper show that the method is promising. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we study tensors with time dimension which arises in many scientific and engineering applications such as time series gene expression analysis and video analysis. In these applications, we are interested in determining a set of components interacting closely and consistently together over a period of time. The main aim of this paper is to develop a numerical method to compute such constrained CANDECOMP/PARAFAC (CP) decompositions. We make use of the total variation regularization to constrain the time dimension factor in the decomposition in order to obtain a piecewise constant function with respect to time points. The components of the other dimensions corresponding to these time points are closely related. For example, in time series gene expression analysis, a set of genes may regulate a biological process together during a specific time period; in video analysis, a set of image pixels may refer to a foreground object in the video frames. We employ ADMM to solve the resulting optimization problem, and study its convergence. Numerical examples on synthetic and real data sets are used to demonstrate that the proposed total variation based CP decomposition model can provide more accurate and interesting results. Copyright © 2015 John Wiley & Sons, Ltd.

We extend the balancing domain decomposition by constraints (BDDC) method to flows in porous media discretised by mixed-hybrid finite elements with combined mesh dimensions. Such discretisations appear when major geological fractures are modelled by one-dimensional or two-dimensional elements inside three-dimensional domains. In this set-up, the global problem and the substructure problems have a symmetric saddle-point structure, containing a ‘penalty’ block due to the combination of meshes. We show that the problem can be reduced by means of iterative substructuring to an interface problem, which is symmetric and positive definite. The interface problem can thus be solved by conjugate gradients with the BDDC method as a preconditioner. A parallel implementation of this algorithm is incorporated into an existing software package for subsurface flow simulations. We study the performance of the iterative solver on several academic and real-world problems. Numerical experiments illustrate its efficiency and scalability. Copyright © 2015 John Wiley & Sons, Ltd.

Users of social media interact with the network and its users. Each interaction creates network-specific data between the engaged users and the chosen social avenue. Over time, these engagements accumulate to describe the user's social fingerprint, a data trail that encapsulates the user's identity on the network. The agglomeration of this information showcases the user's activity on the social network and establishes a traceable social fingerprint. The focus of this study is to examine the accuracy of user identification via semidiscrete decomposition of social network adjacency matrices. The use of semidiscrete decomposition in this application creates separability amongst observed data trails as a means to establish uniqueness amongst the users of the social network. The results presented herein illustrate to what extent a user's social fingerprint can be both quantified and identified on a social network through time. Copyright © 2015 John Wiley & Sons, Ltd.

We present an iterative approach to solve separable nonlinear least squares problems arising in the estimation of wavelength-dependent point spread function parameters for hyperspectral imaging. A variable projection Gauss–Newton method is used to solve the nonlinear least squares problem. An analysis shows that the Jacobian can be potentially very ill conditioned. To deal with this ill conditioning, we use a combination of subset selection and other regularization techniques. Experimental results related to hyperspectral point spread function parameter identification and star spectrum reconstruction illustrate the effectiveness of the resulting numerical scheme. Copyright © 2015 John Wiley & Sons, Ltd.

The null linear discriminant analysis method is a competitive approach for dimensionality reduction. The implementation of this method, however, is computationally expensive. Recently, a fast implementation of null linear discriminant analysis method was proposed in the paper of Sharma *et al.* In this method, a random matrix multiplication with scatter matrix is used to produce the orientation matrix. In this paper, we consider whether the random matrix can be replaced by an arbitrary full rank matrix of appropriate size, and focus on the necessary and sufficient condition to guarantee full column rank of the orientation matrix. We investigate how to choose the matrix properly, such that the two criteria of the null linear discriminant analysis method are satisfied. Furthermore, we give a necessary and sufficient condition to guarantee full column rank of the orientation matrix, and describe the geometric characterization of this condition. Numerical experiments justify our theoretical analysis. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we extend the inexact Uzawa algorithm of Hu and Zou to the nonsymmetric generalized saddle point problem (SPP). The techniques used here are similar with those in the paper of Bramble *et al.*, where the convergence of Uzawa type algorithm for solving nonsymmetric case depends on the spectrum of the preconditioners involved. The main contributions of this paper focus on two new linear Uzawa type algorithms for nonsymmetric generalized saddle point problems and their convergence properties. Our exact Uzawa algorithm can always converge without any prior estimate on the spectrum of two preconditioned subsystems involved. The convergence of our inexact Uzawa algorithm does not need the prior spectrum estimate of Schur complement, which may not be easy to achieve in practical applications. Numerical results of the algorithm on the Navier–Stokes problem are also presented. Copyright © 2015 John Wiley & Sons, Ltd.

We propose some computational methods for extracting static backgrounds from surveillance videos corrupted by noise, blur, or both. The new methods are constructed based on the fact that the matrix representation of a static background consists of identical columns; hence the idea of median filtering is embedded in these methods. These new methods significantly differ from existing methods originating from the robust principal component analysis (RPCA) in that no nuclear-norm term is involved; thus the computation of singular value decomposition can be completely avoided when solving these new models iteratively. This is an important feature because usually the dimensionality of a surveillance video is large and so the involved singular value decomposition (which is inevitable for RPCA-based models) is very expensive computationally. We show that these methods can be easily solved by well-developed operator splitting methods in optimization literature such as the alternating direction method of multipliers. We compare the new methods with their RPCA-based counterparts via testing some synthetic and real videos. Our numerical results show that compared with RPCA-based models, these median filtering-based vaiational models can extract more accurate backgrounds when the background in a surveillance video is static, and numerically, they can be solved much more efficiently. Copyright © 2015 John Wiley & Sons, Ltd.

Our method is a kind of exact interpolation scheme by Achi Brandt *et al*. In the exact interpolation scheme, for **x** being the fine-level approximation of the solution, the coarse-space *V* = *V*(**x**) is constructed so that **x** satisfies **x**∈*V*. We achieve it simply by adding the vector **x** as a first column of the prolongator. (The columns of the prolongator *P* form a computationally relevant basis of the coarse-space *V* = Range(*P*).) The advantages of this construction become obvious when solving non-linear problems. The cost of enriching the coarse-space *V* by the current approximation **x** is a single dense column of the prolongator that has to be updated each iteration. Our method can be used for multilevel acceleration of virtually any iterative method used for solving both linear and non-linear systems. Copyright © 2015 John Wiley & Sons, Ltd.

This paper is contributed to a fast algorithm for Hankel tensor–vector products. First, we explain the necessity of fast algorithms for Hankel and block Hankel tensor–vector products by sketching the algorithm for both one-dimensional and multi-dimensional exponential data fitting. For proposing the fast algorithm, we define and investigate a special class of Hankel tensors that can be diagonalized by the Fourier matrices, which is called anti-circulant tensors. Then, we obtain a fast algorithm for Hankel tensor–vector products by embedding a Hankel tensor into a larger anti-circulant tensor. The computational complexity is about for a square Hankel tensor of order *m* and dimension *n*, and the numerical examples also show the efficiency of this scheme. Moreover, the block version for multi-level block Hankel tensors is discussed. Copyright © 2015 John Wiley & Sons, Ltd.

A fast approximate inversion method is proposed for the block lower triangular Toeplitz with tri-diagonal blocks (BL3TB) matrix. The BL3TB matrix is approximated by a block *ϵ*-circulant matrix, which can be efficiently inverted using the fast Fourier transforms. The error estimation is given to show the high accuracy of the approximation. In applications, the proposed method is employed to solve the fractional sub-diffusion equation whose discretized matrix by a finite difference method is a BL3TB matrix. Numerical experiments are carried out to demonstrate the efficiency of the proposed method. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, the idea of auxiliary space multigrid methods is introduced. The construction is based on a two-level block factorization of local (finite element stiffness) matrices associated with a partitioning of the domain into overlapping or non-overlapping subdomains. The two-level method utilizes a coarse-grid operator obtained from additive Schur complement approximation. Its analysis is carried out in the framework of auxiliary space preconditioning and condition number estimates for both the two-level preconditioner and the additive Schur complement approximation are derived. The two-level method is recursively extended to define the auxiliary space multigrid algorithm. In particular, so-called Krylov cycles are considered. The theoretical results are supported by a representative collection of numerical tests that further demonstrate the efficiency of the new algorithm for multiscale problems. Copyright © 2014 John Wiley & Sons, Ltd.

The following paper discusses the application of a multigrid-in-time scheme to Least Squares Shadowing (LSS), a novel sensitivity analysis method for chaotic dynamical systems. While traditional sensitivity analysis methods break down for chaotic dynamical systems, LSS is able to compute accurate gradients. Multigrid is used because LSS requires solving a very large Karush–Kuhn–Tucker system constructed from the solution of the dynamical system over the entire time interval of interest. Several different multigrid-in-time schemes are examined, and a number of factors were found to heavily influence the convergence rate of multigrid-in-time for LSS. These include the iterative method used for the smoother, how the coarse grid system is formed and how the least squares objective function at the center of LSS is weighted. Copyright © 2014 John Wiley & Sons, Ltd.

No abstract is available for this article.

]]>The main goal of this paper is to approximate the principal *p*th root of a matrix by using a family of high-order iterative methods. We analyse the semi-local convergence and the speed of convergence of these methods. Concerning stability, it is well known that even the simplified Newton method is unstable. Despite it, we present stable versions of our family of algorithms. We test numerically the methods: we check the numerical robustness and stability by considering matrices that are close to be singular and badly conditioned. We find algorithms of the family with better numerical behavior than the Newton and the Halley methods. These two algorithms are basically the iterative methods proposed in the literature to solve this problem. Copyright © 2015 John Wiley & Sons, Ltd.

This article is dedicated to the rapid computation of separable expansions for the approximation of random fields. We consider approaches based on techniques from the approximation of non-local operators on the one hand and based on the pivoted Cholesky decomposition on the other hand. Especially, we provide an *a posteriori* error estimate for the pivoted Cholesky decomposition in terms of the trace. Numerical examples are provided to validate and quantify the presented methods. Copyright © 2015 John Wiley & Sons, Ltd.

Multigrid and related multilevel methods are the approaches of choice for solving linear systems that result from discretization of a wide class of PDEs. A large gap, however, exists between the theoretical analysis of these algorithms and their actual performance. This paper focuses on the extension of the well-known local mode (often local Fourier) analysis approach to a wider class of problems. The semi-algebraic mode analysis (SAMA) proposed here couples standard local Fourier analysis approaches with algebraic computation to enable analysis of a wider class of problems, including those with strong advective character. The predictive nature of SAMA is demonstrated by applying it to the parabolic diffusion equation in one and two space dimensions, elliptic diffusion in layered media, as well as a two-dimensional convection-diffusion problem. These examples show that accounting for boundary conditions and heterogeneity enables accurate predictions of the short-term and asymptotic convergence behavior for multigrid and related multilevel methods. Copyright © 2015 John Wiley & Sons, Ltd.

The use of matchings is a powerful technique for scaling and ordering sparse matrices prior to the solution of a linear system *A**x* = *b*. Traditional methods such as implemented by the HSL software package MC64 use the Hungarian algorithm to solve the maximum weight maximum cardinality matching problem. However, with advances in the algorithms and hardware used by direct methods for the parallelization of the factorization and solve phases, the serial Hungarian algorithm can represent an unacceptably large proportion of the total solution time for such solvers. Recently, auction algorithms and approximation algorithms have been suggested as alternatives for achieving near-optimal solutions for the maximum weight maximum cardinality matching problem. In this paper, the efficacy of auction and approximation algorithms as replacements for the Hungarian algorithm is assessed in the context of sparse symmetric direct solvers when used in problems arising from a range of practical applications. High-cardinality suboptimal matchings are shown to be as effective as optimal matchings for the purposes of scaling. However, matching-based ordering techniques require that matchings are much closer to optimality before they become effective. The auction algorithm is demonstrated to be capable of finding such matchings significantly faster than the Hungarian algorithm, but our
-approximation matching approach fails to consistently achieve a sufficient cardinality. Copyright © 2015 The Authors Numerical Linear Algebra with Applications Published by John Wiley & Sons Ltd.

We present a comparison of different multigrid approaches for the solution of systems arising from high-order continuous finite element discretizations of elliptic partial differential equations on complex geometries. We consider the pointwise Jacobi, the Chebyshev-accelerated Jacobi, and the symmetric successive over-relaxation smoothers, as well as elementwise block Jacobi smoothing. Three approaches for the multigrid hierarchy are compared: (1) high-order *h*-multigrid, which uses high-order interpolation and restriction between geometrically coarsened meshes; (2) *p*-multigrid, in which the polynomial order is reduced while the mesh remains unchanged, and the interpolation and restriction incorporate the different-order basis functions; and (3) a first-order approximation multigrid preconditioner constructed using the nodes of the high-order discretization. This latter approach is often combined with algebraic multigrid for the low-order operator and is attractive for high-order discretizations on unstructured meshes, where geometric coarsening is difficult. Based on a simple performance model, we compare the computational cost of the different approaches. Using scalar test problems in two and three dimensions with constant and varying coefficients, we compare the performance of the different multigrid approaches for polynomial orders up to 16. Overall, both *h*-multigrid and *p*-multigrid work well; the first-order approximation is less efficient. For constant coefficients, all smoothers work well. For variable coefficients, Chebyshev and symmetric successive over-relaxation smoothing outperform Jacobi smoothing. While all of the tested methods converge in a mesh-independent number of iterations, none of them behaves completely independent of the polynomial order. When multigrid is used as a preconditioner in a Krylov method, the iteration number decreases significantly compared with using multigrid as a solver. Copyright © 2015 John Wiley & Sons, Ltd.

A two-grid convergence analysis based on the paper [*Algebraic analysis of aggregation-based multigrid*, by A. Napov and Y. Notay, Numer. Lin. Alg. Appl. 18 (2011), pp. 539–564] is derived for various aggregation schemes applied to a finite element discretization of a rotated anisotropic diffusion equation. As expected, it is shown that the best aggregation scheme is one in which aggregates are aligned with the anisotropy. In practice, however, this is not what automatic aggregation procedures do. We suggest approaches for determining appropriate aggregates based on eigenvectors associated with small eigenvalues of a block splitting matrix or based on minimizing a quantity related to the spectral radius of the iteration matrix. Copyright © 2015 John Wiley & Sons, Ltd.

We analyze inexact fixed-point iterations where the generating function contains an inexact solve of an equation system to answer the question of how tolerances for the inner solves influence the iteration error of the outer fixed-point iteration. Important applications are the Picard iteration and partitioned fluid-structure interaction. For the analysis, the iteration is modeled as a perturbed fixed-point iteration, and existing analysis is extended to the nested case **x** = **F**(**S**(**x**)). We prove that if the iteration converges, it converges to the exact solution irrespective of the tolerance in the inner systems, provided that a nonstandard relative termination criterion is employed, whereas standard relative and absolute criteria do not have this property. Numerical results demonstrate the effectiveness of the approach with the nonstandard termination criterion. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we consider the solution of linear systems of saddle point type by correcting the Uzawa algorithm, which has been proposed in [K. Arrow, L. Hurwicz, H. Uzawa, Studies in nonlinear programming, Stanford University Press, Stanford, CA, 1958]. We call this method as corrected Uzawa (CU) method. The convergence of the CU method is analyzed for solving nonsingular saddle point problem as well as the semi-convergence for the singular case. First, the corrected model for the Uzawa algorithm is established, and the CU algorithm is presented. Then we study the geometric meaning of the CU model. Moreover, we introduce the overall reduction coefficient *α* to measure the effect of the CU process. It is shown that the CU method converges faster than the Uzawa method and several other methods if the overall reduction coefficient *α* satisfies certain conditions. Numerical experiments are presented to illustrate the theoretical results and examine the numerical effectiveness of the CU method. Copyright © 2015 John Wiley & Sons, Ltd.

Symmetric collocation methods with RBFs allow approximation of the solution of a partial differential equation, even if the right-hand side is only known at scattered data points, without needing to generate a grid. However, the benefit of a guaranteed symmetric positive definite block system comes at a high computational cost. This cost can be alleviated somewhat by considering compactly supported RBFs and a multiscale technique. But the condition number and sparsity will still deteriorate with the number of data points. Therefore, we study certain block diagonal and triangular preconditioners. We investigate ideal preconditioners and determine the spectra of the preconditioned matrices before proposing more practical preconditioners based on a restricted additive Schwarz method with coarse grid correction. Numerical results verify the effectiveness of the preconditioners. Copyright © 2015 John Wiley & Sons, Ltd.

For solving the large sparse linear complementarity problems, we establish modified modulus-based matrix splitting iteration methods and present the convergence analysis when the system matrices are *H*_{+}-matrices. The optima of parameters involved under some scopes are also analyzed. Numerical results show that in computing efficiency, our new methods are superior to classical modulus-based matrix splitting iteration methods under suitable conditions. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we present a preconditioned variant of the generalized successive overrelaxation (GSOR) iterative method for solving a broad class of complex symmetric linear systems. We study conditions under which the spectral radius of the iteration matrix of the preconditioned GSOR method is smaller than that of the GSOR method and determine the optimal values of iteration parameters. Numerical experiments are given to verify the validity of the presented theoretical results and the effectiveness of the preconditioned GSOR method. Copyright © 2015 John Wiley & Sons, Ltd.

We revisit the shift-and-invert Arnoldi method proposed in [S. Lee, H. Pang, and H. Sun. *Shift-invert Arnoldi approximation to the Toeplitz matrix exponential*, SIAM J. Sci. Comput., 32: 774–792, 2010] for numerical approximation to the product of Toeplitz matrix exponential with a vector. In this approach, one has to solve two large-scale Toeplitz linear systems in advance. However, if the desired accuracy is high, the cost will be prohibitive. Therefore, it is interesting to investigate how to solve the Toeplitz systems inexactly in this method. The contribution of this paper is in three regards. First, we give a new stability analysis on the Gohberg–Semencul formula (GSF) and define the GSF condition number of a Toeplitz matrix. It is shown that when the size of the Toeplitz matrix is large, our result is sharper than the one given in [M. Gutknecht and M. Hochbruck. *The stability of inversion formulas for Toeplitz matrices*, Linear Algebra Appl., 223/224: 307–324, 1995]. Second, we establish a relation between the error of Toeplitz systems and the residual of Toeplitz matrix exponential. We show that if the GSF condition number of the Toeplitz matrix is medium-sized, then the Toeplitz systems can be solved in a low accuracy. Third, based on this relationship, we present a practical stopping criterion for relaxing the accuracy of the Toeplitz systems and propose an inexact shift-and-invert Arnoldi algorithm for the Toeplitz matrix exponential problem. Numerical experiments illustrate the numerical behavior of the new algorithm and show the effectiveness of our theoretical results. Copyright © 2015 John Wiley & Sons, Ltd.