Spectral approximation with matrices issued from discretized operators

Keywords: AMS classification -2010 Mathematics Subject Classification, (MSC[2010]): 33 Special functions, 33F Computational aspects, 33F05 Numerical approximation and evaluation, 45 Integral Equations, 45C Eigenvalue problems, 45C05 Eigenvalue problems, 45E Singular integral equations, 65 Numerical analysis, 65F Numerical linear algebra, 65F15 Eigenvalues, eigenvectors, 65Y Computer aspects of numerical algorithms, 65Y20 Complexity and performance of numerical algorithms avec classification AMS -2010 Mathematics Subject Classification, (MSC[2010]): 33 Special functions, 33F Computational aspects, 33F05 Numerical approximation and evaluation, 45 Integral Equations, 45C Eigenvalue problems, 45C05 Eigenvalue problems, 45E Singular integral equations, 65 Numerical analysis, 65F Numerical linear algebra, 65F15 Eigenvalues, eigenvectors, 65Y Computer aspects of numerical algorithms, 65Y20 Complexity

In this thesis, we consider the numerical solution of a large eigenvalue problem in which the integral operator comes from a radiative transfer problem.

It is considered the use of hierarchical matrices, an efficient data-sparse representation of matrices, especially useful for large dimensional problems.

It consists on low-rank subblocks leading to low memory requirements as well as cheap computational costs.

We discuss the use of the hierarchical matrix technique in the numerical solution of a large scale eigenvalue problem arising from a finite rank discretization of an integral operator. The operator is of convolution type, it is defined through the first exponential-integral function and hence it is weakly singular.

We access HLIB (Hierarchical matrices LIBrary) that provides, among others, routines for the construction of hierarchical matrix structures and arithmetic algorithms to perform approximative matrix operations. Moreover, it is incorporated the matrix-vector multiply routines from HLIB, as well as LU factorization for preconditioning, into SLEPc (Scalable Library for Eigenvalue Problem Computations) in order to exploit the available algorithms to solve eigenvalue problems. vi It is also developed analytical expressions for the approximate degenerate kernels and deducted error upper bounds for these approximations.

The numerical results obtained with other approaches to solve the problem are used to compare with the ones obtained with this technique, illustrating the efficiency of the techniques developed and implemented in this work.

Acknowledgments

Arriving to this stage, I am glad to complete it by remembering many wonderful people who became involved with this thesis in various ways.

First and foremost I want to thank my thesis supervisors. I am extremely grateful for the guidance and energetic support provided by Professor Paulo Vasconcelos and for all the insightful and always fast feedback whenever required, despite his many other academic commitments. I deeply acknowledge Professor Mario Ahues for all the time he spent providing very important suggestions and comments when it was most necessary. I also thank him for his generous and tireless assistance in my stays in Saint-Etienne, as well as in all my registations at the Université Jean Monnet. It has been an honor to work with both of them, and I am deeply grateful to them for their scientific and general support and encouragement throughout my study as a PhD. student. I express my sincere thanks to Professor J. E. Roman from Universidad Politecnica de Valencia (Spain), for all his support in the joint computacional work held in this thesis. I wish to honor the memory of Professor Alain Largillier, with whom I had some scientific profitable talks about my work and helped to make my first stay in Saint-Etienne a pleasant experience. I am very grateful to my colleagues of Université Jean Monnet, specially to Rola Fares and Taline Boyajian, for the hospitality I experienced in my first stay there.

I also deeply appreciate all the support shown by my friends and IPCA colleagues.

iii iv I wish to express my warm thanks to my beloved husband, Paolo, for his love and understanding in every situation.

Last but not least, a very special thank to my wonderful parents for their endless love, patience and encouragement.

Without you all, I could not have reached this far. Thank you! Ana Luísa

Aproximação espectral com matrizes e operadores discretizados Resumo

Nesta dissertação considera-se a solução numérica de um problema de valores próprios de grandes dimensões, no qual o operador provém de um problema de transferência radiativa.

Procede-se ao estudo do uso de matrizes hierárquicas, uma representação eficiente de matrizes, bastante interessante para o uso em problemas de grandes dimensões. Matrizes hierárquicas são representações eficientes de estruturas de dados esparsas de matrizes densamente povoadas, sendo a ideia básica a de dividir uma determinada matriz numa hierarquia de blocos e aproximar determinados blocos por uma matriz de característica pequena. A sua utilização vem permitir, para além da diminuição da memória requerida, a redução dos custos computacionais.

A aplicação do uso das matrizes hierárquicas é analisada no contexto da solução numérica de um problema de valores próprios de grandes dimensões que resulta da discretização de um operador integral. O operador é de convolução e é definido através da primeira função exponencial-integral sendo, desta forma, fracamente singular.

Para o cálculo computacional, acede-se à HLIB (Hierarchical matrices LIBrary) que fornece rotinas para a construção da estrutura das matrizes hierárquicas, bem como algoritmos para operações aproximadas com estas viii matrizes. Acrescenta-se que se incorporam algumas rotinas da HLIB, como multiplicação matriz-vector ou a decomposição LU, na SLEPc (Hierarchical matrices LIBrary) de forma a explorar os algoritmos existentes para a resolução de problemas de valores próprios.

Desenvolvem-se ainda expressões analíticas para a aproximação dos núcleos degenerados utilizados na tese e deduzem-se também limites superiores de erros para estas aproximações.

Os resultados numéricos obtidos com outras abordagens para solucionar o problema em questão são utilizados para comparação com os obtidos com a nova técnica, vindo ilustrar a eficiência desta última.

Approximation spectrale de matrices issues d'opérateurs discrétisés Résumé

Cette thèse considère la solution numérique d'un problème aux valeurs propres de grandes dimensions, dans lequel l'opérateur est dérivé d'un problème de transfert radiatif.

Ainsi, cette thèse étudie l'utilisation de matrices hiérarchiques, une représentation efficace de tableaux, très intéressante pour une utilisation avec des problèmes de grandes dimensions. Les matrices sont des représentations hiérarchiques de structures de données efficaces pour les matrices denses, l'idée de base étant la division d'une matrice en une hiérarchie de blocs et l´approximation de certains blocs par une matrice de petite caractéristique. Son utilisation permet de diminuer la mémoire nécessaire tout en réduisant les coûts informatiques.

L'application de l'utilisation de matrices hiérarchique est analysée dans le contexte de la solution numérique d'un problème aux valeurs propres de grandes dimensions résultant de la discrétisation d'un opérateur intégral.

L'opérateur est de convolution et est défini par la première fonction exponentielle intégrale, donc faiblement singulière.

Pour le calcul informatique, nous avons accès à HLIB (Hierarchical matrices LIBrary) qui fournit des routines pour la construction de la structure hiérarchique des matrices et des algorithmes pour les opérations approxima-

List of Symbols

We now describe the main notations used in this thesis.

T integral operator

E i exponential-integral

Introduction

Eigenvalue problems are considered to be a very important subject of linear algebra, as they appear in many pratical applications in science and engineering. For example, in stability and control, eigenvalues give important information about damping, phase and magnitude of oscillations [START_REF] Makarov | A general method for small signal stability analysis[END_REF][START_REF] Makarov | Calculation of oscillatory stability margins in the space of power system controlled parameters[END_REF].

In [START_REF] Zhang | Towards an automatic and application-based eigensolver selection[END_REF] we can perceive some applications where the computation of eigenvalues is of crucial importance:

"The computation of eigenvalues and eigenvectors is an important and often time-consuming phase in computer simulations. In the last decades, very effective methods have been devised to compute some or all the eigenvalues of a matrix, such as the QR method for the latter case, a standard method for dense matrices of moderate size. With this sophisticated technique, the complete set of eigenvalues and eigenvectors is computed. Considering a real nonsingular square matrix A, an orthogonal matrix Q and an upper triangular matrix R, the QR method consists, briefly, in iterating the following steps: transform A into a hessenberg (tridiagonal if symmetric) matrix H (using e.g. Householder reflections), decompose H in QR, multiply Q and R together in reverse order (i.e., RQ) to form a new matrix H and the diagonal of H will converge to the eigenvalues, for more see [START_REF] Van Der | Computational methods for large eigenvalue problems[END_REF][START_REF] Golub | Matrix computations[END_REF].

Without

In this context, we refer to LAPACK (Linear Algebra PACKage) [ABB99],

a library with subroutines to solve the most frequent problems that come about in numerical linear algebra: linear equations, linear least squares problems, singular value decomposition and eigenvalue problems. LAPACK uses for computation calls as much as possible from the Basic Linear Algebra Subprograms (BLAS), that is a library for vector and matrix operations, adapted to the hierarchical memory of today's computers. LAPACK, along with general (dense) structures, has special implementations of the algorithms to deal with special structures, e.g. band matrices.

For large, usually sparse matrices, very effective iterative methods have been devised to approximate the eigenspace associated to any part of the spectrum. Examples of such techniques are the restarted Krylov methods and preconditioned eigensolvers such as Jacobi-Davidson [BDD + 00]. Also, these methods are progressively taking shape as high-quality implementations in software libraries such as SLEPc, the Scalable Library for Eigenvalue

Problem Computations [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF], thus enabling application programmers to cope with challenging problems coming from a wide range of applications.

Libraries such as SLEPc try to make problems computationally tractable by combining two main ingredients: (i) exploiting sparsity of the matrices, and

(ii) exploiting parallelism.

The sparsity of matrices is a desirable property that appears, for instance, in the context of partial differential equations with standard discretization techniques such as the finite element method. This situation is very common in practice, and allows iterative methods to be competitive by benefiting from the cheap, linear-cost matrix-vector products. However, there are cases in which the problem is formulated as an integral equation, either from the very nature of the problem or from some special treatment of PDEs such as the boundary element method.

Discretisation schemes like Galerkin or collocation methods are very popular techniques used for solving integral equations numerically. The integral equation is solved by these methods using a system of linear equations, but sparsity of matrices is not guaranteed, so in principle full (dense) storage must be used, making impossible the numerical treatment of large dimensional problems (consequent blow-up in computational cost). Furthermore, large memory requirements and large time consuming are expected even for moderate size ones. Some methods have been developed in order to avoid dealing with discretized dense matrices, for example, using compactly supported orthonormal wavelets to represent the integral operator [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF] or hierarchical matrices, H-matrices for brief, introduced by Hackbush and his collaborators [START_REF] Hackbusch | A sparse matrix arithmetic based on Hmatrices. Part I: Introduction to H-matrices[END_REF][START_REF] Hackbusch | A sparse H-matrix arithmetic. Part II: application to multi-dimensional problems[END_REF], being this last approach our choice in this thesis.

HLIB (Hierarchical matrices LIBrary) is a software that provides, among others, routines for the construction of hierarchical matrix structures and arithmetic algorithms to perform approximate matrix operations. The idea of using H-matrices, is that they provide an inexpensive representation of large densely populated matrices, by means of a decomposition into a hierarchy of blocks, in which most of those, usually far from the diagonal, are computed using a low-rank approximation in factorized form, resulting in a data-sparse storage scheme. It is meant by data-sparse that few data is needed for the matrix representation, that is, for an n × n-matrix, instead of n 2 entries, only O(nk log n) data may be required for a good approximation, where k is the rank of the new representation (it also determines the accuracy of the approximation) [START_REF] Börm | Construction of data sparse H 2 -matrices by hierarchical compression[END_REF].

The low-rank blocks' approximation is based on geometrical considerations, i.e., on an admissibility condition involving the notion of distance between two subsets of an index set; testing this condition, it is decided if the blocks of the adjacency matrix are to be approximated, which significantly reduces not only storage requirements but also computational costs.

The admissibility condition results naturally from the derivation of the error bound in order to produce convergent approximations.

Being SLEPc already an effcient tool in parallel computing for large and sparse eigenvalue problems, being sparsity of matrices a desirable property and H-matrices an excellent way to achieve that, this thesis accesses all the power of HLIB and incorporates the matrix-vector multiply routines from HLIB into SLEPc.

The main goal of this thesis is to solve an eigenvalue problem arising from an integral formulation of a radiative transfer problem in stellar atmospheres.

The kernel of the integral operator under study is weakly singular and the discretization procedure gives rise to large dimensional problems, either due to large integration limits or to fine grids. In the latter, the eigenproblem becomes very hard to solve since the eigenvalues of the discretized operator tend to be clustered. This thesis also aims to show the use of data-sparse structures provided by HLIB in the computation of eigenpairs through the SLEPc library. This is the first time such an integration is done and tested.

The thesis is divided into three parts, each one with two chapters and it is organized as follows. In Chapter 1 is presented the eigenvalue problem where approximate solution is sought, as well as a brief discription of hierarchical matrices along with the state-of-the-art numerical methods to solve large and sparse eigenvalue problems. Chapter 2 considers all the basics about the

Theoretical aspects

In this chapter we introduce the eigenvalue problem, its solution and the concept of hierarchical matrices.

The problem and its discretization

We consider an eigenvalue problem, arising from an integral formulation of a radiative transfer problem in stellar atmospheres [START_REF] Rutily | The finite Laplace transform for solving a weakly singular integral equation occurring in transfer theory[END_REF]. The spectral decomposition of an integral operator T has a numerical interest in the integral equation

x = f + T x, (1.1)
where x is the source function of the problem and f describes the distribution of internal and external sources [START_REF] Rutily | Multiple scattering theory and integral equations[END_REF]. T is a compact Fredholm integral operator from some Banach space X into itself, defined by

(T x) (τ ) := τ 0 g(τ, σ)x(σ)dσ, τ ∈ [0, τ ] , (1.2)
where τ < ∞ is the optical thickness of the atmosphere and g a weakly singular kernel. In particular, the kernel we will deal with is of the form 28 CHAPTER 1. THEORETICAL ASPECTS g(τ, σ) = k(|τ -σ|) and is weakly singular in the following sense [AdL + 02]:

lim τ →σ k(|τ -σ|) = +∞; k ∈ C 0 (]0, τ ]) ∩ L 1 ([0, τ ]); sup τ ∈[0,τ ] τ 0 k(|τ -σ|)dσ < +∞; k(|τ -σ|) > 0 for all τ, σ ∈]0, τ ];
k is a decreasing function on ]0, τ ].

We remark that for simplification, in what follows we will use the notation g(τ, σ) for the operator's kernel, and the one we will work with is defined as

