We present a fast direct algorithm for the solution of linear systems arising from elliptic equations. We extend the work of Xia et al. (2009) on combining the multifrontal method with hierarchical matrices. We offer a more geometric interpretation of that approach, extend it in two dimensions to the unstructured mesh case, and detail an adaptive decomposition procedure for selectively refined meshes. Linear time complexity is shown for a quasi-uniform grid and demonstrated via numerical results for the adaptive algorithm. We also provide an extension to three dimensions with proven linear complexity but a more practical variant with slightly worse scaling is also described. / text

H-matrix techniques use a data-sparse tree structure to represent a dense or a sparse matrix. The leaves of the tree store matrix sub-blocks that are represented in full-matrix format or Rk-matrix (low rank matrix) format. H-matrix arithmetic is defined over the H-matrix representation, which includes operations such as addition, multiplication, inversion, and LU factorization. These H-matrix operations approximate results with almost optimal computational complexity. Based on the properties of H-matrices, the H-matrix preconditioner technique has been introduced. It uses H-matrix operations to construct preconditioners, which are used in iterative methods to speed up the solution of large systems of linear equations (Ax = b). To apply the H-matrix preconditioner technique, the first step is to represent a problem in H-matrix format. The approaches to construct an H-matrix can be divided into two categories: geometric approaches and algebraic approaches. In this thesis, we present our contributions to algebraic H-matrix construction approaches and H-matrix preconditioner technique. We have developed a new algebraic H-matrix construction approach based on matrix graphs and multilevel graph clustering approaches. Based on the new construction approach, we have also developed a scheme to build algebraic H-matrix preconditioners for systems of saddle point type. To verify the effectiveness of our new construction approach and H-matrix preconditioner scheme, we have applied them to solve various systems of linear equations arising from finite element methods and meshfree methods. The experimental results show that our preconditioners are competitive to other H-matrix preconditioners based on domain decomposition and existing preconditioners such as JOR and AMG preconditioners. Our H-matrix construction approach and preconditioner technique provide an alternative effective way to solve large systems of linear equations.

This thesis studies the use of the Source Reconstruction Method (SRM) to characterize antennas. The SRM calculates equivalent sources/currents on an arbitrarily shaped reconstruction surface to represent the original antenna. This is done by enforcing that the original antenna and equivalent currents radiate the same field at user selected measurement locations. These equivalent currents spatially characterize the original antenna because they can be used in direct radiation problems to obtain field estimates anywhere outside the reconstruction surface, including the far-field. First a spherical SRM algorithm is implemented and the diagnostic capabilities of the SRM are also synthetically shown through an example with an array of elementary dipoles. It is then shown that the SRM compares well to pre-existing commercial antenna software over different frequencies and can also be used successfully with a partial dataset. It is demonstrated that the equivalent currents can also provide meaningful information with experimental data. Next the hierarchical matrix framework is studied in conjunction with the SRM to decrease the algorithm's memory requirement and increase the speed of execution. It is shown that it is beneficial to use the hierarchical matrix framework with the SRM when using Love's condition or with measured data on a surface very close to the reconstruction surface. The SRM is then used to obtain incident field estimates in microwave imaging systems. Using a 2D transverse magnetic framework, we show that even with the limited data available in typical microwave tomography setups the SRM can produce incident field estimates in the imaging domain. These estimates are then used along with an MR-GNI algorithm to image synthetic and experimental objects with uncalibrated measured data. / October 2016

Hierarchical matrix (H-matrix) techniques can be used to efficiently treat dense matrices. With an H-matrix, the storage requirements and performing all fundamental operations, namely matrix-vector multiplication, matrix-matrix multiplication and matrix inversion can be done in almost linear complexity. In this work, we tried to gain even further speedup for the H-matrix arithmetic by utilizing multiple processors. Our approach towards an H-matrix distribution relies on the splitting of the index set. The main results achieved in this work based on the index-wise H-distribution are: A highly scalable algorithm for the H-matrix truncation and matrix-vector multiplication, a scalable algorithm for the H-matrix matrix multiplication, a limited scalable algorithm for the H-matrix inversion for a large number of processors.

