Thesis (MSc)--Stellenbosch University, 2011. / ENGLISH ABSTRACT: In several numerical approaches to PDEs shifted linear systems of the form (zI - A)x = b, need to be solved for several values of the complex scalar z. Often, these linear systems are large and sparse. This thesis investigates efficient numerical methods for these systems that arise from a contour integral approximation to PDEs and compares these methods with direct solvers. In the first part, we present three model PDEs and discuss numerical approaches to solve them. We use the first problem to demonstrate computations with a dense matrix, the second problem to demonstrate computations with a sparse symmetric matrix and the third problem for a sparse but nonsymmetric matrix. To solve the model PDEs numerically we apply two space discrerization methods, namely the finite difference method and the Chebyshev collocation method. The contour integral method mentioned above is used to integrate with respect to the time variable. In the second part, we study a Hessenberg reduction method for solving shifted linear systems with a dense matrix and present numerical comparison of it with the built-in direct linear system solver in SciPy. Since both are direct methods, in the absence of roundoff errors, they give the same result. However, we find that the Hessenberg reduction method is more efficient in CPU-time than the direct solver. As application we solve a one-dimensional version of the heat equation. In the third part, we present efficient techniques for solving shifted systems with a sparse matrix by Krylov subspace methods. Because of their shift-invariance property, the Krylov methods allow one to obtain approximate solutions for all values of the parameter, by generating a single approximation space. Krylov methods applied to the linear systems are generally slowly convergent and hence preconditioning is necessary to improve the convergence. The use of shift-invert preconditioning is discussed and numerical comparisons with a direct sparse solver are presented. As an application we solve a two-dimensional version of the heat equation with and without a convection term. Our numerical experiments show that the preconditioned Krylov methods are efficient in both computational time and memory space as compared to the direct sparse solver. / AFRIKAANSE OPSOMMING: In verskeie numeriese metodes vir PDVs moet geskuifde lineêre stelsels van die vorm (zI − A)x = b, opgelos word vir verskeie waardes van die komplekse skalaar z. Hierdie stelsels is dikwels groot en yl. Hierdie tesis ondersoek numeriese metodes vir sulke stelsels wat voorkom in kontoerintegraalbenaderings vir PDVs en vergelyk hierdie metodes met direkte metodes vir oplossing. In die eerste gedeelte beskou ons drie model PDVs en bespreek numeriese benaderings om hulle op te los. Die eerste probleem word gebruik om berekenings met ’n vol matriks te demonstreer, die tweede probleem word gebruik om berekenings met yl, simmetriese matrikse te demonstreer en die derde probleem vir yl, onsimmetriese matrikse. Om die model PDVs numeries op te los beskou ons twee ruimte-diskretisasie metodes, naamlik die eindige-verskilmetode en die Chebyshev kollokasie-metode. Die kontoerintegraalmetode waarna hierbo verwys is word gebruik om met betrekking tot die tydveranderlike te integreer. In die tweede gedeelte bestudeer ons ’n Hessenberg ontbindingsmetode om geskuifde lineêre stelsels met ’n vol matriks op te los, en ons rapporteer numeriese vergelykings daarvan met die ingeboude direkte oplosser vir lineêre stelsels in SciPy. Aangesien beide metodes direk is lewer hulle dieselfde resultate in die afwesigheid van afrondingsfoute. Ons het egter bevind dat die Hessenberg ontbindingsmetode meer effektief is in terme van rekenaartyd in vergelyking met die direkte oplosser. As toepassing los ons ’n een-dimensionele weergawe van die hittevergelyking op. In die derde gedeelte beskou ons effektiewe tegnieke om geskuifde stelsels met ’n yl matriks op te los, met Krylov subruimte-metodes. As gevolg van hul skuifinvariansie eienskap, laat die Krylov metodes mens toe om benaderde oplossings te verkry vir alle waardes van die parameter, deur slegs een benaderingsruimte voort te bring. Krylov metodes toegepas op lineêre stelsels is in die algemeen stadig konvergerend, en gevolglik is prekondisionering nodig om die konvergensie te verbeter. Die gebruik van prekondisionering gebasseer op skuif-en-omkeer word bespreek en numeriese vergelykings met direkte oplossers word aangebied. As toepassing los ons ’n twee-dimensionele weergawe van die hittevergelyking op, met ’n konveksie term en daarsonder. Ons numeriese eksperimente dui aan dat die Krylov metodes met prekondisionering effektief is, beide in terme van berekeningstyd en rekenaargeheue, in vergelyking met die direkte metodes.