g(τ, σ) := 2 E 1 (|τ -σ|) , (1.3)
where ∈ ]0, 1[ is the albedo, assumed as a constant. The kernel depends on E 1 , the first of a family of functions E ν , the exponential-integral functions [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], defined by

E ν (τ ) := ∞ 1 exp(-τ µ) µ ν dµ, τ > 0, ν ≥ 0, (1.4)
and

X := L 1 ([0, τ ]).
The application of the eigenvalues in physics is under study, since the absence of internal and external sources leads to the nullity of the radiation field (f = 0 results in x = 0 if τ < ∞, [I.W60]). The numerical application of the eigenvalues is present in the approximated sum of the Neumann series that produces the solution x of (1.1). The largest eigenvalue in magnitude, λ 1 , informs us about the speed of convergence of that series and the next greatest eigenvalue, λ 2 , about how quickly, starting in a certain term, can be held the replacement of part of the series terms by those of a geometric series of ratio λ 1 ; these ideas can be found in [dH80].

We consider the problem of finding (λ, x) ∈ C × X, such that

T x = λx = 0.
(1.5)

THE PROBLEM AND ITS DISCRETIZATION

29

Solving integral equations is of extreme importance, as they are recursively used to model some physical problems. In order to obtain approximate solutions to these kind of problems, the integral operator T is replaced by a finite rank operator T n . So, in our case, the eigenvalue problem for the integral operator is replaced by a matrix eigenvalue problem, which can be solved computationally.

Popular schemes have been around for some time and are used to solve integral equations, for instance, Nyström method or projection methods as Galerkin, Sloan or Kantorovich methods.

Before focusing in our problem, just a brief glimpse to these methods.

In the Nyström method, the integral of the operator is replaced by a numerical quadrature formula; this method is, in principle, applied for solving equations of the second kind, see [START_REF] Kress | Linear integral equations[END_REF]. Projection methods extend the previous method in the sense that the problem may be set in a wide class of Banach spaces such as L p spaces (1 ≤ p ≤ +∞), while Nyström method context is at least that of continuous functions. As examples we have:

• Galerkin Method: T n = π n T π n • Kantorovich Method: T n = π n T • Sloan Method: T n = T π n
where π n : X → X is a projection operator with finite dimensional range

X n ⊂ X and, for each x ∈ X, as n → ∞ π n x → x.
A more detailed explanation of the methods, their convergence and error analysis, may be found e.g. in [START_REF] Atkinson | The numerical solution of integral equations of the second kind[END_REF][START_REF] Kress | Linear integral equations[END_REF].

In this work, we apply a projection method, on a finite dimensional subspace X n := Span {e n,j : j = 1, ..., n} built as follows: let (τ n,j ) n j=0 on [0, τ ] CHAPTER 1. THEORETICAL ASPECTS be a family of grids, and e n,j (τ ) := 1 if τ n,j-1 < τ ≤ τ n,j , and e n,j (τ ) := 0 otherwise, and for x ∈ X define the linear forms

x, e * n,j := 1 τ n,j -τ n,j-1 τ n,j τ n,j-1

x(σ)dσ.

The bounded n-rank projection π n onto the subspace X n is defined by

π n x := n j=1
x, e * n,j e n,j , and, using the Kantorovich method, the finite rank approximation T n of T is defined by

T n x = π n T x = n j=1
T x, e * n,j e n,j .

The spectral problem for the finite rank operator T n can be solved through an auxiliary n × n matrix eigenvalue problem

A n x n = λ n x n = 0, (1.6)
where A n (i, j) := T e n,j , e * n,i , see [START_REF] Ahues | Defect correction for spectral computations for a singular integral operator[END_REF], and (λ n , x n ) is an approximation of (λ, x) for some x properly normalized.

The entries of this matrix are computed explicitly using (1.7), leading to a dense storage, and for large n the matrix's generation has a high computational cost.

A n (i, j) = 2(τ n,i -τ n,i-1 ) τ n,i τ n,i-1 τ n,j τ n,j-1 E 1 (|τ -σ|) dσdτ =        2(τ n,i -τ n,i-1 ) (-E 3 (|τ n,i -τ n,j |) + E 3 (|τ n,i-1 -τ n,j |)+ +E 3 (|τ n,i -τ n,j-1 |) -E 3 (|τ n,i-1 -τ n,j-1 |)), i = j [1 + 1 τ n,i -τ n,i-1 (E 3 (τ n,i -τ n,i-1 ) -1 2 )], i = j (1.7)
For each (i, j), four evaluations of the function E 3 are required. There is a clear decay in magnitude away from the diagonal, depending on τ and on n: for constant values of the former, smaller values of the latter imply faster

EIGENSOLVERS

decay from the diagonal. The idea of zeroing out all entries with magnitude less than a certain tolerance to avoid working with dense matrix storage, and a theoretical treatment validating this approach for a required precision, was done respectively in [START_REF] Almeida | A numerical study of iterative refinement schemes for weakly singular integral equations[END_REF] and [START_REF] Titaud | Reduction of computation in the numerical resolution of a second kind weakly singular Fredholm equation[END_REF]. Nevertheless, this strategy requires the computation of every matrix entry in order to evaluate its magnitude, since an a priori determination of the maximum bandwidth is not possible.

Besides the high generation cost, the main drawback of dense (or banded)

storage is that operations such as matrix-vector product are expensive. Datasparse representation, such as the general scheme of H-matrices, may be useful to tackle this kind of problems.

Eigensolvers

In this work we are concerned with the computation of a few eigenvalues and eigenvectors of matrix A n , that is, to obtain a partial solution of (1.6), by means of iterative eigensolvers. In the following subsections, we describe the methods very briefly, focusing on the required matrix operations that must be available in the implementation of H-matrices.

Iterative eigensolvers

There exist a lot of iterative methods for the partial solution of eigenvalue problems, that is, for computing a subset of the eigenvalues. A detailed discussion can be found in [BDD + 00]. Here we restrict our discussion to two families of methods, namely Krylov methods (e.g. Lanczos, Arnoldi or Krylov-Shur methods) and Davidson methods (e.g. Generalized Davidson or Jacobi-Davidson methods).

Given the eigenproblem

Ax = λx = 0, (1.8)
which has n eigenvalues λ counting multiplicities, the goal is to find a subset of the eigenvalues in a given region of the spectrum (for the moment, we CHAPTER 1. THEORETICAL ASPECTS consider the simplest case where we seek the largest magnitude eigenvalues).

Iterative eigensolvers are based on iteratively improving a subspace V in such a way that it eventually contains a good approximation of the eigenspace associated to the wanted eigenvalues.

Let V ∈ R n×m be a basis of V such that V T V = I. Then the Rayleigh-Ritz projection method computes H = V T AV and uses its eigendecomposition

HY = Y Θ to obtain approximate eigenpairs (θ i , x i = V y i ) of A. A more sophisticated alternative is the harmonic Rayleigh-Ritz method [BDD + 00],
that may provide better approximations in the case of interior eigenvalues.

Krylov methods use so-called Krylov subspaces associated with matrix A and a given initial vector v 1 ,

K m (A, v 1 ) = span{v 1 , Av 1 , A 2 v 1 , . . . , A m-1 v 1 } , (1.9)
where without loss of generality we assume that v 1 has unit length and is the first column of V .

The method of Arnoldi is an elegant algorithm that computes an orthonormal basis of the Krylov subspace and at the same time computes the projected matrix H, all this in an efficient and numerically stable way.

In brief, the Arnoldi algorithm computes the m columns of V sequentially

V m = [v 1 v 2 ... v m ],
where column v j+1 is the result of orthogonalizing Av j with respect to previous columns, and normalizing. The orthogonalization is carried out by means of a Gram-Schmidt procedure or other with better numerical properties, that removes all the components in the directions of v 1 , . . . , v j . The computed quantities satisfy a relation of the form

AV m = V m H m + βv m+1 e T m , (1.10)
where H m is an upper Hessenberg matrix, i.e., h ij = 0 for i > j + 1. The last term of the Arnoldi relation is the residual and gives an indication of how close is K m (A, v 1 ) to an invariant subspace. In particular, β is used to assess the accuracy of the computed Ritz pairs. See [BDD + 00] for additional details.

The Lanczos method is related to the Arnoldi method in the sense that Lanczos can be seen as a particular case of Arnoldi when the matrix is symmetric. In this case, the projected matrix is tridiagonal. For more detail see e.g. [BDD + 00, vdV02, Par98].

With Arnoldi, Ritz pairs will converge very fast provided that the initial vector v 1 is rich in the direction of the wanted eigenvectors. However, this is usually not the case and consequently many iterations will be required, but this cannot be allowed in a practical implementation in order to keep the storage requirements and the computational cost per iteration bounded. A workaround is to do a restart of the algorithm, that is, stop after m iterations and rerun the algorithm with a new v 1 computed from the recently obtained spectral approximations. An added benefit of this strategy is that it can be useful for driving convergence of the eigensolver towards a part of the spectrum different from the one targeted naturally by the method.

A very effective and elegant restart mechanism is the Krylov-Schur method [START_REF] Stewart | A Krylov-Schur algorithm for large eigenproblems[END_REF]. It is defined by generalizing the Arnoldi decomposition of order m,

(1.10), to a Krylov decomposition of order m,

AV m = V m B m + v m+1 b T m+1 , (1.11)
in which matrix B m is not restricted to be upper Hessenberg and b m+1 is an arbitrary vector. Krylov decompositions are invariant under orthogonal similarity transformations, so that

AV m Q = V m Q(Q T B m Q) + v m+1 b T m+1 Q, (1.12)
where

Q T Q = I, is also a Krylov decomposition. In particular, one can choose Q in such a way that S m = Q T B m Q is in (real) Schur form, that is,
upper (quasi-)triangular with the eigenvalues in the 1 × 1 or 2 × 2 diagonal blocks. This particular class of relation, called Krylov-Schur decomposition, can be written in block form as

A Ṽ1 Ṽ2 = Ṽ1 Ṽ2 S 11 S 12 0 S 22 + v m+1 bT 1 bT 2 , (1.13)
and has the nice feature that it can be truncated, resulting in a smaller Krylov-Schur decomposition,

A Ṽ1 = Ṽ1 S 11 + v m+1 bT 1 , (1.14)
that can be extended again to order m. The crux of the Krylov-Schur eigensolver is to carry out this truncation-extension process repeatedly, always keeping the wanted eigenvalues in the leading principal submatrix S 11 . Therefore, the strategy that is used for sorting the eigenvalues of S m will have an impact on which part of the spectrum will be approximated by the Krylov-Schur eigensolver.

So far, the only operation required with matrix A is the matrix-vector product, which can be carried out very efficiently in the H-matrix representation.

Computation of eigenvalues around a given target

The Krylov-Schur method could in principle be used to compute any part of the spectrum, by keeping the wanted eigenvalues in the truncated factorization. Discarding the rest of the factorization has the effect of filtering out the information associated to the unwanted eigenvectors. However, when computing eigenvalues in the interior of the spectrum, this filter is not powerful enough, and components associated to extreme eigenvalues keep on appearing, thus hindering convergence to the wanted ones.

The simplest solution to compute eigenvalues closest to a given target, σ, is to use a spectral transformation, in such a way that eigenvalues are mapped to a different position while eigenvectors remain unchanged. One such transformation is the shift-and-invert technique, that solves the problem

(A -σI) -1 x = θx, (1.15)
where the transformed eigenvalues satisfy the simple relation θ = (λ -σ) -1 .

In this transformation, the eigenvalues θ of the operator that are largest in 1.2. EIGENSOLVERS magnitude correspond, in the original problem, to the eigenvalues λ that are close, in absolute values, to σ, see Figure 1.1. As eigenvalues λ closest to the target become dominant in the transformed spectrum, so Krylov methods will have a fast convergence. This can be implemented by simply replacing the action of A by the action of (A-σI) -1 in the Krylov subspace expansion, that is, by solving linear systems with coefficient matrix A -σI whenever a matrix-vector product is required. These linear systems must be solved very accurately, since Krylov methods can be very sensitive to numerical errors introduced in the computation of the Krylov subspace, so in most applications a direct linear solver will be required, rather than an iterative method. It is generally claimed that the main drawback of the shift-andinvert technique is the high cost associated to direct linear solvers, since the memory requirements and computational effort can be very high for large, sparse matrices. In the case of the H-matrix representation, this downside disappears because computing the factorization has much smaller cost, both in terms of storage and operations, as well as the corresponding triangular solves.

An alternative to the spectral transformation is the use of a preconditioned eigensolver such as Jacobi-Davidson. These methods expand the subspace in a different way, attempting to make the whole computation more robust with respect to numerical error in the application of the operator. This allows to use iterative linear solvers such as GMRES (Generalized Minimal Residual) in the so-called correction equation.

Jacobi-Davidson method for linear problems was proposed by Sleijpen and van der Vorst [START_REF] Sleijpen | A Jacobi-Davidson iteration method for linear eigenvalue problems[END_REF] combining the Davidson method to expand the subspace, in which eigenvector approximations are constructed, with the Jacobi's idea, of looking for the orthogonal complement of a given eigenvector approximation. The method has been further developed and adapted to generalized eigenproblems [START_REF] Sleijpen | Jacobi-Davidson type methods for generalized eigenproblems and polynomial eigenproblems[END_REF]. This combined use is fundamental to enhance Davidson's idea. Whereas in Davidson's method accurate preconditioners lead to slow convergence (or even to stagnation), the Jacobi-Davidson This method was motivated by the fact that standard iterative eigensolvers often require an expensive factorization of the matrix when interior eigenvalues are desired (e.g., shift-and-invert Arnoldi with a direct linear solver). Jacobi-Davidson tries to reduce the cost by solving approximately linear systems, generally using iterative methods, without affecting the robustness.

This method usually is aiming at a particular eigenvalue, but if one is interested in more than one (near a specified target), a scheme presented in [START_REF] Fokkema | Jacobi-Davidson style QR and QZ algorithms for the reduction of matrix pencils[END_REF] can be used.

Nevertheless, this does not seem to be the best approach in the context of H-matrix representation, in view of the efficiency of matrix factorization.

From a practical perspective, to implement this kind of methods it is required to be able to build a preconditioner for matrix A.

H-MATRICES

H-matrices

Special mathematical and numerical treatments are required to find a suitable representation of the operator and to improve performance. As previously pointed out, H-matrices are an excellent way to tackle these requirements and overcome the difficulties.

Let I = {1, ..., n} denote the set of the indices of the basis functions e n,i , and t and s two subsets of I where α = i∈t supp e n,i and β = j∈s supp e n,j

are the corresponding domains.

The notions of a cluster tree and of a block cluster tree are important components in the construction of H-matrices. While the first describes a hierarchical partitioning over the index I giving us the candidates for checking an admissibility condition for low-rank approximation, the second contains the H-matrix's structure.

A cluster tree corresponding to the index set I, T I , satisfies the following properties, [BGH03b]:

• Each node of T I is a subset of I.

• The root of T I is the index set I.

• A leaf consists of a minor number of indices. A leaf is not more subdivided and this happens when the cardinality of a node is less than or equal to a certain threshold defined a priori.

• A node that is not a leaf is subdivided into two sons and is equal to their disjoint union.

The representation of H-matrices uses a tree structure, named block cluster tree and represented by T I×I . It describes a hierarchical block partitioning of a matrix. Its root represents the entire matrix and the inner nodes are the matrix sub-blocks that are being partitioned possibly more on the succeeding level, i.e., these blocks can be further subdivisible or not. In the last

CHAPTER 1. THEORETICAL ASPECTS {0, 1, 2, 3} y y % % {0, 1} | | % % {2, 3} y y " " {0} {1} {2} {3} Figure 1.2: Example of a T I .
case, the blocks, named leaves, are represented by either low-rank matrices or full matrices.

In order to ilustrate better these two structures, Figure 1.2 gives an example of a cluster tree T I with I = {0, 1, 2, 3} and, based upon T I of this figure, a block cluster tree is constructed in Figure 1.3.

Starting with I × I, the root of the so-called cluster tree T I×I , a splitting process begins through which each block is subdivided into four successors, until either an admissibility condition is satisfied or the block is already sufficiently small to be still subdivided. This last situation occurs when the refinement has arrived to the leaves, which possesses a certain minimal amount of elements defined a priori, as was previoulsy described.

By means of the following admissibility condition, which comes from error bounding reasons as will become apparent later in this thesis, we test if the

domain α × β is admissible diam(α) < ηdist(α, β), (1.16)
for η > 0 fixed, that is, if the corresponding block of indices t×s is admissible.

Let M | t×s be the corresponding submatrix. It will be approximated by a low- Following the previous figures of the cluster tree and block cluster tree with I = {0, 1, 2, 3}, Figure 1.4 appears to illustrate the respective H-matrix structure, where the red blocks correspond to the inadmissible leaves, while the green ones correspond to the admissible blocks. If the rank k is smaller than t and s, corresponding to the matrix M size, considerable savings in the storage and work complexity of a low-rank matrix compared to a full matrix are obtained [START_REF] Grasedyck | Construction and arithmetics of H-matrices[END_REF].

rank matrix M | t×s = AB T , where M | t×s ∈ R t×s , A ∈ R t×k , B ∈ R
For the low-rank approximation, the kernel of the integral operator is replaced by a degenerate approximation g(τ, σ), such that the integration with respect to the diferent variables is segregated,

g(τ, σ) := k-1 ρ=0 q ρ (τ )p ρ (σ).
(1.17) Furthermore, the approximantion g has, naturally, to converge fast to the kernel function g. In this work, the process of building the cluster tree is

H-MATRICES

undertaken by HLIB [START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF][START_REF] Börm | Hierarchical Matrices[END_REF] and will be detailed in Subsection 2.1 of this thesis.

The matrix entries

M ij = 1 h n,i τ 0 τ 0 e n,i (
τ )g (τ, σ) e n,j (σ)dσdτ for (i, j) ∈ t × s are then approximately given by

M | t×s ij = 1 h n,i τ 0 τ 0 e n,i (τ )g (τ, σ) e n,j (σ)dσdτ = 1 h n,i τ 0 τ 0 e n,i (τ ) k-1 ρ=0 q ρ (τ )p ρ (σ)e n,j (σ)dσdτ = k-1 ρ=0 1 h n,i τ 0 e n,i (τ )q ρ (τ )dτ τ 0 e n,j (σ)p ρ (σ)dσ, (1.18)
where

h n,i = τ n,i -τ n,i-1 .
We remark the importance of using (1.17) regarding the separation of the variables τ and σ, since it allows in (1.18) the possibility to write the double integral as a product of two single integrals, which are the entries of the matrices A and B of the factorized submatrix with rank at most k. So, the entries of those two matrices are defined by

A iρ := 1 h n,i τ n,i τ n,i-1 q ρ (τ )dτ, A = (A iρ ) ∈ R t×{0,...,k-1}
and (1.19)

B jρ := τ n,j τ n,j-1 p ρ (σ)dσ, B = (B jρ ) ∈ R s×{0,...,k-1} .
Chapter 2

Numerical libraries

This thesis makes use of a considerable set of software libraries which includes the most recent contributions of software developments based on novel mathematical approaches to numerical computation of eigenvalue problems.

These libraries are tuned to explore the modern architecture of today's computers. So, in this chapter we give a short presentation of each one of the main libraries used.

HLIB library

HLIB [START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF][START_REF] Börm | Hierarchical Matrices[END_REF] is a library for hierarchical matrices written in C programming language using BLAS and LAPACK libraries to perform lower-level algebraic operations (e.g. dense matrix-matrix multiplication). It includes functions for H matrix arithmetics, the treatment of partial differential equations and a number of integral operators and in addition support routines for the creation of cluster trees, visualization and numerical quadrature.

In the implementation of an H-matrix, called supermatrix, rkmatrices are a straightforward representation of low-rank submatrices while fullmatrices represent the submatrices corresponding to inadmissible leaves (usually dense matrices although not having to bear any specific structure). This implemen-

44

CHAPTER 2. NUMERICAL LIBRARIES tation may be found e.g. in [START_REF] Börm | Hierarchical Matrices[END_REF], and it is done in such a way that its structure is similar to the block cluster tree's structure.

This library provides routines for building the structure of an hierarchical matrix, that is, of cluster trees, block cluster trees, low-rank matrices and block matrices; discretization functions that fill these structures by approximations of operators, arithmetic algorithms that produce approximative matrix operations (e.g. addition, multiplication, vector-matrix multiplication, factorizations, inversion); conversion routines that turn sparse and dense matrices into H-matrices; service functions matrix structures, for instance to plot or to handle files. Moreover, based upon basic linear algebra subroutines, it also provides algorithms to compute LU-decompostion and Cholesky decomposition.

So, HLIB provide algorithms that perform matrix operations in the hierarchical matrix format efficiently and the actual proof for the efficiency, namely the complexity estimates, can be found, besides [START_REF] Börm | Hierarchical Matrices[END_REF], in [START_REF] Grasedyck | Construction and arithmetics of H-matrices[END_REF] where the authors analyse the complexity (storage, addition, multiplication and inversion) of the H-matrix arithmetics.

In this work, the process of building the cluster tree is undertaken by HLIB, which, as already mentioned, enables matrix operations of almost linear complexity, being therefore particularly adequate for large dimensional problems. Alternatively, AHMED (Another software library on Hierarchical Matrices for Elliptic Differential equations) could be used. A complete reference for H-matrices as well as for this library is [START_REF] Bebendorf | Hierarchical matrices: a means to efficiently solve elliptic boundary value problems[END_REF]. With the H-matrix representation, it is possible to realize common computations with linear-polylogarithmic complexity rather than quadratic or cubic cost.

For instance, the multiplication of an n × n H-matrix by a vector requires about 4n log 2 n floating-point operations, the (approximate) LU decomposition about 6n log 2 2n operations, and 2n log 2 n operations for the back solves.

See [START_REF] Hackbusch | A sparse matrix arithmetic based on Hmatrices. Part I: Introduction to H-matrices[END_REF][START_REF] Hackbusch | A sparse H-matrix arithmetic. Part II: application to multi-dimensional problems[END_REF] for additional details on H-matrix arithmetic.

PETSC AND SLEPC LIBRARIES

45

PETSc and SLEPc libraries

Working with sparse representation is not standard, as it happens in dense representation of an operator matrix. This comes from the fact that sparse storage is more complicated, as there may be more variations, so less standardized, and in this way more complex to work. An error bound for the computed solution in hermitian problems is available in literature,

λ -λ ≤ r 2 (2.1)
where r = Ax-λx is the residual vector, being λ, x the computed eigenpair and λ the exact eigenvalue. For the non-hermitian case, such simple relation as (2.1) is not available.

SLEPc is built on top of PETSc [BBE + 10], a parallel framework for the numerical solution of partial differential equations, whose approach is to encapsulate mathematical algorithms using object-oriented programming techniques in order to be able to manage the complexity of efficient numerical message-passing codes. It uses primarily the basic data structures such as those for representing vectors and matrices.

PETSc is object-oriented in the sense that all the code is built around a set of data structures and algorithmic objects. The application programmer works directly with these objects rather than concentrating on the underlying data structures. The three basic abstract data objects are index sets, vectors and matrices. Built on top of this foundation are various classes of solver objects, including linear, non-linear and time-stepping solvers.

SLEPc inherits all the good properties of PETSc, including portability to a wide range of parallel platforms, scalability to a large number of processors, and run-time flexibility giving full control over the solution process (one can for instance specify the solver at run time, or change relevant parameters

BLAS AND LAPACK LIBRARIES

such as the tolerance or the size of the subspace basis). Also, since PETSc provides a uniform interface with all of its linear solvers and a large family of preconditioners, it is possible to compare diverse combinations of method and preconditioner, just by their specification at execution time.

The solvers in PETSc (and SLEPc) have a data-structure neutral implementation. This means that the computation can be done with different matrix storage formats, and also even with a matrix that is not stored explicitly. By default, a matrix in PETSc is stored in a parallel compressed-row sparse format (called aij), where each processor stores a subset of rows.

Other formats include the symmetric variant (sbaij), where only the upper triangular part is stored, as well as the dense storage (both sequential and parallel).

For implementing a matrix-free solver with so-called shell matrices, the application programmer has to create one of such matrices and define its operations, by binding a user-defined subroutine for each operation. Only the operations required by the actual computation need to be set, so in the simplest case it is sufficient to implement the matrix-vector product.

For more advanced functionality, e.g., preconditioning, other operations are required as well. We use this feature to interface our code to HLIB.

BLAS and LAPACK libraries

Among the main numerical linear algebra packages, is certainly the LAPACK and BLAS libraries.

LAPACK (Linear Algebra PACKage) [ABB99] is a library with subroutines for solving the most usual problems appearing in numerical linear algebra, having been conceived to be efficient on various modern highperformance computers. It provides the following routines:

• Linear equations,

• Linear least squares problems,

• Eigenvalue problems,

• Singular value decomposition.

LAPACK can also manage many associated computations such as matrix factorizations or estimating condition numbers.

This library contains: driver routines for solving standard types of problems, each solving a complete problem, e.g., solving a system of linear equations; computational routines each performing a different computational task, e.g. LU factorization; and auxiliary routines to perform a certain subtask or common low-level computation, e.g. computing a matrix-norm. Note that each driver routine calls a sequence of computational routines.

LAPACK routines are design so that the computation is performed by calls, as much as possible, to the Basic Linear Algebra Subprograms (BLAS).

BLAS is a library that provides standard building blocks for performing basic vector and matrix operations. The Level 1 BLAS perform scalar, vector and vector-vector operations, the Level 2 BLAS perform matrix-vector operations, and the Level 3 BLAS perform matrix-matrix operations. A large set of research works, as well as industrial ones, rely on top of this library.

The first published paper on the library was [START_REF] Lawson | Basic linear algebra subprograms for fortran usage[END_REF] and, maybe, one of the most important references is [START_REF] Duff | An overview of the sparse basic linear algebra subprograms: The new standard from the BLAS technical forum[END_REF].

Highly efficient machine-specific implementations of the BLAS are available for many modern high-performance computers. The BLAS, being efficient, portable, and widely available, enable LAPACK routines to achieve high performance with portable code, allowing this to be high quality linear algebra software.

OpenMP library

The OpenMP standard was formulated in 1997 as an Application Program Interface (API) for writing portable, multithreaded applications, shared memory parallelism (note that a shared memory process consists of multiple

OPENMP LIBRARY

threads) [START_REF] Chapman | Using OpenMP: portable shared memory parallel programming[END_REF]. It is a specification for a set of compiler directives, runtime library routines, and environment variables that can be used to specify shared memory parallelism in Fortran, C and C++ programs, on a wide variety of architectures.

OpenMP provides a portable (most major platforms have been implemented including Unix/Linux platforms and Windows NT), scalable model for developers of shared memory parallel applications, in which the most of its major features includes its various constructs and directives for specifying parallel regions, work sharing, synchronization and data environment.

This library is the result of efforts for some standardization, providing a standard among a variety of shared memory architectures/platforms, seting up a simple and limited set of directives for programming shared memory machines, be ease to use, have capability to incrementally parallelize a serial program as well as capability to implement both coarse-grain and fine-grain parallelism, and portability.

The OpenMP is based on the existence of multiple threads in the shared memory programming paradigm, being an explicit (i.e. not automatic) programming model allowing the programmer full control over parallelization.

An OpenMP application begins with a single thread (master thread). As the program executes, the application may be faced with parallel regions in which the master thread gives rise to thread teams, including the master thread. At the end of a parallel region, the thread teams are "left" and the master thread continues execution. Inside a parallel region there can be nested parallel regions where each thread of the original parallel region becomes the master of its own thread team. Nested parallelism can continue to additionally nest other parallel regions.

Moreover, OpenMP provides the distinction between data that is shared from data that is private, an important issue in parallel programming (performance). In this context, shared variables are shared by all the threads from the thread team (e.g. a change of the shared variable in one thread may become visible to another thread in the parallel region), while in private Chapter 3

Low-rank approximations

Previously, in (1.7), it was given a dense representation for the matrix A n .

In this chapter, we discuss several strategies for representing A n by means of data-sparse representation, which can be implemented in various flavours.

Degenerate approximations, generally given by (1.17), can be built in different ways, depending on the properties of the kernel. For instance, they can be based on polynomial interpolation or Taylor series expansion. Another approach consists in computing a low-rank approximation from an explicitly built matrix block by means of a singular value decomposition. We next describe these three approaches in turn.

Taylor approximation

Using truncated Taylor series, the degenerate kernel expression (1.17) is specifically defined by q ρ (τ ) = (τ -τ 0 ) ρ and p ρ (σ) = 1 ρ! ∂ ρ τ g(τ 0 , σ), where τ 0 is taken as the midpoint of α.

The entries of the matrices A and B in (1.19) are analytically determined explicitly as follows.

Proposition 2. The elements of the matrix A are given by

A iρ = 1 h n,i τ n,i τ n,i-1 q ρ (τ )dτ = 1 h n,i (τ n,i -τ 0 ) ρ+1 ρ + 1 - (τ n,i-1 -τ 0 ) ρ+1 ρ + 1 . (3.1)

CHAPTER 3. LOW-RANK APPROXIMATIONS

Proof. Immediate.

In order to get the expression for g and for the entries B jρ , a preceding step must be made, to derive the analytical expression for ∂ ρ τ g(τ 0 , σ).

Lemma 3. The partial derivatives with respect to the variable τ of the kernel g defined in (1.3), are given by

∂ ρ τ g(τ 0 , σ) =    -2 ρ k=1 (-1) ρ-1 (ρ-1)! (ρ-k)! e -τ 0 +σ (τ 0 -σ) k , σ < τ 0 -2 ρ k=1 (-1) k-1 (ρ-1)! (ρ-k)! e τ 0 -σ (τ 0 -σ) k , σ > τ 0 . Proof. Considering the following results present in [AS60], E 1 (z) = ∞ 1 e -zt t dt and d dz E 1 (z) = -E 0 (z) = - e -z z ,
the first derivative of the kernel with respect to τ is given by

∂ τ g(τ, σ) = 2 ∂ τ E 1 (|τ -σ|) = 2 ∂ τ (E 1 • f (τ, σ)) where f (τ, σ) = |τ -σ| = 2 - e -f (τ,σ) f (τ, σ) × ∂ τ f (τ, σ) = 2 - e -|τ -σ| |τ -σ| × ∂ τ f (τ, σ) = -2 e -(τ -σ) τ -σ , σ < τ -2 e -(-τ +σ) τ -σ , σ > τ = - 2 e -|τ -σ| τ -σ , for τ = σ.
1. For σ > τ , after some initial calculations, we obtain

∂ 2 τ g(τ, σ) = - 2 ∂ 1 τ e τ -σ τ -σ = - 2 e τ -σ (τ -σ) -e τ -σ (τ -σ) 2 = - 2 e τ -σ τ -σ - e τ -σ (τ -σ) 2 , 3.1. TAYLOR APPROXIMATION 55 ∂ 3 τ g(τ, σ) = - 2 ∂ 1 τ e τ -σ τ -σ - e τ -σ (τ -σ) 2 = - 2 e τ -σ τ -σ - e τ -σ (τ -σ) 2 -∂ 1 τ e τ -σ (τ -σ) 2 = - 2 e τ -σ τ -σ - 2e τ -σ (τ -σ) 2 + 2e τ -σ (τ -σ) 3 ,
and, generalizing, we reach the expression for any value of ρ

∂ ρ τ g(τ, σ) = - 2 ρ k=1 (-1) k-1 × (ρ -1)! (ρ -k)! × e τ -σ (τ -σ) k .
2. For σ < τ , analogously to the previous case, the following general recursive rule is obtained

∂ ρ τ g(τ, σ) = - 2 ρ k=1 (-1) ρ-1 (ρ -1)! (ρ -k)! e -τ +σ (τ -σ) k .
Using the previous proposition, the expression for the truncated Taylor series kernel is explicitly given by

g(τ, σ) = k-1 ρ=0 q ρ (τ )p ρ (σ) =    -2 k-1 ρ=0 (τ -τ 0 ) ρ 1 ρ ρ i=1 (-1) ρ-1 1 (ρ-i)! e -τ 0 +σ (τ 0 -σ) i , σ < τ 0 -2 k-1 ρ=0 (τ -τ 0 ) ρ 1 ρ ρ i=1 (-1) i-1 1 (ρ-i)! e τ 0 -σ (τ 0 -σ) i , σ > τ 0 . (3.2)
With the result of Lemma 3, we can readily obtain the expression for the entries of matrix B, by evaluating the derivatives at τ = τ 0 .

Proposition 4. The entries of matrix B can be computed as

B jρ =                  -2 ρ k=1 (-1) ρ-1 ρ(ρ-k)! [Γ (1 -k, τ 0 -τ n,j ) -Γ (1 -k, τ 0 -τ n,j-1 )] , τ n,j-1 < τ 0 , τ n,j < τ 0 2 ρ k=1 1 ρ(ρ-k)! [Γ (1 -k, -τ 0 + τ n,j-1 ) -Γ (1 -k, -τ 0 + τ n,j )] , τ n,j-1 > τ 0 , τ n,j > τ 0 A n (i, j) 1 , τ j-1 < τ 0 < τ j 56 CHAPTER 3. LOW-RANK APPROXIMATIONS Proof. First recall the incomplete Gamma function Γ (a, x) = ∞ x
e -t t 1-a dt. 1. For τ n,j-1 > τ 0 , τ n,j > τ 0 :

B jρ = τ n,j τ n,j-1 1 ρ! ∂ ρ τ g(τ 0 , σ)dσ = τ n,j τ n,j-1 1 ρ! - 2 ρ k=1 (-1) k-1 (ρ -1)! (ρ -k)! e τ 0 -σ (τ 0 -σ) k dσ = - 2 ρ k=1 (-1) k-1 ρ (ρ -k)! τ n,j τ n,j-1 e τ 0 -σ (-1) k (-τ 0 + σ) k dσ
and using integration by substitution (u = -τ 0 + σ),

B jρ = - 2 ρ k=1 (-1) k-1 ρ (ρ -k)! × 1 (-1) k -τ 0 +τ n,j -τ 0 +τ n,j-1 e -u u k du = 2 ρ k=1 1 ρ (ρ -k)! × ∞ -τ 0 +τ n,j-1 e -u u k du - ∞ -τ 0 +τ n,j e -u u k du = 2 ρ k=1 1 ρ (ρ -k)! × [Γ (1 -k, -τ 0 + τ n,j-1 ) -Γ (1 -k, -τ 0 + τ n,j )] .
2. For τ n,j-1 < τ 0 , τ n,j < τ 0 :

B jρ = τ n,j τ n,j-1 1 ρ! ∂ ρ τ g(τ 0 , σ)dσ = τ n,j τ n,j-1 1 ρ! - 2 ρ k=1 (-1) ρ-1 (ρ -1)! (ρ -k)! e -τ 0 +σ (τ 0 -σ) k dσ = - 2 ρ k=1 (-1) ρ-1 ρ (ρ -k)! τ n,j τ n,j-1 e -(τ 0 -σ) (τ 0 -σ) k dσ doing the substitution u = τ 0 -σ, B jρ = - 2 ρ k=1 (-1) ρ-1 ρ (ρ -k)! - τ 0 -τ n,j τ 0 -τ n,j-1
e -u u k du , and as

τ 0 -τ n,j-1 > τ 0 -τ n,j , = - 2 ρ k=1 (-1) ρ-1 ρ (ρ -k)! τ 0 -τ n,j-1 τ 0 -τ n,j e -u u k du = - 2 ρ k=1 (-1) ρ-1 ρ (ρ -k)! [Γ (1 -k, τ 0 -τ n,j ) -Γ (1 -k, τ 0 -τ n,j-1 )] .

ERROR BOUNDS FOR TAYLOR APPROXIMATION

57 3. It is compulsory to comment about the entries satisfying τ n,j-1 < τ 0 and τ n,j > τ 0 . Although in this situation the previous integral is divergent, those entries do not satisfy the admissibility condition, which implies that they must be computed using (1.7).

Proposition 4 can be written using exponential integral functions instead of the incomplete Gamma function, as it is next explained.

Proposition 5.

B jρ =                -2 ρ k=1 (-1) ρ-1 ρ(ρ-k)! (τ 0 -τ n,j ) 1-k E k (τ 0 -τ n,j ) - -(τ 0 -τ n,j-1 ) 1-k E k (τ 0 -τ n,j-1 ) , τ n,j-1 < τ 0 , τ n,j < τ 0 2 ρ k=1 1 ρ(ρ-k)! (-τ 0 + τ n,j-1 ) 1-k E k (-τ 0 + τ n,j-1 ) - -(-τ 0 + τ n,j ) 1-k E k (-τ 0 + τ n,j ) , τ n,j-1 > τ 0 , τ n,j > τ 0
Proof. Immediate, using the next equality present in [AS60]:

E n (x) := x n-1 Γ (1 -n, x) .

Error bounds for Taylor approximation

In order to develop error bounding proofs, we introduce some technical constants:

Definition 6. Let 1. κ be any constant such that ∀τ ∈ α, |τ -τ 0 | ≤ diam(α) ≤ κ < 1. (3.3)
We recall that the Taylor expansion is done around τ 0 .

CHAPTER 3. LOW-RANK APPROXIMATIONS

2. γ 0 be given by

γ 0 := diam(α) dist(α, β) < 1. (3.4)
We recall that

∀τ ∈ α, ∀σ ∈ β, dist(α, β) ≤ |τ -σ|. (3.5)
Result (3.4) can be rewritten as diam(α) < dist(α, β) which is, actually, the admissibility condition given in (1.16) (with η = 1). The parameter η controlls the speed of convergence, i.e., the quality of the approximation, and a usual choice is η = 1, see e.g. [START_REF] Börm | Introduction to hierarchical matrices with applications[END_REF][START_REF] Hackbusch | A sparse matrix arithmetic based on Hmatrices. Part I: Introduction to H-matrices[END_REF][START_REF] Hackbusch | Hierarchical matrices based on a weak admissibility criterion[END_REF].

In terms of notation, consider k = g(τ, σ) -g(τ, σ) to be the error when using g as an approximation of g.

In the next proposition we propose an upper bound for the error when g is taken as the truncated Taylor expansion of the kernel.

Proposition 7. The error of using the truncated Taylor series g, defined by

(3.2) as an approximation of the kernel g, given in (1.3), can be estimated by

| k | ≤        C 0 κ k 1-κ , |τ 0 -σ| > 1 C 1 κ k+1 (1-κ k ) 2 + kκ k 1-κ k , |τ 0 -σ| = 1 C 1 γ k 0 (k+γ 0 -kγ 0 ) (1-γ 0 ) 2 , |τ 0 -σ| < 1
where k ∈ N, κ is defined in (3.3), C 0 := 2 e -|τ 0 -σ| (|τ 0 -σ| -1) -1 and

C 1 := 2 e -|τ 0 -σ| .
Proof. Bearing in mind (3.2), we can write,

| k | ≤ ∞ ρ=k (τ -τ 0 ) ρ - 2 ρ i=1 1 (ρ -i)! e -|τ 0 -σ| (τ 0 -σ) i ≤ ∞ ρ=k 2 |τ -τ 0 | ρ ρ i=1 e -|τ 0 -σ| (τ 0 -σ) i = ∞ ρ=k 2 |τ -τ 0 | ρ 1 e |τ 0 -σ| ρ i=1 1 |τ 0 -σ| i
Now three different cases must be studied:

1. If |τ 0 -σ| > 1 | k | = ∞ ρ=k 2 |τ -τ 0 | ρ 1 e |τ 0 -σ| 1 - 1 |τ 0 -σ| ρ+1 1 -1 |τ 0 -σ| -1 ≤ ∞ ρ=k 2 |τ -τ 0 | ρ 1 e |τ 0 -σ| 1 1 -1 |τ 0 -σ| -1 = ∞ ρ=k 2 |τ -τ 0 | ρ e -|τ 0 -σ| (|τ 0 -σ| -1) -1
and considering C 0 := 2 e -|τ 0 -σ| (|τ 0 -σ| -1) -1 it results

| k | = C 0 ∞ ρ=k |τ -τ 0 | ρ = C 0 ∞ ρ=0 |τ -τ 0 | ρ - k-1 ρ=0 |τ -τ 0 | ρ = C 0 1 1 -|τ -τ 0 | - 1 -|τ -τ 0 | k 1 -|τ -τ 0 | = C 0 |τ -τ 0 | k 1 -|τ -τ 0 | ≤ C 0 κ k 1 -κ 2. If |τ 0 -σ| = 1 | k | ≤ ∞ ρ=k 2 |τ -τ 0 | ρ 1 e |τ 0 -σ| ρ = 2e |τ 0 -σ| ∞ ρ=k ρ|τ -τ 0 | ρ ≤ 2e |τ 0 -σ| ∞ ρ=k ρκ ρ
Bearing in mind the two following geometric series and respective sums,

• S(κ) = ∞ ρ=0 κ ρ = 1 1-κ • S (κ) = ∞ ρ=1 ρκ ρ-1 = 1 (1-κ) 2 and remarking that ∞ ρ=1 ρκ ρ-1 = ∞ ρ=0 ρκ ρ-1
we take ρ = k + j,

| k | ≤ 2e |τ 0 -σ| ∞ j=0 (k + j)κ k+j = 2e |τ 0 -σ| ∞ j=0 jκ k+j + ∞ j=0 kκ k+j = 2e |τ 0 -σ| κ k+1 ∞ j=0 jκ j-1 + kκ k ∞ j=0 κ j = 2e |τ 0 -σ| κ k+1 (1 -κ) 2 + kκ k 1 -κ = C 1 κ k+1 (1 -κ) 2 + kκ k 1 -κ where C 1 := 2 e -|τ 0 -σ| . 3. If |τ 0 -σ| < 1 | k | ≤ 2e |τ 0 -σ| ∞ ρ=k |τ -τ 0 | ρ ρ i=0 1 |τ 0 -σ| i -1 = C 1 ∞ ρ=k |τ -τ 0 | ρ ρ i=1 1 |τ 0 -σ| i ≤ C 1 ∞ ρ=k |τ -τ 0 | ρ ρ |τ 0 -σ| ρ = C 1 ∞ ρ=k ρ |τ -τ 0 | |τ 0 -σ| ρ using (3.5) and (3.3), | k | ≤ C 1 ∞ ρ=k ρ diam(α) dist(α, β)
| k | = C 1 ∞ ρ=k ργ ρ 0 = C 1 ∞ ρ=0 ργ ρ 0 - k-1 ρ=0 ργ ρ 0 = C 1 γ 0 ∞ ρ=0 ργ ρ-1 0 -γ 0 k-1 ρ=0 ργ ρ-1 0 = C 1 γ 0 ∞ ρ=0 d dγ 0 γ ρ 0 -γ 0 k-1 ρ=0 d dγ 0 γ ρ 0 = C 1 γ 0 d dγ 0 ∞ ρ=0 γ ρ 0 -γ 0 d dγ 0 k-1 ρ=0 γ ρ 0 = C 1 γ 0 d dγ 0 1 1 -γ 0 -γ 0 d dγ 0 1 -γ k 0 1 -γ 0 = C 1 γ 0 1 (1 -γ 0 ) 2 -γ 0 -kγ k-1 0 + kγ k 0 + 1 -γ k 0 (1 -γ 0 ) 2 = C 1 γ 0 (1 -γ 0 ) 2 - γ 0 -kγ k 0 + (k -1)γ k+1 0 (1 -γ 0 ) 2 = C 1 γ k 0 (k + γ 0 -kγ 0 ) (1 -γ 0 ) 2 .
We remark that the three error bounds have exponential convergence.

Interpolation approximation

Since Taylor expansion involves the computation of the derivative of the kernel function and a recursive rule for it, only seldom can be developed.

Polynomial interpolation appears as a good alternative; it requires only the evaluations of the kernel function.

Using polynomial interpolation, the degenerate kernel expression (1.17) is now given by q ρ (τ ) = L ρ (τ ) and p ρ (σ) = g (x ρ , σ).

Consider (x ρ ) k-1 ρ=0 a family of interpolation points in the subdomain α (for β the procedure is similar) and (L ρ ) k-1 ρ=0 the corresponding Lagrange 62 CHAPTER 3. LOW-RANK APPROXIMATIONS polynomials:

L ρ (x) = k-1 θ=0,θ =ρ x -x θ x ρ -x θ ,
for all x in an interval [a 1 , a 2 ]. For this interval, the set of interpolant points (x ρ ) k ρ=0 can be chosen as the Chebyshev points of order k -1:

x ρ := a 2 + a 1 2 + a 2 -a 1 2 cos 2ρ + 1 2k π .
To minimize the error in the approximation, for the reason to be soon explained, the degenerate kernel is chosen as

g (τ, σ) = k-1 ρ=0 g (x ρ , σ) L ρ (τ ) , if diam(α) ≤ diam(β), k-1 ρ=0 g (τ, x ρ ) L ρ (σ) , otherwise, (3.6)
and the admissibility condition (1.16) is adapted to

min{diam(α), diam(β)} < ηdist(α, β). (3.7)
This implies that the computation of the entries of matrices A and B from (1.19) is now done as follows: if diam(α) ≤ diam(β),

A iρ = 1 h n,i τ n,i τ n,i-1 L ρ (τ ) dτ and B jρ = τ n,j τ n,j-1 g (x ρ , σ) dσ, while if diam(β) ≤ diam(α), A iρ = 1 h n,i τ n,i τ n,i-1 g (τ, x ρ ) dτ and B jρ = τ n,j τ n,j-1 L ρ (σ) dσ.
All integrals can be obtained using numerical quadrature formulae, nevertheless, it is possible to have a probable more accurate result by computing analytically the integrals involving the kernel. For that, first recall the result in [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF],

dE n (z) dz = -E n-1 (z) for n ∈ {1, 2, 3, ...} leading immediately to E n-1 (z)dz = -E n (z). (3.8)
Proposition 8. When diam(β) ≤ diam(α), the entries A iρ are explicitly given by

A iρ =        -2h n,i [E 2 (τ n,i -x ρ ) -E 2 (τ n,i-1 -x ρ )] , τ n,i-1 ≥ x ρ , τ n,i ≥ x ρ 2h n,i [E 2 (-τ n,i + x ρ ) -E 2 (-τ n,i-1 + x ρ )] , τ n,i-1 < x ρ , τ n,i < x ρ 2h n,i [-E 2 (-τ n,i-1 + x ρ ) -E 2 (τ n,i -x ρ ) + 2] , τ n,i-1 < x ρ , τ n,i ≥ x ρ Proof. Recalling that g (τ, x ρ ) = 2 E 1 (|τ -x ρ |) =        2 E 1 (τ -x ρ ) , τ ≥ x ρ 2 E 1 (-τ + x ρ ) , τ < x ρ using (3.8
) and integration by substitution, we get

1. if τ n,i-1 ≥ x ρ and τ n,i ≥ x ρ , A iρ = 2h n,i τ n,i τ n,i-1 E 1 (τ -x ρ ) dτ = - 2h n,i [E 2 (τ n,i -x ρ ) -E 2 (τ n,i-1 -x ρ )] 2. if τ n,i-1 < x ρ and τ n,i < x ρ , A iρ = 2h n,i τ n,i τ n,i-1 E 1 (-τ + x ρ ) dτ = 2h n,i [E 2 (-τ n,i + x ρ ) -E 2 (-τ n,i-1 + x ρ )] 3. if τ n,i-1 < x ρ and τ n,i ≥ x ρ , A iρ = 2h n,i xρ τ n,i-1 E 1 (-τ + x ρ ) dτ + τ n,i xρ E 1 (τ -x ρ ) dτ = 2h n,i [-E 2 (-τ n,i-1 + x ρ ) -E 2 (τ n,i -x ρ ) + 2] CHAPTER 3. LOW-RANK APPROXIMATIONS
The entries of the matrix B are given in Proposition 9.

Proposition 9. When diam(α) ≤ diam(β), the entries of the matrix B are explicitly given by

B jρ =        2 [E 2 (x ρ -τ n,j ) -E 2 (x ρ -τ n,j-1 )] , τ n,j-1 ≤ x ρ , τ n,j ≤ x ρ -2 [E 2 (-x ρ + τ n,j ) -E 2 (-x ρ + τ n,j-1 )] , τ n,j-1 > x ρ , τ n,j > x ρ 2 [-E 2 (x ρ -τ n,j-1 ) -E 2 (-x ρ + τ n,j ) + 2] , τ n,j-1 ≤ x ρ , τ n,j > x ρ
Proof. Using a similar procedure to the one done in Proposition 8, we obtain

1. if τ n,j-1 ≤ x ρ and τ n,j ≤ x ρ , B jρ = 2 τ n,j τ n,j-1 E 1 (x ρ -σ) dσ taking u = x ρ -σ, = 2 xρ-τ n,j xρ-τ n,j-1 E 1 (u) (-du) and as x ρ -τ n,j ≤ x ρ -τ n,j-1 , = 2 xρ-τ n,j-1 xρ-τ n,j E 1 (u) du = 2 [-E 2 (x ρ -τ n,j-1 ) + E 2 (x ρ -τ n,j )] = 2 [E 2 (x ρ -τ n,j ) -E 2 (x ρ -τ n,j-1 )] 2. if τ n,j-1 > x ρ and τ n,j > x ρ , B jρ = 2 τ n,j τ n,j-1 E 1 (-x ρ + σ) dσ = - 2 [E 2 (-x ρ + τ n,j ) -E 2 (-x ρ + τ n,j-1 )] 3. if τ n,j-1 ≤ x ρ and τ n,j > x ρ , B jρ = 2 xρ τ n,j-1 E 1 (x ρ -σ) dσ + τ n,j xρ E 1 (-x ρ + σ) dσ = 2 lim x→xρ x τ n,j-1 E 1 (x ρ -σ) dσ + lim y→xρ τ n,j y E 1 (-x ρ + σ) dσ 3.4. ERROR BOUNDS FOR INTERPOLATION APPROXIMATION 65 taking u = x ρ -σ and v = -x ρ + σ, B jρ = 2 lim x→xρ xρ-x xρ-τ n,j-1 E 1 (u) (-du) + lim y→xρ τ n,j -xρ y-xρ E 1 (v) dv and as x ρ -x ≤ x ρ -τ n,j-1 , B jρ = 2 lim x→xρ xρ-τ n,j-1 xρ-x E 1 (u) du + lim y→xρ τ n,j -xρ y-xρ E 1 (v) dv = 2 lim x→xρ [-E 2 (u)] xρ-τ n,j-1 xρ-x + lim y→xρ [-E 2 (v)] τ n,j -xρ y-xρ = 2 [-E 2 (x ρ -τ n,j-1 ) + E 2 (0) -E 2 (-x ρ + τ n,j ) + E 2 (0)] = 2 [-E 2 (x ρ -τ n,j-1 ) -E 2 (-x ρ + τ n,j ) + 2]

Error bounds for interpolation approximation

Inspired by [START_REF] Börm | Low-rank approximation of integral operators by interpolation[END_REF], we give upper bounds for the error when polynomial interpolation is chosen to obtain the degenerate approximation of the kernel function.

The next Lemma introduces an upper bound for the derivative of the kernel of the integral operator we are dealing with.

Lemma 10. The derivative of the kernel function (1.3) with respect to the first argument of the kernel is bounded by

∂ ρ g ∂τ ρ (τ, σ) ≤ 2 (ρ -1)! 1 e dist(α,β) ρ k=1 1 [dist(α, β)] k .
(3.9) CHAPTER 3. LOW-RANK APPROXIMATIONS Proof. From Lemma 3 and using (3.5)

∂ ρ g ∂τ ρ (τ, σ) ≤ - 2 ρ k=1 (ρ -1)! (ρ -k)! e -|τ -σ| (τ -σ) k ≤ 2 (ρ -1)! 1 e |τ -σ| ρ k=1 1 |τ -σ| k ≤ 2 (ρ -1)! 1 e dist(α,β) ρ k=1 1 [dist(α, β)] k
In the next Lemma, an expression for the error bound when the interpolation points are the zeros of the Chebyshev polynomial is given: Lemma 11. In [START_REF] Süli | An introduction to numerical analysis[END_REF], for f ∈ C k+1 [a, b] and using Chebyshev interpolation, a bound for the approximation error is given by

k ∞ ≤ (b -a) k+1 2 2k+1 (k + 1)! f (k+1) ∞ .
(3.10)

Applying interpolation distinctly to the different arguments of the kernel, the error bound for the approximation of the kernel g by its interpolant g, when interpolation is done with respect to the first argument of the kernel (given by the following proposition (12)), is majorized.

Proposition 12.

(k-1)α ≤    e dist(α,β) k4 k 1-diam(α) k 1-diam(α) , diam(α) < 1 or 1 < diam(α) < 4 e dist(α,β) 4 k , diam(α) = 1

ERROR BOUNDS FOR INTERPOLATION APPROXIMATION 67

Proof. Using (3.10), (3.9) and (1.16)

(k-1)α ≤ diam(α) k 2 2k-1 k! ∂ k g ∂τ k (τ, σ) ≤ diam(α) k 2 2k-1 k! (k -1)! 2e dist(α,β) k i=1 1 [dist(α, β)] i = diam(α) k 4 k k 1 e dist(α,β) k i=1 1 [dist(α, β)] i = e dist(α,β) k4 k k i=1 diam(α) k [dist(α, β)] i ≤ e dist(α,β) k4 k k-1 i=0 diam(α) i =    e dist(α,β) k4 k 1-diam(α) k 1-diam(α) , diam(α) = 1 e dist(α,β) 4 k , diam (α) = 1
But to ensure convergence when k tends to infinity, in the case diam(α) = 1, regarding the presence of the factors 4 k and diam(α) k , respectively in the denominator and numerator, we must impose the additional condition diam(α) < 4.

Similarly, when interpolation is done with respect to the second argument of the kernel, we get the following upper bound to the approximation error.

Proposition 13.

(k-1) β ≤    e dist(α,β) k4 k 1-diam(β) k 1-diam(β) , diam(β) < 1 or 1 < diam(β) < 4 e dist(α,β) 4 k , diam(β) = 1
Proof. Analogous to Proposition 12.

It is clear now that, in order to minimize the approximating error, interpolation should be applied distinctively to both arguments of the kernel g, taking into account the diameter of α and β, that is why for polynomial interpolation the kernel is written as (3.6). Still related with the two previous propositions, we can state CHAPTER 3. LOW-RANK APPROXIMATIONS Proposition 14. For D = min{diam(α), diam(β)}, we get

| k-1 | ≤    e dist(α,β) k4 k 1-D k 1-D , D < 1 or 1 < D < 4 e dist(α,β) 4 k , D = 1
Proof. Take the minimum of diam(α) and diam(β) to minimize the error bound.

This justifies the changing of the left hand side of the admissibility condition (1.16) into (3.7) in the present approach, that is, for polynomial interpolation; the replacement with D, defined in Proposition 14, takes place for minimization of the error bound.

SVD approximation

This approach consists in computing a low-rank approximation from an explicitly built matrix block by means of a singular value decomposition (SVD).

A SVD on each admissible block, with already generated entries, can be used to preserve the most valuable information and discard the rest. The resulting rank-k approximation can be expressed as

AB T = k ρ=1 u ρ σ ρ v T ρ , (3.11)
where σ ρ are the singular values (in descending order of magnitude), and u ρ and v ρ are the corresponding left and right singular vectors, respectively.

Note that in the SVD, the eigenvectors of (AB T ) T (AB T ) are the right singular vectors of AB T , the eigenvectors of (AB T )(AB T ) T are the left singular vectors of AB T and the eigenvalues of (AB T ) T (AB T ) (or (AB T )(AB T ) T ) are the squares of the singular values of AB T .

The rank, k, can be chosen to be a fixed value, or alternatively to be set dynamically on each block, based on a prescribed tolerance . In the latter case, the condition σ k > ≥ σ k+1 holds.

SVD APPROXIMATION 69

However, in computational terms, the cost to obtain the H-matrix is high since the entries of every admissible block must be explicitly generated first, in the present case through (1.7), prior to the singular value decomposition.

Nevertheless, all the computations done afterwards can be performed cheaply.

This approach may only be of pratical use if data can be stored, after being generated, for subsequent (and repetitive) use.

the rank and minimum size of admissible blocks (degree=6 and bound=80, respectively), as well as η = 1 in (1.16).

For building the H-matrix we have used HLIB's supermatrix data structure, and recursively computed the cluster tree and populated it with admissible or inadmissible blocks, as described in Section 1.3.

In our code, we have implemented three operations of the shell matrix: matrix-vector multiplication, shift of origin A := A + σI, and extraction of the diagonal. These are simply calls to the corresponding HLIB functions, appropriately wrapped according to PETSc convention. In addition, we have also implemented a shell spectral transformation in SLEPc that similarly implements a specialized version of the shift-and-invert technique of (1.15) by means of HLIB's LU decomposition.

A final note about the implementation is that both the exponentialintegral and incomplete Gamma functions are available via GSL, the GNU Scientific Library.

In Table 4.1 we show the CPU time required for the matrix generation phase with dimension n varying from 4000 to 256000. With a uniform grid the resulting matrix is symmetric and the code takes this fact into consideration: the generation time for the symmetric case is almost half of the time required for the unsymmetric counterpart. The H-matrix approach, either with Taylor or Lagrange for computing the admissible blocks, represents a significant gain in generation time compared with the version with conventional sparse storage (note that in the sparse version we compute all matrix elements and then decide whether they are too small to be stored). The time reported for the SVD version includes the computation of matrix elements as in the sparse version as well as the time required for low-rank approximation through SVD decomposition (with LAPACK). Although this variant is the most expensive one, if the problem is to be solved several times one may consider this approach since it allows both for a fixed rank-k and for a rank satisfying a prescribed tolerance, as mentioned in Section 3.5. For large values of n the CPU time required to compute the entries is prohibitive for the sparse implementation. The slight differences reported for Taylor and Lagrange result from the fact that in the latter case we implemented numerical quadrature while for Taylor we used the formulae presented in Section 3.1;

the computation of the incomplete Gamma function at the required points is skewing the results a bit. For increasing problem size there is a constant growth factor less than three for these two approaches while the problem size is quadrupling. The growth factor respects the estimated n log(n) asymptotic cost, in contrast with the sparse version that follows n 2 . Some values were not reported due to their high value.

Table 4.2 complements the previous comments, showing the number of stored elements for all approaches. Note that the actual memory requirements for sparse storage are quite larger, since the space needed for indices is considerable, while for the H-matrix representation the overhead in negligible. The last column shows the compression factor for H-matrix format as a percentage of the full (dense) storage, revealing noteworthy gains for increasing values of n.

Table 4.3 reports on the CPU time for the solution phase using Taylor, Lagrange and SVD data-sparse representation as well as the sparse approach. values with relative tolerance on the residual of 10 -7 ), the shift-and-invert technique is required to enable convergence of the Krylov-Schur method. In the following, an LU factorization on the H-matrix representation is used in the linear solver required in the application of the shift-and-invert operator.