This research introduces a row compression and nested product decomposition of an nxn hierarchical representation of a rank structured matrix A, which extends the compression and nested product decomposition of a quasiseparable matrix. The hierarchical parameter extraction algorithm of a quasiseparable matrix is efficient, requiring only O(nlog(n))operations, and is proven backward stable. The row compression is comprised of a sequence of small Householder transformations that are formed from the low-rank, lower triangular, off-diagonal blocks of the hierarchical representation. The row compression forms a factorization of matrix A, where A = QC, Q is the product of the Householder transformations, and C preserves the low-rank structure in both the lower and upper triangular parts of matrix A. The nested product decomposition is accomplished by applying a sequence of orthogonal transformations to the low-rank, upper triangular, off-diagonal blocks of the compressed matrix C. Both the compression and decomposition algorithms are stable, and require O(nlog(n)) operations. At this point, the matrix-vector product and solver algorithms are the only ones fully proven to be backward stable for quasiseparable matrices. By combining the fast matrix-vector product and system solver, linear systems involving the hierarchical representation to nested product decomposition are directly solved with linear complexity and unconditional stability. Applications in image deblurring and compression, that capitalize on the concepts from the row compression and nested product decomposition algorithms, will be shown.

In dieser Arbeit entwickeln wir einen Löser für elliptische Randwertprobleme. Dazu diskretisieren wir ein Randwertproblem mittels der Methode der finiten Elemente und erhalten ein Gleichungssystem. Mittels Gebietszerlegungstechniken unterteilen wir das Gebiet der Differentialgleichung und können Teilprobleme des Randwertproblems definieren. Durch die Gebietszerlegung wird eine Hierarchie von Zerlegungen definiert, die wir mittels eines Gebietszerlegungsbaumes festhalten. Anhand dieses Baumes definieren wir nun einen Löser für das Randwertproblem. Dabei berechnen wir die verschiedenen Matrizen des Lösers durch den sogenannten HDD-Algorithmus (engl. hierarchical domain decomposition). Die meisten der zu erstellenden Matrizen sind dabei vollbesetzt, weshalb wir sie mittels Hierarchischer Matrizen approximieren. Mit Hilfe der Hierarchischen Matrizen können wir die Matrizen mit einem fast linearen Aufwand erstellen und speichern. Der Aufwand der Matrixoperationen ist ebenfalls fast linear. Damit wir die Hierarchischen Matrizen für den HDD-Algorithmus verwenden können, müssen wir die Technik der Hierarchischen Matrizen erweitern. Unter anderem führen wir eingeschränkte Clusterbäume, eingeschränkte Blockclusterbäume und die verallgemeinerte Addition für Hierarchische Matrizen ein. Zusätzlich führen wir eine neue Clusterbaum-Konstruktion ein, die auf den HDD-Algorithmus zugeschnitten ist. Die Kombination des HDD-Algorithmus mit Hierarchischen Matrizen liefert einen Löser, den wir mit einem fast linearen Aufwand berechnen können. Der Aufwand zur Berechnung einer Lösung sowie der Speicheraufwand ist ebenfalls fast linear. Des Weiteren geben wir noch einige Modifizierungen des HDD-Algorithmus für weitere Anwendungsmöglichkeiten an. Zusätzlich diskutieren wir die Möglichkeiten der Parallelisierung, denn durch die Verwendung der Gebietszerlegung wird das Randwertproblem in unabhängige Teilprobleme aufgeteilt, die sich sehr gut parallelisieren lassen. Wir schließen die Arbeit mit numerischen Tests ab, die die theoretischen Aussagen bestätigen.