In this thesis we study methods for solving systems of linear algebraic equati- ons with multiple right hand sides. Specifically we focus on block Krylov subspace methods and global Krylov subspace methods, which can be derived by various approaches to generalization of methods GMRES and LSQR for solving systems of linear equations with single right hand side. We describe the difference in construction of orthonormal basis in block methods and F-orthonormal basis in global methods, in detail. Finally, we provide numerical experiments for the deri- ved algorithms in MATLAB enviroment. On carefully selected test problems we compare convergence properties of the methods. 1

The SHIRA method of Mehrmann and Watkins belongs among the structure preserving Krylov subspace methods for solving skew-Hamiltonian eigenvalue problems. It can also be applied to Hamiltonian eigenproblems by considering a suitable transformation. Structure induced shift-and-invert techniques are employed to steer the algorithm towards the interesting region of the spectrum. However, the shift cannot be altered in the middle of the computation without discarding the information that has been accumulated so far. This paper shows how SHIRA can be combined with ideas from Ruhe's Rational Krylov algorithm to yield a method that permits an adjustment of shift after every step of the computation, adding greatly to the flexibility of the algorithm. We call this new method rational SHIRA. A numerical example is presented to demonstrate its efficiency.

Partial differential equations are the chief means of providing mathematical models in science, engineering and other fields. Optimal control of partial differential equations (PDEs) has tremendous applications in engineering and science, such as shape optimization, image processing, fluid dynamics, and chemical processes. In this thesis, we develop and analyze several efficient numerical methods for the optimal control problems governed by elliptic PDE, parabolic PDE, and wave PDE, respectively. The thesis consists of six chapters. In Chapter 1, we briefly introduce a few motivating applications and summarize some theoretical and computational foundations of our following developed approaches. In Chapter 2, we establish a new multigrid algorithm to accelerate the semi-smooth Newton method that is applied to the first-order necessary optimality system arising from semi-linear control-constrained elliptic optimal control problems. Under suitable assumptions, the discretized Jacobian matrix is proved to have a uniformly bounded inverse with respect to mesh size. Different from current available approaches, a new strategy that leads to a robust multigrid solver is employed to define the coarse grid operator. Numerical simulations are provided to illustrate the efficiency of the proposed method, which shows to be computationally more efficient than the popular full approximation storage (FAS) multigrid method. In particular, our proposed approach achieves a mesh-independent convergence and its performance is highly robust with respect to the regularization parameter. In Chaper 3, we present a new second-order leapfrog finite difference scheme in time for solving the first-order necessary optimality system of the linear parabolic optimal control problems. The new leapfrog scheme is shown to be unconditionally stable and it provides a second-order accuracy, while the classical leapfrog scheme usually is well-known to be unstable. A mathematical proof for the convergence of the proposed scheme is provided under a suitable norm. Moreover, the proposed leapfrog scheme gives a favorable structure that leads to an effective implementation of a fast solver under the multigrid framework. Numerical examples show that the proposed scheme significantly outperforms the widely used second-order backward time differentiation approach, and the resultant fast solver demonstrates a mesh-independent convergence as well as a linear time complexity. In Chapter 4, we develop a new semi-smooth Newton multigrid algorithm for solving the discretized first-order necessary optimality system that characterizes the optimal solution of semi-linear parabolic PDE optimal control problems with control constraints. A new leapfrog discretization scheme in time associated with the standard five-point stencil in space is established to achieve a second-order accuracy. The convergence (or unconditional stability) of the proposed scheme is proved when time-periodic solutions are considered. Moreover, the derived well-structured discretized Jacobian matrices greatly facilitate the development of an effective smoother in our multigrid algorithm. Numerical simulations are provided to illustrate the effectiveness of the proposed method, which validates the second-order accuracy in solution approximations as well as the optimal linear complexity of computing time. In Chapter 5, we offer a new implicit finite difference scheme in time for solving the first-order necessary optimality system arising in optimal control of wave equations. With a five-point central finite difference scheme in space, the full discretization is proved to be unconditionally convergent with a second-order accuracy, which is not restricted by the classical Courant-Friedrichs-Lewy (CFL) stability condition on the spatial and temporal step sizes. Moreover, based on its advantageous developed structure, an efficient preconditioned Krylov subspace method is provided and analyzed for solving the discretized sparse linear system. Numerical examples are presented to confirm our theoretical conclusions and demonstrate the promising performance of proposed preconditioned iterative solver. Finally, brief summaries and future research perspectives are given in Chapter 6.