The factorization is the most costly operation but is performed only once, while triangular solves are required at each iteration of the eigensolver. In Table 4.3 we show the factorization time as well as the total solution time.

For these tests, we used a Krylov basis of dimension 16, and with this size the method does not need to restart (except for the matrix of n = 256000

where two restarts are required).

As expected, the computation of eigenpairs with the implementations of the data-sparse representation is very fast compared to the sparse storage, which shows a fast degradation in performance for increasing dimension. The SVD approach is competitive with Lagrange and Taylor approximations, and results for SVD on large values of n are not reported only due to the high generation time. As mentioned above, the present problem is hard to solve since for incresing values of n, and for fixed τ , the eigenvalues tend to become more and more clustered. For problems with better separation of the spectrum, the shift-and-invert step can be avoided and consequently its

Sensitivity analysis

We consider now a smaller problem, τ = 1000 and n = 2000, to assess the sensitivity of the proposed resolution method regarding to some of the parameters. The two main parameters under review are k, the degree on the Taylor or Lagrange approximations, and the minimum size for the blocks, named bound. The former is already well understood, larger values of k should improve the quality of the approximation but at a higher cost. The latter is a numerical way to stop the splitting process before reaching sets containing only one element; it defines which sets are small enough.

In Figure 4 For the problem under study, Lagrange approximations seem to be slightly better than Taylor ones, allowing best approximations either for fixed values of bound or degree. This can be justified by the fact that the former is a global approximation while the latter is of local nature. Furthermore, Lagrange is more flexible since it does not require the computation of derivatives. As stated previously, SVD decomposition allows for the best approximation but it is very expensive to compute. ) for Taylor and of O(10 -9 ) for Lagrange.