This dissertation introduces a novel fast direct solver and preconditioner for the solution of block tridiagonal linear systems that arise from the discretization of elliptic partial differential equations on a Cartesian product mesh, such as the variable-coefficient Poisson equation, the convection-diffusion equation, and the wave Helmholtz equation in heterogeneous media. The algorithm extends the traditional cyclic reduction method with hierarchical matrix techniques. The resulting method exposes substantial concurrency, and its arithmetic operations and memory consumption grow only log-linearly with problem size, assuming bounded rank of off-diagonal matrix blocks, even for problems with arbitrary coefficient structure. The method can be used as a standalone direct solver with tunable accuracy, or as a black-box preconditioner in conjunction with Krylov methods. The challenges that distinguish this work from other thrusts in this active field are the hybrid distributed-shared parallelism that can demonstrate the algorithm at large-scale, full three-dimensionality, and the three stressors of the current state-of-the-art multigrid technology: high wavenumber Helmholtz (indefiniteness), high Reynolds convection (nonsymmetry), and high contrast diffusion (inhomogeneity). Numerical experiments corroborate the robustness, accuracy, and complexity claims and provide a baseline of the performance and memory footprint by comparisons with competing approaches such as the multigrid solver hypre, and the STRUMPACK implementation of the multifrontal factorization with hierarchically semi-separable matrices. The companion implementation can utilize many thousands of cores of Shaheen, KAUST's Haswell-based Cray XC-40 supercomputer, and compares favorably with other implementations of hierarchical solvers in terms of time-to-solution and memory consumption.

Large dense matrices are ubiquitous in scientific computing, arising from the discretization of integral operators associated with elliptic pdes, Schur complement methods, covariances in spatial statistics, kernel-based machine learning, and numerical optimization problems. Hierarchical matrices are an efficient way for storing the dense matrices of very large dimension that appear in these and related settings. They exploit the fact that the underlying matrices, while formally dense, are data sparse. They have a structure consisting of blocks many of which can be well-approximated by low rank factorizations. A hierarchical organization of the blocks avoids superlinear growth in memory requirements to store n × n dense matrices in a scalable manner, requiring O(n) units of storage with a constant depending on a representative rank k for the low rank blocks. The asymptotically optimal storage requirement of the resulting hierarchical matrices is a critical advantage, particularly in extreme computing environments, characterized by low memory per processing core. The challenge then becomes to develop the parallel linear algebra operations that can be performed directly on this compressed representation. In this dissertation, I implement a set of hierarchical basic linear algebra subroutines (HBLAS) optimized for GPUs, including hierarchical matrix vector multiplication, orthogonalization, compression, low rank updates, and matrix multiplication. I develop a library of open source batched kernel operations previously missing on GPUs for the high performance implementation of the H2 operations, while relying wherever possible on existing open source and vendor kernels to ride future improvements in the technology. Fast marshaling routines extract the batch operation data from an efficient representation of the trees that compose the hierarchical matrices. The methods developed for GPUs extend to CPUs using the same code base with simple abstractions around the batched routine execution. To demonstrate the scalability of the hierarchical operations I implement a distributed memory multi-GPU hierarchical matrix vector product that focuses on reducing communication volume and hiding communication overhead and areas of low GPU utilization using low priority streams. Two demonstrations involving Hessians of inverse problems governed by pdes and space-fractional diffusion equations show the effectiveness of the hierarchical operations in realistic applications.