No description available.

The purpose of this thesis is to study how to overcome difficulties one typically encounters when solving non-negative inverse problems by standard Krylov subspace methods. We first give a theoretical background to the non-negative inverse problems. Then we concentrate on selected modifications of Krylov subspace methods known to improve the solution significantly. We describe their properties, provide their implementation and propose an improvement for one of them. After that, numerical experiments are presented giving a comparison of the methods and analyzing the influence of the present parameters on the behavior of the solvers. It is clearly demonstrated, that the methods imposing nonnegativity perform better than the unconstrained methods. Moreover, our improvement leads in some cases to a certain reduction of the number of iterations and consequently to savings of the computational time while preserving a good quality of the approximation.

We propose a new framework for understanding block Krylov subspace methods, which hinges on a matrix-valued inner product. We can recast the ``classical" block Krylov methods, such as O'Leary's block conjugate gradients, global methods, and loop-interchange methods, within this framework. Leveraging the generality of the framework, we develop an efficient restart procedure and error bounds for the shifted block full orthogonalization method (Sh-BFOM(m)). Regarding BFOM as the prototypical block Krylov subspace method, we propose another formalism, which we call modified BFOM, and show that block GMRES and the new block Radau-Lanczos method can be regarded as modified BFOM. In analogy to Sh-BFOM(m), we develop an efficient restart procedure for shifted BGMRES with restarts (Sh-BGMRES(m)), as well as error bounds. Using this framework and shifted block Krylov methods with restarts as a foundation, we formulate block Krylov subspace methods with restarts for matrix functions acting on multiple vectors f(A)B. We obtain convergence bounds for \bfomfom (BFOM for Functions Of Matrices) and block harmonic methods (i.e., BGMRES-like methods) for matrix functions. With various numerical examples, we illustrate our theoretical results on Sh-BFOM and Sh-BGMRES. We also analyze the matrix polynomials associated to the residuals of these methods. Through a variety of real-life applications, we demonstrate the robustness and versatility of B(FOM)^2 and block harmonic methods for matrix functions. A particularly interesting example is the tensor t-function, our proposed definition for the function of a tensor in the tensor t-product formalism. Despite the lack of convergence theory, we also show that the block Radau-Lanczos modification can reduce the number of cycles required to converge for both linear systems and matrix functions. / Mathematics