One should bear in mind that the required precision for the solution was 10 -7 , so both solutions should be considered correct under the precision of data. For the same bound but reducing the degree to 3, both approximations deliver a relative error of O(10 -7 ). For larger values of bound, the number of correct digits of the approximate solutions with respect to the reference case increases, so results are insensitive to this parameter. To explore further this aspect, one should produce more refined eigenvalue approximate solutions.

On the other hand, smaller values of bound, 10 or less, begin to produce solutions under the required precision, so the relative error, as mesured in this section, becomes larger.

For SVD, the only case where solutions don't match at least at O(10 -7 ) is for degree 9 and bound 5 (indeed, as explained, this combination of parameters is no sense). Except for the previous case, singular value decomposition captures extremely well the relevant information for the eingenvalue computation.

To better clarify the impact of the bound parameter on the H-matrix, we plot in Figures 4.2 to 4.5 the structure for a small size problem, τ = n = 100, and for a fixed degree equal to 6. It is clear, Figure 4.5, that a larger bound, is not adequate for such a small size problem. At the other extreme, a small bound, Figure 4.2, produces very small inadmissible (and admissible) blocks, giving rise to important information, green shapes, in blocks very far from the diagonal; for such a small bound, maybe a higher degree would be recommended. Part III

Conclusions and future research

