Matrix exponential and Krylov subspaces for fast time domain computations: recent advances

We show how finite difference (or finite element) time domain computations can be accelerated by employing recent advances in the matrix exponential time integration and Krylov subspace techniques

A block Krylov subspace time-exact solution method for linear ODE systems

We propose a time-exact Krylov-subspace-based method for solving linear ODE (ordinary differential equation) systems of the form $y'=-Ay + g(t)$ and $y''=-Ay + g(t)$, where $y(t)$ is the unknown function. The method consists of two stages. The first stage is an accurate piecewise polynomial approximation of the source term $g(t)$, constructed with the help of the truncated SVD (singular value decomposition). The second stage is a special residual-based block Krylov subspace method. The accuracy of the method is only restricted by the accuracy of the piecewise polynomial approximation and by the error of the block Krylov process. Since both errors can, in principle, be made arbitrarily small, this yields, at some costs, a time-exact method. Numerical experiments are presented to demonstrate efficiency of the new method, as compared to an exponential time integrator with Krylov subspace matrix function evaluations

Residual, restarting and Richardson iteration for the matrix exponential

A well-known problem in computing some matrix functions iteratively is a lack of a clear, commonly accepted residual notion. An important matrix function for which this is the case is the matrix exponential. Assume, the matrix exponential of a given matrix times a given vector has to be computed. We interpret the sought after vector as a value of a vector function satisfying the linear system of ordinary differential equations (ODE), whose coefficients form the given matrix. The residual is then defined with respect to the initial-value problem for this ODE system. The residual introduced in this way can be seen as a backward error. We show how the residual can efficiently be computed within several iterative methods for the matrix exponential. This completely resolves the question of reliable stopping criteria for these methods. Furthermore, we show that the residual concept can be used to construct new residual-based iterative methods. In particular, a variant of the Richardson method for the new residual appears to provide an efficient way to restart Krylov subspace methods for evaluating the matrix exponential.\u

Improving approximate matrix factorizations for implicit time integration in air pollution modelling

For a long time operator splitting was the only computationally feasible way of implicit time integration in large scale Air Pollution Models. A recently proposed attractive alternative is Rosenbrock schemes combined with Approximate Matrix Factorization (AMF). With AMF, linear systems arising in implicit time stepping are solved approximately in such a way that the overall computational costs per time step are not higher than those of splitting methods. We propose and discuss two new variants of AMF. The first one is aimed at yet a further reduction of costs as compared with conventional AMF. The second variant of AMF provides in certain circumstances a better approximation to the inverse of the linear system matrix than standard AMF and requires the same computational work

Experiments with MRAI time stepping schemes on a distributed memory parallel environment

Implicit time stepping is often difficult to parallelize. The recently proposed Minimal Residual Approximate Implicit (MRAI) schemes are specially designed as a cheaper and parallelizable alternative for implicit time stepping. A several GMRES iterations are performed to solve approximately the implicit scheme of interest, and the step size is adjusted to guarantee stability. A natural way to apply the approach is to modify a given implicit scheme in which one is interested. Here, we present numerical results for two parallel implementations of MRAI schemes. One is based on the simple Euler Backward scheme, and the other is the MRAI-modified multistep ODE solver LSODE. On the Cray T3E and IBM SP2 platforms, the MRAI codes exhibit parallelism of explicit schemes. The model problem under consideration is the 3D spatially discretized heat equation. The speed-up results for the Cray T3E and IBM SP2 are reported

On Skew-symmetric Preconditioning for Strongly Non-symmetric Linear Systems

To solve iteratively linear system $Au=b$ with large sparse strongly non-symmetric matrix $A$ we propose preconditioning $\hat A \hat u = \hat b$, $\hat A=(I+\tau L_1)^{-1} A (I+\tau U_1)^{-1},\; \tau>0$ where respectively lower and upper triangular matrices $L_1$ and $U_1$ are so that $L_1+U_1=1/2(A-A^*)$. Such preconditioning technique may be treated as a variant of ILU-factorization, and we call it MSSILU --- Modified Skew-Symmetric ILU. \ud We investigate and optimize (with respect to $\tau$) convergence of preconditioned Richardson method (RM) of the following special form: ${\hat x}^{m+1}=(I-\tau \hat A){\hat x}^m+\tau {\hat b},\; m\geq 0$, where $\tau$ is the same as in $\hat A$. For this method we give an estimate for rate of convergence in relevant Euclidean norm for the case of positivereal matrix $A$. \ud Numerical experiments have included solving linear systems arising from 5-point FD approximation of convection--diffusion equation with dominated convection by MSSILU+RM, MSSILU+GMRES(2) and MSSILU+GMRES(10).\u

Parallel processing and non-uniform grids in global air quality modeling

A large-scale global air quality model, running efficiently on a single vector processor, is enhanced to make more realistic and more long-term simulations feasible. Two strategies are combined: non-uniform grids and parallel processing. The communication through the hierarchy of non-uniform grids interferes with the inter-processor communication. We discuss load balance in the decomposition of the domain, I/O, and inter-processor communication. A model shows that the communication overhead for both techniques is very low, whence non-uniform grids allow for large speed-ups and high speed-up can be expected from parallelization. The implementation is in progress, and results of experiments will be reported elsewhere