Science and engineering problems frequently require solving a sequence of single linear systems or a sequence of dual linear systems. We develop algorithms that recycle Krylov subspaces and preconditioners from one system (or pair of systems) in the sequence to the next, leading to efficient solutions. Besides the benefit of only having to store few Lanczos vectors, using BiConjugate Gradients (BiCG) to solve dual linear systems may have application-specific advantages. For example, using BiCG to solve the dual linear systems arising in interpolatory model reduction provides a backward error formulation in the model reduction framework. Using BiCG to evaluate bilinear forms -- for example, in the variational Monte Carlo (VMC) algorithm for electronic structure calculations -- leads to a quadratic error bound. Since one of our focus areas is sequences of dual linear systems, we introduce recycling BiCG, a BiCG method that recycles two Krylov subspaces from one pair of dual linear systems to the next pair. The derivation of recycling BiCG also builds the foundation for developing recycling variants of other bi-Lanczos based methods like CGS, BiCGSTAB, BiCGSTAB2, BiCGSTAB(l), QMR, and TFQMR. We develop a generalized bi-Lanczos algorithm, where the two matrices of the bi-Lanczos procedure are not each other's conjugate transpose but satisfy this relation over the generated Krylov subspaces. This is sufficient for a short term recurrence. Next, we derive an augmented bi-Lanczos algorithm with recycling and show that this algorithm is a special case of generalized bi-Lanczos. The Petrov-Galerkin approximation that includes recycling in the iteration leads to modified two-term recurrences for the solution and residual updates. We generalize and extend the framework of our recycling BiCG to CGS, BiCGSTAB and BiCGSTAB2. We perform extensive numerical experiments and analyze the generated recycle space. We test all of our recycling algorithms on a discretized partial differential equation (PDE) of convection-diffusion type. This PDE problem provides well-known test cases that are easy to analyze further. We use recycling BiCG in the Iterative Rational Krylov Algorithm (IRKA) for interpolatory model reduction and in the VMC algorithm. For a model reduction problem, we show up to 70% savings in iterations, and we also demonstrate that solving the problem without recycling leads to (about) a 50% increase in runtime. Experiments with recycling BiCG for VMC gives promising results. We also present an algorithm that recycles preconditioners, leading to a dramatic reduction in the cost of VMC for large(r) systems. The main cost of the VMC method is in constructing a sequence of Slater matrices and computing the ratios of determinants for successive Slater matrices. Recent work has improved the scaling of constructing Slater matrices for insulators, so that the cost of constructing Slater matrices in these systems is now linear in the number of particles. However, the cost of computing determinant ratios remains cubic in the number of particles. With the long term aim of simulating much larger systems, we improve the scaling of computing determinant ratios in the VMC method for simulating insulators by using preconditioned iterative solvers. The main contribution here is the development of a method to efficiently compute for the Slater matrices a sequence of preconditioners that make the iterative solver converge rapidly. This involves cheap preconditioner updates, an effective reordering strategy, and a cheap method to monitor instability of ILUTP preconditioners. Using the resulting preconditioned iterative solvers to compute determinant ratios of consecutive Slater matrices reduces the scaling of the VMC algorithm from O(n^3) per sweep to roughly O(n^2), where n is the number of particles, and a sweep is a sequence of n steps, each attempting to move a distinct particle. We demonstrate experimentally that we can achieve the improved scaling without increasing statistical errors. / Ph. D.

Les travaux de ce doctorat concernent le développement de méthodes itératives pour la résolution de systèmes linéaires creux de grande taille comportant de nombreux seconds membres. L’application visée est la résolution d’un problème inverse en géophysique visant à reconstruire la vitesse de propagation des ondes dans le sous-sol terrestre. Lorsque de nombreuses sources émettrices sont utilisées, ce problème inverse nécessite la résolution de systèmes linéaires complexes non symétriques non hermitiens comportant des milliers de seconds membres. Dans le cas tridimensionnel ces systèmes linéaires sont reconnus comme difficiles à résoudre plus particulièrement lorsque des fréquences élevées sont considérées. Le principal objectif de cette thèse est donc d’étendre les développements existants concernant les méthodes de Krylov par bloc. Nous étudions plus particulièrement les techniques de déflation dans le cas multiples seconds membres et recyclage de sous-espace dans le cas simple second membre. Des gains substantiels sont obtenus en terme de temps de calcul par rapport aux méthodes existantes sur des applications réalistes dans un environnement parallèle distribué. / This PhD thesis concerns the development of flexible Krylov subspace iterative solvers for the solution of large sparse linear systems of equations with multiple right-hand sides. Our target application is the solution of the acoustic full waveform inversion problem in geophysics associated with the phenomena of wave propagation through an heterogeneous model simulating the subsurface of Earth. When multiple wave sources are being used, this problem gives raise to large sparse complex non-Hermitian and nonsymmetric linear systems with thousands of right-hand sides. Specially in the three-dimensional case and at high frequencies, this problem is known to be difficult. The purpose of this thesis is to develop a flexible block Krylov iterative method which extends and improves techniques already available in the current literature to the multiple right-hand sides scenario. We exploit the relations between each right-hand side to accelerate the convergence of the overall iterative method. We study both block deflation and single right-hand side subspace recycling techniques obtaining substantial gains in terms of computational time when compared to other strategies published in the literature, on realistic applications performed in a parallel environment.