Chapter 5

Conclusions

The H-matrix technique has its origin in the panel clustering [START_REF] Hackbusch | On the fast matrix multiplication in the boundary element method by panel clustering[END_REF][START_REF] Hackbusch | The panel clustering algorithm[END_REF] where the kernel of the operator is approximated by a degenerate one. However, the panel clustering matrices cannot cheaply be multiplied or inverted, and we arrive to one great advantage of using H-matrices: the arithmetic, that is, matrix-vector multiplication, matrix-matrix multiplication, inversion, matrix functions, can be performed in almost linear complexity. So, new applications as the solution of matrix equations or just the construction of very efficient preconditioners for linear systems can be handled with efficacy using this technique, [START_REF] Börm | Construction of data sparse H 2 -matrices by hierarchical compression[END_REF].

In this work we solve a large eigenvalue problem issued from a radiative transfer equation in stellar atmospheres. The computation is time intensive both in memory requirements, since large dimensional cases are to be treated, and in the solution phase, due to the clustering of the eigenvalues.

Sophisticated numerical algorithm implementations for the solution phase are thus required. The H-matrix data representation provided by HLIB was integrated in the SLEPc and PETSc frameworks to tackle those two difficulties.

We report on the numerical low-rank approximations developed, on the details of the integration of the libraries under use, as well as presented CHAPTER 5. CONCLUSIONS a brief explanation of the problem and its interesting characteristics. The combined use of numerical solution methods optimized for high-end computer environments and of the clever data storage representations, although not with ease, provide an efficient and fast answer capable of solving even larger dimensional problems, maybe requiring the use of distributed processing.