La méthode PEEC classique repose sur une méthode intégrale semi-analytique pour permettre la détermination d'un schéma électrique équivalent à l'aide de constantes localisées. Cette méthode est particulièrement bien adaptée pour la modélisation de régions conductrices du type filaire. S'il est possible actuellement de prendre en compte dans cette méthode des régions minces conductrices, cette approche demeure limitée et insatisfaisante. En effet, des contraintes très fortes pèsent sur les maillages qu'il est possible de traiter (discrétisation des géométries en quadrangles) et l'approche est limitée en fréquence (pas d'effet capacitif). L'objectif de cette thèse est d'introduire les effets capacitifs mais aussi magnétiques dans la méthode PEEC afin d'accéder à un outil général, performant et utilisable au niveau industriel. En particulier, la généralité de la formulation et sa flexibilité devrait permettre une utilisation simple à l'utilisateur du logiciel InCa3D non expert en méthode numérique. Le travail consistera donc à : • Consolider les travaux précédents par l'élargissement de l'approche, notamment en introduisant les effets capacitifs et magnétiques dans les formulations. • Proposer des méthodes de compressions matricielles adaptées pour limiter les temps de calcul et sauvegarder de la mémoire. / The classical method PEEC is based on a semi-analytical integral method to construct an equivalent electric circuit using lumped components. This method is particularly well-suited to modelise filiform conductors. It is actually possible to consider thin conductive regions with this method, but it's still limited and unsatisfactory. In fact, the meshes that can be used are very constrained (geometrically discretized by quadrangles) and the frequency approach is limited (capacitive effect is neglected). The aim of this thesis is to introduce the capacitive and magnetic effects into the method PEEC to get a general tool, efficient and usable at the industry level. Particularly, the generality of the formulation and its flexibility should enable a simple use of the software InCad3D for non-expert user on numerical methods. The work consists in : • Improving the last works by introducing the capacitive and magnetic effects in the formulations. • Suggesting some methods of matricial compression to improve the efficiency of the computation, and to lower the needed memory.

In dieser Arbeit entwickeln wir einen Löser für elliptische Randwertprobleme. Dazu diskretisieren wir ein Randwertproblem mittels der Methode der finiten Elemente und erhalten ein Gleichungssystem. Mittels Gebietszerlegungstechniken unterteilen wir das Gebiet der Differentialgleichung und können Teilprobleme des Randwertproblems definieren. Durch die Gebietszerlegung wird eine Hierarchie von Zerlegungen definiert, die wir mittels eines Gebietszerlegungsbaumes festhalten. Anhand dieses Baumes definieren wir nun einen Löser für das Randwertproblem. Dabei berechnen wir die verschiedenen Matrizen des Lösers durch den sogenannten HDD-Algorithmus (engl. hierarchical domain decomposition). Die meisten der zu erstellenden Matrizen sind dabei vollbesetzt, weshalb wir sie mittels Hierarchischer Matrizen approximieren. Mit Hilfe der Hierarchischen Matrizen können wir die Matrizen mit einem fast linearen Aufwand erstellen und speichern. Der Aufwand der Matrixoperationen ist ebenfalls fast linear. Damit wir die Hierarchischen Matrizen für den HDD-Algorithmus verwenden können, müssen wir die Technik der Hierarchischen Matrizen erweitern. Unter anderem führen wir eingeschränkte Clusterbäume, eingeschränkte Blockclusterbäume und die verallgemeinerte Addition für Hierarchische Matrizen ein. Zusätzlich führen wir eine neue Clusterbaum-Konstruktion ein, die auf den HDD-Algorithmus zugeschnitten ist. Die Kombination des HDD-Algorithmus mit Hierarchischen Matrizen liefert einen Löser, den wir mit einem fast linearen Aufwand berechnen können. Der Aufwand zur Berechnung einer Lösung sowie der Speicheraufwand ist ebenfalls fast linear. Des Weiteren geben wir noch einige Modifizierungen des HDD-Algorithmus für weitere Anwendungsmöglichkeiten an. Zusätzlich diskutieren wir die Möglichkeiten der Parallelisierung, denn durch die Verwendung der Gebietszerlegung wird das Randwertproblem in unabhängige Teilprobleme aufgeteilt, die sich sehr gut parallelisieren lassen. Wir schließen die Arbeit mit numerischen Tests ab, die die theoretischen Aussagen bestätigen.