La solution numérique des systèmes dynamiques est un moyen efficace pour étudier des phénomènes physiques complexes. Cependant, dans un cadre à grande échelle, la dimension du système rend les calculs infaisables en raison des limites de mémoire et de temps, ainsi que le mauvais conditionnement. La solution de ce problème est la réduction de modèles. Cette thèse porte sur les méthodes de projection pour construire efficacement des modèles d'ordre inférieur à partir des systèmes linéaires dynamiques de grande taille. En particulier, nous nous intéressons à la projection sur la réunion de plusieurs sous-espaces de Krylov standard qui conduit à une classe de modèles d'ordre réduit. Cette méthode est connue par l'interpolation rationnelle. En se basant sur ce cadre théorique qui relie la projection de Krylov à l'interpolation rationnelle, quatre algorithmes de type Lanczos rationnel pour la réduction de modèles sont proposés. Dans un premier temps, nous avons introduit une méthode adaptative de type Lanczos rationnel par block pour réduire l'ordre des systèmes linéaires dynamiques de grande taille, cette méthode est basée sur l'algorithme de Lanczos rationnel par block et une méthode adaptative pour choisir les points d'interpolation. Une généralisation de ce premier algorithme est également donnée, où différentes multiplicités sont considérées pour chaque point d'interpolation. Ensuite, nous avons proposé une autre extension de la méthode du sous-espace de Krylov standard pour les systèmes à plusieurs-entrées plusieurs-sorties, qui est le sous-espace de Krylov global. Nous avons obtenu des équations qui décrivent cette procédure. Finalement, nous avons proposé une méthode de Lanczos étendu par block et nous avons établi de nouvelles propriétés algébriques pour cet algorithme. L'efficacité et la précision de tous les algorithmes proposés, appliqués sur des problèmes de réduction de modèles, sont testées dans plusieurs exemples numériques. / Numerical solution of dynamical systems have been a successful means for studying complex physical phenomena. However, in large-scale setting, the system dimension makes the computations infeasible due to memory and time limitations, and ill-conditioning. The remedy of this problem is model reductions. This dissertations focuses on projection methods to efficiently construct reduced order models for large linear dynamical systems. Especially, we are interesting by projection onto unions of Krylov subspaces which lead to a class of reduced order models known as rational interpolation. Based on this theoretical framework that relate Krylov projection to rational interpolation, four rational Lanczos-type algorithms for model reduction are proposed. At first, an adaptative rational block Lanczos-type method for reducing the order of large scale dynamical systems is introduced, based on a rational block Lanczos algorithm and an adaptive approach for choosing the interpolation points. A generalization of the first algorithm is also given where different multiplicities are consider for each interpolation point. Next, we proposed another extension of the standard Krylov subspace method for Multiple-Input Multiple-Output (MIMO) systems, which is the global Krylov subspace, and we obtained also some equations that describe this process. Finally, an extended block Lanczos method is introduced and new algebraic properties for this algorithm are also given. The accuracy and the efficiency of all proposed algorithms when applied to model order reduction problem are tested by means of different numerical experiments that use a collection of well known benchmark examples.