The numerical tests illustrate the success of the proposed solutions.

Moreover, error bounds for the proposed low-rank approximations were developed, ensuring that under the stated conditions, these approximations are mathematically supported.

This thesis also presents a useful summary on numerical state-of-art methods for the solution of algebraic eigenvalue problems and on some of today's most important software libraries for high-end scientific computing.

Intentionally, some important computational skills to fullfil the aims of this work were hidden to keep the reading as much clear and fluent as possible.

Indeed, dealing with such huge number of software libraries, knowing in depth the merits and limitations of the numerical algorithms involved, requires a high level of specialization. Furthermore, being familiar with computer languages is crucial.

Chapter 6

Future research

In this chapter, we make some comments on future research directions that may amplify the application of the H-matrix technique to other cases, which would provide the next steps along the path to a more wide application of it.

Moreover, in Section 6.2 it is described some state-of-art about some topics in nonlinear eigenproblems.

Some limitations on H-matrix technique

We have started out dealing with an integral operator defined from L 1 ([0, 1]) into itself, and with kernel

f (τ, σ) = ln (1 -cos 2π (τ -σ)) , (6.1)
aiming to compute some few eigenvalues and eigenvectors of the discretized matrix of the operator, with an analogous aproach as it was done with the operator T in (1.2) with kernel (1.3). The kernel (6.1) was purposely chosen as it is a wealky singular kernel for which the exact eigenvalues are known:

• simple eigenvalue: λ 0 = -ln 2,

• double eigenvalues:

λ k = -1 k , k ≥ 1.

CHAPTER 6. FUTURE RESEARCH

It would be an interesting work to compute a few eigenvalues using the technique used in the thesis, followed by the comparsion between them and the corresponding exact ones.

Nevertheless, the study done with this kernel showed us the apparent restriction of the use of the H-matrix technique, similiar to what was done with the kernel (1.3) in this thesis. We have encountered some difficulties: we did not manage to find a solution to the computation of the double integral of the kernel (6.1), necessary to achieve the formula for the inadmissible leaves of the H-matrix structure; we were unable to find a recursive rule for the derivative of the kernel, essential to the construction of the Taylor series expansion of the kernel, one of the options for the low-rank construction; we did not manage to write the kernel with both variables separated, prerequisite to construct a degenerate approximation of the kernel, see (1.17), for the lowrank approximation.

Despite the previous facts, polynomial interpolation is an optimal option for the low-rank approximation in this case, using a similar procedure to what was done in Section 3.3, except for the cases of low-rank approximation close to singularities, besides that the question remains how to compute the inadmissible leaves. These two points appear to be a good point of research in an immediate future work, and we leave here a hint of the study we are having in hands.

One possibility is to approximate the new operator T by Tn obtained from T by means of a particular quadrature approximation. First we rewrite the kernel as a function of one variable, f (τ, σ) = f (r), with r = |τ -σ|.

The essence of this procedure is using a truncated kernel function fn instead of the kernel itself. The truncation is done close enough to the singularity, that is, given a point c n of the same order of the weights of the quadrature formula (e.g. for trapezoidal quadrature rule it is of order 1 n ), for the points of the domain located between c n and the singularity, the kernel is evaluated as fn (c n ); in the other points of the domain, fn is defined to be f . These ideas may be found in [START_REF] Largillier | A note on a collectively compact approximation for weakly singular integral operators[END_REF].

FIRST GLANCE ON NONLINEAR EIGENVALUE PROBLEMS 91

6.2 First glance on nonlinear eigenvalue problems At the begining of the research project, we had thought to make some contributions in nonlinear eigenvalue problems, but in fact the research followed other directions. Nevertheless, some investigation about this subject was done, speaking in terms of state-of-the-art, and in this section we present a summary of it as we consider nonlinear eigenproblems one topic for research in the future.

Regarding the importance of spectral theory in the case of linear operators, for instance, spectra of differential operators applied to the theory of elliptic boundary value problems or the application in classical quantum mechanics, it is quite natural that several efforts have been made to study and define spectra in the nonlinear case.

In order to recall the definition of spectrum in the linear case, let X and Y be two Banach spaces over K = R or K = C and an operator L ∈ BL(X), where BL(X) is the algebra of all bounded linear operators from X into itself. The resolvent set of L is defined by re (L) = {λ ∈ K : λI -L is a bijection} (6.2) and the spectrum of the operator by

sp (L) = K\re(L) = {λ ∈ K : λI -L is not a bijection} , (6.3)
where I is the identity operator.

The study of nonlinear eigenvalues is older than the nonlinear spectral theory, remarking that initially, the word spectrum was used in the sense of the set of eigenvalues, the point spectrum. Considering F : X -→ X a continuous nonlinear operator, we note that the definition of the point spectrum of F coincides with the one of the linear case sp p (F ) = {λ ∈ K : F (x) = λx for some x = 0} . (6.4) CHAPTER 6. FUTURE RESEARCH It was just in the late 60's that Kachurovskij and Neuberger independently studied spectral values that were not necessarily eigenvalues. Then, serveral different spectral "theories" started to appear, aiming to preserve the properties considered useful regarding the linear case, as well as trying to extend them to a possible wider range of nonlinear problems.

In [START_REF] Appell | Nonlinear spectral theory[END_REF], the authors describe that it was assumed that a spectrum of a continuous nonlinear operator should satisfy some essential requirements, as reducing to the familiar spectrum in case of linear operators; keeping some of the usual properties of the linear spectrum (e.g. compactness, nonemptiness), containing the eigenvalues of the operator considered, as well as having applications (e.g. existence, uniqueness).

On the contrary of what could be possibly expected, nonlinear spectral theory have faced some "disappointing" facts. For example, in oposition to the linear case: the spectrum of a nonlinear operator has almost no information about the operator itself; in general the properties of boundedness, closedness, or nonemptiness fail; the spectra can be disjoint from the set of eigenvalues.

Next we give a glimpse over some spectra that have appeared over the last forty years:

• Rhodius spectrum [Rho77],
Consider F : X -→ X a continuous operator such that F (0) = 0. The is the Kachurovskij spectrum of F .

• Dörfner spectrum [Dör97],

Consider F : X -→ X linearly bounded on X, i.e.

[F ] B = sup Finally, F is said to be FMV-regular if F is stably solvable, [F ] α > 0 and [F ] q > 0. Now we have all the notations needed to define the Furi-Martelli-Vignoli spectrum (FMV spectrum):

x =0 F (x) x < ∞.
sp F M V (F ) = K\ {λ ∈ K : λI -F is F M V -regular} . • Feng spectrum [Fen97],
Consider F : X -→ Y continuous on X. We give the following notations to achive the definition of the Feng spectrum: x in a constricted neighborhood of p. Moreover, F is said to be y-epi at p if it is y-admissible at p and y-epi on any small enough neighborhood of p (the definition of F y-epi can be seen, e.g. in [START_REF] Calamai | An overview on spectral theory for nonlinear operators[END_REF]). Considering F, J : X -→ Y continuous operators, the point spectrum (6.4) can be written in a more general way sp p (F, J) = {λ ∈ K : F (x) = λJ(x) for some x = 0} . (6.8)

FMV

Anyway, the first definition is a particular case of (6.8), that is, sp p (F ) = sp p (F, I), with X = Y .

Recalling that, in the linear case, the approximate point spectrum of L is

CHAPTER 6. FUTURE RESEARCH

As was mentioned before, the FMV-spectrum has asymptotic characteristics, and so (6.8) should be replaced by the asymptotic point spectrum sp q (F, J) = {λ ∈ K : [λJ -F ] q = 0} , that is always in sp F M V (F ).

In [START_REF] Santucci | On the definition of eigenvalues for nonlinear operators[END_REF][START_REF] Appell | Nonlinear spectral theory[END_REF] it is possible to find point phantom of (F, J) defined as φ p (F, J) = {λ ∈ K : λ connected eigenvalue of(F, J)}

where by connected eigenvalue of (F, J) is meant the nullset N (λJ -F ) = {x ∈ X : F (x) = λJ(x)} containing an unbounded connected set C with 0 ∈ C.

In [START_REF] Calamai | A new spectrum for nonlinear operators in Banach spaces[END_REF] it is defined the approximate point spectrum of F at p, a subset of sp CF V (F ) sp CF V ap (F, p) = λ ∈ K : [F ] A p = 0 and [F ] q p = 0 that coincides with the usual approximate point spectrum in the linear case.

Combining some interesting summary schemes in [App03, APV00, ADPV04],

relating some possible inclusions between the spectra, we present also the Table 6.2 as a possible summary, for F, J : X -→ Y nonlinear operators and Table 6.2: Possible inclusions between some nonlinear spectra φ(F, J) ⊆ sp F M V (F, J) ⊆ sp F (F, J) ⊆ ⊆ ⊆ φ p (F, J) ⊆ sp q (F, J) sp p (F, J)

when reducing to the linear case, L ∈ BL(X) with J = I, we have the Table 6.3.

What seems it should be done, when we want to apply spectral theory to a particular nonlinear problem, is to choose with care a spectrum having at least some of the needed characteristics. In [START_REF] Appell | Nonlinear spectral theory[END_REF] the authors state that it is not the intrinsic structure of the spectrum itself that results in interesting applications, but considered a tool for solving nonlinear equations. It is also referred that in nonlinear analysis, two techniques have been effectively used for the study of nonlinear eigenproblems: topological methods (fixed point theorems, degree theory) and variational methods (critical points, nonlinear functionals).

The knowledge of the spectrum of a nonlinear operator seems to be useful in solvability results for nonlinear equations, moreover we emphasise that bifurcation theory, one of the most important fields of nonlinear functional analysis, is closely related to nonlinear eigenvalue problems. Other more applications, as e.g. to boundary value problems, may be found in literature.

Nonlinear eigenvalue problems arise in many applications, e.g. structural dynamics, acoustics, fluid mechanics, control theory and quantum physics.

Usually the procedure is to discretize the problem (e.g. Galerkin schemes) and to the resulting matrix problem an iterative projection method may be applied .

It can also be held what is called linearization, for example, in polynomial or rational eigenvalue problems in which first it is performed a linearization, resulting in a larger linear eigenvalue problem with the same eigenvalues and then apply a method for linear eigenvalue problems. However, in [START_REF] Mehrmann | Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods[END_REF] the authors state that this is not ideal, since it makes the problem larger and moreover may considerably increase the conditioning of the problem. Even though, this approach seems to be appropriate for this kind of problems.

  Chave com classificação MAS -2010 Mathematics Subject Classification, (MSC[2010]): 33 Special functions; 33F Computational aspects; 33F05 Numerical approximation and evaluation; 45 Integral Equations; 45C Eigenvalue problems; 45C05 Eigenvalue problems; 45E Singular integral equations; 65 Numerical analysis; 65F Numerical linear algebra; 65F15 Eigenvalues, eigenvectors; 65Y Computer aspects of numerical algorithms; 65Y20 Complexity and performance of numerical algorithms.

x 1 . 1 1

 111 tive avec ces matrices. Nous incorporons certaines routines comme la multiplication matrice-vecteur ou la décomposition LU, en SLEPc (Hierarchical matrices LIBrary) pour explorer les algorithmes existants afin de résoudre les problèmes de valeur propre. Nous développons aussi des expressions analytiques pour l'approximation des noyaux dégénérés utilisés dans la thèse et déduire ainsi les limites supérieures d'erreur pour ces approximations. Les résultats numériques obtenus avec d'autres techniques pour résoudre le problème en question sont utilisés pour la comparaison avec ceux obtenus avec la nouvelle technique, illustrant l'efficacité de ce dernier. The problem and its discretization . . . . . . . . . . . . . . . 1.2 Eigensolvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Iterative eigensolvers . . . . . . . . . . . . . . . . . . . 1.2.2 Computation of eigenvalues around a given target . . . 1.3 H-matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Numerical libraries 2.1 HLIB library . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 PETSc and SLEPc libraries . . . . . . . . . . . . . . . . . . . 2.3 BLAS and LAPACK libraries . . . . . . . . . . . . . . . . . . 2.4 OpenMP library . . . . . . . . . . . . . . . . . . . . . . . . . . II Development of numerical approximations 3 Low-rank approximations 3.1 Taylor approximation . . . . . . . . . . . . . . . . . . . . . . . 3.2 Error bounds for Taylor approximation . . . . . . . . . . . . . 3.3 Interpolation approximation . . . . . . . . . . . . . . . . . . . xi xii CONTENTS 3.4 Error bounds for interpolation approximation . . . . . . . . . 3.5 SVD approximation . . . . . . . . . . . . . . . . . . . . . . . . 4 Numerical results 4.1 Hardware . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Serial approach . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Parallel approach (threads) . . . . . . . . . . . . . . . . . . . . 4.4 Sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . . Some limitations on H-matrix technique . . . . . . . . . . . . 6.2 First glance on nonlinear eigenvalue problems . . . . . . . . . References A Research projects List of Figures 1.1 Shift-and-invert transformation . . . . . . . . . . . . . . . . . 1.2 Example of a T I . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Example of a T I×I . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Example of an H-matrix . . . . . . . . . . . . . . . . . . . . . 4.1 Norm of the difference between matrix A n and its H-matrix representation, for several values of the degree k and minimum bound. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 H-matrix structure for Lagrange approximation with degree 6 and bound 5 (τ = 100, n = 100). . . . . . . . . . . . . . . . 4.3 H-matrix structure for Lagrange approximation with degree 6 and bound 10 (τ = 100, n = 100). . . . . . . . . . . . . . . 4.4 H-matrix structure for Lagrange approximation with degree 6 and bound 20 (τ = 100, n = 100). . . . . . . . . . . . . . . 4.5 H-matrix structure for Lagrange approximation with degree 6 and bound 40 (τ = 100, n = 100). . . . . . . . . . . . . . . xiii xiv LIST OF FIGURES List of Tables 4.1 CPU time for the generation phase . . . . . . . . . . . . . . . 4.2 Number of stored elements in the case of dense, sparse and H-matrix representations . . . . . . . . . . . . . . . . . . . . . 4.3 CPU time for the solution phase . . . . . . . . . . . . . . . . . 4.4 Computed eigenvalues . . . . . . . . . . . . . . . . . . . . . . 4.5 Execution time for the matrix generation in parallel computation 4.6 List of the 5 largest eigenvalues produced by the sparse version and Taylor and Lagrange approximations, using degree 6, bound 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 A general view of the main properties of some nonlinear spectra 6.2 Possible inclusions between some nonlinear spectra . . . . . . 6.3 Possible inclusions of different spectra applied to the linear case xv xvi

A

  coefficient matrix of the system of linear equations A i,j matrix entry in the row i and column j A T transpose of matrix A A -1 inverse of matrix A H-matrix hierarchical matrix I index set of the basis function operator L sp R (F ) Rhodius spectrum of the operator F xvii sp N (F ) Neuberger spectrum of the operator F sp K (F ) Kachurovskij spectrum of the operator F sp D (F ) Dörfner spectrum of the operator F sp F M V (F ) Furi-Martelli-Vignoli spectrum of the operator F

  being exhaustive, eigenvalues and eigenvectors are used in the study of nuclear reactor dynamics (stability of neutron fluxes [...]), in finite element dynamic analysis of structural models (e.g., seismic simulations of civil infrastructure [...]), in the design of the next generation of particle accelerators [...], in the definition of a set of eigenfaces in biometric-based identification systems [...], in the solution of Schrödinger's equation in chemistry and physics [...], in the design of microelectromechanical systems (MEMS [...]), and in the study of conformational changes of proteins [...]. Because of the need for higher levels of simulation details and accuracy, the size and complexity of the computations grow as fast as the advancement of the computer hardware. In order to cope with the increasing need of solving eigenvalue problems, various useful numerical algorithms that are suitable for solving large-scale eigenvalue problems are developed."

  software libraries used. These chapters constitute the first part of the thesis, designated by Framework and guidelines. The development of the numerical approximation is done in Part II; in Chapter 3, three possible strategies for obtaining a low-rank approximation are discussed, using H-matrix representation, while on Chapter 4 results of some numerical experiments along with specific details of the implementation are shown. The last two chapters, which embodies Part III, are devoted to conclusions and possible extensions for future research.

Figure 1 . 1 :

 11 Figure 1.1: Shift-and-invert transformation

Figure 1 . 3 :Definition 1 .

 131 Figure 1.3: Example of a T I×I based upon T I of Figure 1.2.

Figure 1 . 4 :

 14 Figure 1.4: Example of an H-matrix structure based upon T I×I of Figure 1.3.

  SLEPc, the Scalable Library for Eigenvalue Problem Computations[START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF][START_REF] Hernandez | SLEPc users manual[END_REF], is a software package for the solution of large-scale eigenvalue problems on parallel computers. It can be used to solve standard and generalized eigenvalue problems, as well as other types of related problems such as the quadratic eigenvalue problem or the singular value decomposition. SLEPc can work with either real or complex arithmetic, in single or double precision, and it is not restricted to symmetric (hermitian) problems. It can be used from code written in C, C++, and FORTRAN. Most of the methods in this library are projection methods, including different variants of Krylov and Davidson iterations.SLEPc is able to cope with different problem types such as hermitian, non-hermitian (the default case), Generalized hermitian, Generalized nonhermitian, and Generalized non-hermitian with positive (semi-)definite B (a generalized eigenproblem is commonly written as Ax = λBx). The library is very flexible and it is possible to specify how many eigenvalues/eigenvectors to compute. In relation to the eigenvalues of interest, it is possible, in real symmetric problems, to compute the largest or smallest eigenvalues in magnitude, the leftmost or rightmost ones and even those closest to a given target value. For the complex case, more options are available. SLEPc provides a collection of eigensolvers: Power Iteration with deflation, Subspace Iteration with Rayleigh-Ritz projection, Arnoldi method, Lanczos method, Krylov-Schur and Davidson-type solvers. Most of the eigensolvers are based on the subspace projection paradigm, in particular, it includes a robust and efficient parallel implementation of Krylov-Schur method described in Section 1.2, which is the default solver in SLEPc. Sev-CHAPTER 2. NUMERICAL LIBRARIES eral Davidson-type solvers are included as well, in particular Generalized Davidson and Jacobi-Davidson, with various possibilities for the computation of the correction vector. In these solvers, the user can easily select which preconditioner to use. Apart from eigensolvers, some spectral transformations such as the shiftand-invert technique of (1.15) are available, where the user can compute interior eigenvalues with the help of linear solvers included in PETSc (Portable, Extensible Toolkit for Scientific Computation).

Figure 4 . 1 :

 41 Figure 4.1: Norm of the difference between matrix A n and its H-matrix representation, for several values of the degree k and minimum bound.

  Figure 4.4 illustrates how this type of approximation can be used as a clever and cheap way to treat the problem at a cost of a block-band matrix eigenvalue computation. Probably the best choice would be Figure 4.3, combining a moderate size for the inadmissible blocks and setting the admissible blocks neither in an atomized form nor too much aggregated.

Figure 4 . 2 :

 42 Figure 4.2: H-matrix structure for Lagrange approximation with degree 6 and bound 5 (τ = 100, n = 100).

  set sp R (F ) = K\ {λ ∈ K : λI -F homeomorphism on X} = K\ λ ∈ K : λI -F is a bijection and (λI -F ) -1 continuous is the Rhodius spectrum of F . • Neuberger spectrum [Neu69],Consider F : X -→ X a continuously differentiable operator such thatF (0) = 0. The set sp N (F ) = K\ {λ ∈ K : λI -F diffeomorphism on X} = K\ λ ∈ K : λI -F is a bijection and (λI -F ) -1 of class C 1 6.2. FIRST GLANCE ON NONLINEAR EIGENVALUE PROBLEMS 93 is the Neuberger spectrum of F . • Kachurovskij spectrum [Kac69],Consider F : X -→ X belonging to the class Lip(X) such thatF (0) = 0. F ∈ Lip(X) means that F is Lipschitz continuous on X, that is, [F ] Lip = sup x =y F (x) -F (y) x -y < ∞.The set sp K (F ) = K\ {λ ∈ K : λI -F lipeomorphism on X} = K\ λ ∈ K : λI -F is a bijection and (λI -F ) -1 ∈ Lip(X)

  The set sp D (F ) =K\ λ ∈ K : λI -F is a bijection and (λI -F ) -1 is linearly bounded on X} is the Dörfner spectrum of F .• Furi-Martelli-Vignoli spectrum [FMV78],Before the definition of this spectrum, some other definitions have to be presented. The (Karatowski) measure of compactness α(A) of a bounded set A ⊂ X is defined as the infimum of real numbers δ > 0 such that A admites a finite covering by sets of diameter less than δ. Considering F : X -→ Y continuous, the following numbers are defined [F ] α = sup {k > 0 : α (F (A)) ≥ kα (A) for every bounded A ⊂ X} , operator F is said to be stably solvable if, given any compact operator G in X with [G] Q = 0, the equation F (x) = G(x) has a solution x ∈ X.

[

  F ] A = inf {k > 0 : α (F (A)) ≤ kα (A) for every bounded A ⊂ X} ,Considering G : B r -→ Y a continuous operator, with B r = {x ∈ X : x ≤ r}, it is defined ν r (F ) = inf {k > 0, there exists G satisfying [G] A ≤ 0, G | ∂Br ≡ 0, and F (x) = G(x), ∀x ∈ B r } and ν(F ) = inf {ν r (F ), r > 0} .

Finally,F•

  is said to be F-regular if F is stably solvable, [F ] α > 0 and[F ] b > 0 and ν(F ) > 0. Now we have all the notations needed to define the Feng spectrum:sp F (F ) = K\ {λ ∈ K : λI -F is F -regular} .6.2. FIRST GLANCE ON NONLINEAR EIGENVALUE PROBLEMS 95 Väth phantom [SV00],Let U ⊂ X be open, bounded, connected and with 0 ∈ U . ConsiderF : U -→ Y continuous. F is epi on U if F (x) = 0 on ∂U and the equation F (x) = G(x) has a solution in U for any compact operator G : U -→ Y that satisfies G(x) ≡ 0 on ∂U . F is strictly epi on U if inf { F (x) : x ∈ ∂U } > 0,and there exists some k > 0 such that for any operatorG : U -→ Y that satisfies G(x) ≡ 0 on ∂U and [G] A ≤ k, the equation F (x) = G(x)has a solution in U .Finally, F : X -→ Y is said to be v-regular if it is strictly epi on some U . Now we have all the notations needed to define the Väth phantom:φ(F ) = {λ ∈ K : λI -F is not v-regular} . • Calamai-Furi-Vignoli spectrum [CFV10], Consider U ⊂ X open, the continuous operator F : U -→ Y and p ∈ U . Let U p denote the open neighborhood {x ∈ X : p + x ∈ U } of 0 ∈ X and consider F p : U p -→ Y continuous with F p (x) = F (p + x) -F (p).Recalling (6.7) and (6.6), considering B(p, r) ⊂ X an open ball with center at p and radius r > 0, B(p, r) ⊂ U , it is defined now[F ] A p = lim r→0 [F | B(p,r) ] A , [F ] q p = lim r→0 [F | B(p,r) ] q ,and, for y ∈ Y , F is y-admissible at p if F (p) = y and F (x) = y for every

  defined by sp ap (L) = {λ ∈ K : exists a sequence x n , x n → 0 such that λx n -Lx n → 0 as n → ∞} in some of the nonlinear spectra a corresponding notion of eigenvalues (point spectrum or approximate point spectrum) was introduced; note that sp p (L) ⊆ sp ap (L) and sp(L) = sp ap (L) ∪ sp p (L * ) = sp p (L) ∪ sp ap (L * ), being L * ∈ BL(X * ) the adjoint of L.

6. 2 .

 2 FIRST GLANCE ON NONLINEAR EIGENVALUE PROBLEMS 99 Table 6.3: Possible inclusions of different spectra applied to the linear case sp(L) = φ(L, I) = sp F M V (L, I) = sp F (L, I) ⊆ ⊆ ⊆ φ p (L, I) ⊆ sp q (L, I) ⊇ sp p (L, I)

Table 4 .

 4 1: CPU time (in seconds) for the generation phase for τ = 4000 and

	varying n, with a uniform (left) and non-uniform (right) grid. The results
	correspond to bound=80 and degree=6.
			Uniform			Non-uniform
	n Taylor Lagrange SVD Sparse Taylor Lagrange SVD Sparse
	4000	2.1	2.0	20.1	15.3	3.9	3.5	38.6	30.9
	8000	6.1	5.7	99.6	61.8	10.4	9.4	182.7	123.3
	16000	16.7	15.6	536.0	245.4	29.5	26.9	927.2	496.4
	32000	49.6	47.0	3505.0 982.1	85.4	78.9	5630.7 1972.0
	64000 140.6	133.5	-	-	245.4	227.2	-	-
	128000 394.4	374.4	-	-	690.3	637.9	-	-
	256000 1019.8	958.3	-	-	1783.4	1616.1	-	-

Table 4 .

 4 2: Number of stored elements in the case of dense, sparse and Hmatrix representation, for the non-uniform grid case in Table4.1.

	n	Dense	Sparse H-matrix Compression
	4000 1.6 • 10 7 2.9 • 10 5 1.9 • 10 6	12.0%
	8000 6.4 • 10 7 1.2 • 10 6 4.2 • 10 6	6.5%
	16000 2.6 • 10 8 4.5 • 10 6 9.1 • 10 6	3.5%
	32000 1.0 • 10 9 1.8 • 10 7 1.9 • 10 7	1.9%
	64000 4.1 • 10 9	-	4.1 • 10 7	1.0%
	128000 1.6 • 10 10	-	8.6 • 10 7	0.5%
	256000 6.6 • 10 10	-	1.9 • 10 8	0.28%
	Since the spectrum is tightly clustered (see Table 4.4 for the five largest eigen-

Table 4

 4

	.3: CPU time (in seconds) for the solution phase for τ = 4000
	and varying n, with a non-uniform grid (non-symmetric case). The results
	correspond to bound=80 and degree=6.
			Solution			Factorization
	n Taylor Lagrange SVD Sparse Taylor Lagrange
	4000	0.3	0.3	0.3	0.2	0.3	0.3
	8000	0.5	0.5	0.5	1.1	0.4	0.4
	16000	1.0	1.1	1.2	8.2	0.8	0.9
	32000	2.6	2.7	2.9	66.3	2.1	2.3
	64000	6.5	6.8	-	-	5.6	5.9
	128000	16.4	17.3	-	-	14.5	15.4
	256000	47.5	49.3	-	-	39.3	41.3
	Table 4.4: Computed eigenvalues for the case of a uniform grid with n =
	16000 and τ = 4000.
	Eigenvalue		Taylor	Lagrange		SVD
	λ 1	0.749999843422 0.749999843459 0.749999843598
	λ 2	0.749999374216 0.749999374253 0.749999374391
	λ 3	0.749998592208 0.749998592245 0.749998592383
	λ 4	0.749997497401 0.749997497438 0.749997497576
	λ 5	0.749996089801 0.749996089838 0.749996089976
	computational cost.
	As mentioned at the end of Section 1.2, Davidson methods do not seem
	appropriate in the context of hierarchical matrices, since the shift-and-invert
	technique is very cheap in this case. However, we wanted to do some ex-
	periments. With Jacobi (diagonal) preconditioning, we were able to solve
	the problem (although after many iterations) by tuning the parameters of
	SLEPc's Davidson solver. For instance, for n = 8000 with uniform grid, the
	response time is 15.5 seconds, as opposed to 0.5 seconds with shift-and-invert
	Krylov-Schur. A much powerful preconditioner is to use the LU factoriza-

  .1 we show the distance from the H-matrix approximation to the original one, A n . For most part of the combinations of both parameters, List of the 5 largest eigenvalues produced by the sparse version and Taylor and Lagrange approximations, using degree 6, bound 20 the distance between the matrices is sufficiently small, providing good perspectives for the success of the numerical approximations. It is important to note that the eigenvalues are computed with a tolerance of 10 -7 . Smaller values for bound reduce the minimum sizes of the blocks thus increasing the number of blocks, particularly the low-rank approximation ones, and the Hmatrix becomes sparser. With the decrease of bound the admissibility condition becomes eager. Similarly, lower values of degree may produce worst but

	Eigenvalue	Sparse	Taylor	Lagrange
	λ 1	0.749997399534 0.749997056849 0.749997402836
	λ 2	0.749989598317 0.749989290172 0.749989601309
	λ 3	0.749976596889 0.749976284686 0.749976599928
	λ 4	0.749958396151 0.749958079593 0.749958399195
	λ 5	0.749934997365 0.749934670108 0.749935000520
	Table 4.6: cheaper approximations. These two parameters must be considered together,
	since it may not make sense, for instance, to consider a degree greater than a
	bound. Very low values of either parameters may lead to wrong answers (val-
	ues reported in red and orange color show a significative difference between
	A n and its H-matrix approximation). Moreover, this analysis is, naturally,
	problem dependent, although there are geometrical considerations of general
	nature.

Table 4 .

 4 6 shows the computed eigenvalues for the Sparse, Taylor and Lagrange approximations, for degree 6 and bound 20. The relative error (distance from the approximation solution to the one provided by the Sparse version (reference)) is of O(10 -7

these entries are given by(1.7)

variables this doesn't occur, having private copies made for each thread in the thread team, and in this way, changes made in one thread are not visible in the private variables in the other threads.

Finally, OpenMP provides multiple types of synchronization to help in many different situations (synchronization when multiple threads are running at the same time and synchroniztion between them is needed).

Part II

Development of numerical approximations Chapter 4 Numerical results

In this chapter we begin by a concise summary of the machine that has been used and present some numerical experiments that aim at illustrating the benefits of the H-matrix representation with respect to conventional storage, both in terms of performance and memory requirements.

Hardware

The tests in Sections 4.2 and 4.3 have been executed on a Linux workstation with an Intel Core i7 950 processor at 3,06 GHz with 8 MB of L3 cache memory and 8 GB of main memory. This processor has 4 cores with hyperthreading technology (a total of 8 virtual processors).

The less computer intensive tests, in Section 4.4 , have been performed on a Linux Ubuntu 11.04 desktop edition with an Intel Core i5 M450 at 2.40 GHz with 4GB of main memmory (4 cores).

Serial approach

For all the tests, in Sections 4.2 and 4.3, we chose to use a fixed value of the τ parameter, in particular τ = 4000. We also set a constant value for CHAPTER 4. NUMERICAL RESULTS

Parallel approach (threads)

We have also developed a straightforward parallel version of the H-matrix generation, based on the OpenMP API for shared-memory parallel programming. In particular, we follow a tasking approach with the OpenMP task directive [ACD + 09], where each recursive call constitutes a new task.

Regarding parallelization of the generation phase, the multi-threaded version was analyzed up to 8 threads.

Table 4.5 shows the measured execution times along with the achieved speedups, for the two longest times in Table 4.1 (non-uniform grid, Taylor with n = 256000 and SVD with n = 32000). In the case of the Taylor approximation, speedup is virtually ideal up to 4 threads and decays significantly later. This can be attributed in part to the fact that only 4 physical cores are available. In the case of SVD generation, speedup is much worse, thus revealing a problem with load imbalance, due to the fact that the high cost of

SENSITIVITY ANALYSIS CHAPTER 6. FUTURE RESEARCH

Finally, F is said to be regular at p if [F ] A p > 0 and [F ] q p > 0 and F p is 0-epi at 0. Now we are able to define the spectrum of the map F at the point p, that is, the Calamai-Furi-Vignoli spectrum:

Besides the previous references to the papers where the respective spectrum was introduced, in [ADPV04, APV00, App03] just to mention a few, it is possible to find not only the definitions, but also the comparison between some of the spectra previously mentioned. Next, not aiming to be exhaustives, we illustrate some differences between these spectra.

The Rhodius, Neuberg, Kachurovskij and Dörfner spectra contain the eigenvalues of the operator, and 0 if the operator is compact in a Banach space of infinite dimension (this is analogous to the linear case). It can be seen that in these spectra, what was done in their definitions was simply the replacement of BL(X) in (6.2) and (6.3) by other classes of continuous nonlinear operators, in this four cases, respectively by the class of continuous, continuously differentiable, Lipschitz continuous, and linearly bounded operators.

Apparently, the asymptotic spectrum ( x → ∞) introduced by Furi, Martelli and Vignoli, the global Feng spectrum (x ∈ X), and the local spectrum x ∈ U given by Väth are more present in literature, probably because, even though being quite different in nature, are identical for homogeneous operators and when applied to linear operators they give exactly the usual spectrum of the linear theory. Note that usually the word "spectrum" in the case of Väth is replaced by "phantom" because its construction is quite far from what is usually called spectrum.

The fact that FMV spectrum do not contain the set of eigenvalues, has induced Feng to define a new spectrum having the same topological properties of the FMV spectrum, but containing the classical eigenvalues, as happens in the linear case. The Väth phantom is closed and can be in some cases bounded, hence compact.

FIRST GLANCE ON NONLINEAR EIGENVALUE PROBLEMS 97

In Table 6.1, inspired in [START_REF] Calamai | An overview on spectral theory for nonlinear operators[END_REF], we give a summary of the comparison of the different spectra described before, from the viewpoint of topological properties (having in mind the properties of the linear spectra).

Table 6.1: A general view of the main properties of some nonlinear spectra Spectrum nonempty closed bounded compact

CHAPTER 6. FUTURE RESEARCH

There is a vast literature on numerical methods for nonlinear eigenvalue problems. Generally, they are split into dense and large sparse problems. For the last ones, projection methods that work directly for general nonlinear eigenproblems, are a good option. In this case the search spaces have to be expanded by directions that have a high approximation potential for the desired eigenvector.

We have seen in Section 1.2 that for sparse linear eigenproblems iterative projection methods, like e.g. Arnoldi or Jacobi-Davidson, are a tool for the computation of a few eigenpairs of the spectrum. As previously observed, Jacobi-Davidson method for linear problems was proposed by Sleijpen and van der Vorst; the method was later extended to quadratic eigenproblems [START_REF] Betcke | A Jacobi-Davidson-type projection method for nonlinear eigenvalue problems[END_REF]. Voss considered also the use of this method for the nonlinear eigenproblem (nonsymmetric case) [START_REF] Voss | A Jacobi-Davidson method for nonlinear and nonsymmetric eigenproblems[END_REF].

Generally two approaches to subspace expansion may be found in the literature: Jacobi-Davidson [START_REF] Betcke | A Jacobi-Davidson-type projection method for nonlinear eigenvalue problems[END_REF][START_REF] Voss | A Jacobi-Davidson method for nonlinear eigenproblems[END_REF] and nonlinear Arnoldi [START_REF] Voss | An Arnoldi method for nonlinear eigenvalue problems[END_REF] type expansion. Both methods approximate inverse iteration, which is known to provide a direction of a high approximating potential to the eigenpair it is aiming for.

A more detailed study of these issues may be found in [START_REF] Mehrmann | Nonlinear eigenvalue problems: A challenge for modern eigenvalue methods[END_REF][START_REF] Voss | Numerical methods for sparse nonlinear eigenvalue problems[END_REF].

Appendix A

Research projects

This thesis is a cotutelle project developed within two Universities: the Uni-