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.

In this master thesis, a 2D device simulator run on a hybrid classical-quantum computer was developed. The simulator was developed to treat statistical quantum effects such as quantum tunneling and quantum confinement in nanoscale transistors. The simulation scheme is based on a self-consistent solution of the coupled non-linear 2D SchrödingerPoisson equations. The Open Boundary Condition (OBC) of the Schrödinger equation, which allows for electrons to pass through the device between the leads (source and drain), are modeled with the QuantumTransmitting Boundary Method (QTBM). The differential equations are discretized with the finite-element method, using rectangular mesh elements. The self-consistent loop is a very time-consuming process, mainly due to the solution of the discretized OBC Schrödinger equation. To accelerate the solution time of the Schrödinger equation, a quantum assisted domain decomposition method is implemented. The domain decomposition method of choice is the Block Cyclic Reduction (BCR) method. The BCR method is at least 15 times faster (CPU time) than solving the whole linear system of equations with the Python solver numpy.linalg.solve, based on the LAPACK routine _gesv. In the project, we also propose an alternative approach of the BCR method called the "extra layer BCR" that shows an improved accuracy for certain types of solutions. In a quantum assisted version, the matrix inverse solver as a step in the BCR method was computed on the D-Wave quantum annealer chip ADVANTAGE_SYSTEM4.1 [4]. Two alternative methods to solve the matrix inverses on a quantum annealer were compared. One is called the "unit vector" approach, based on work by Rogers and Singleton [5], and the other is called the "whole matrix" approach which was developed in the thesis. Because of the limited amount of qubits available on the quantum annealer, the "unit vector" approach was more suitable for adaption in the BCR method. Comparing the quantum annealer to the Python inverse solver numpy.linalg.inv, also based on LAPACK, it was found that an accurate solution can be achieved, but the simulation time (CPU time) is at best 500 times slower than numpy.linalg.inv.

Cette thèse vise à la conception de nouveaux algorithmes rapides en calcul numérique via les matrices de Toeplitz. Tout d'abord, nous avons introduit un algorithme rapide sur le calcul de l'inverse d'une matrice triangulaire de Toeplitz en se basant sur des notions d'interpolation polynomiale. Cet algorithme nécessitant uniquement deux FFT(2n) est manifestement efficace par rapport à ses prédécésseurs. ensuite, nous avons introduit un algorithme rapide pour la résolution d'un système linéaire de Toeplitz bande. Cette approche est basée sur l'extension de la matrice donnée par plusieurs lignes en dessus, de plusieurs colonnes à droite et d'attribuer des zéros et des constantes non nulles dans chacune de ces lignes et de ces colonnes de telle façon que la matrice augmentée à la structure d'une matrice triangulaire inférieure de Toeplitz. La stabilité de l'algorithme a été discutée et son efficacité a été aussi justifiée. Finalement, nous avons abordé la résolution d'un système de Toeplitz bandes par blocs bandes de Toeplitz. Ceci étant primordial pour établir la connexion de nos algorithmes à des applications en restauration d'images, un domaine phare en mathématiques appliquées. / This thesis aims to design new fast algorithms for numerical computation via the Toeplitz matrices. First, we introduced a fast algorithm to compute the inverse of a triangular Toeplitz matrix with real and/or complex numbers based on polynomial interpolation techniques. This algorithm requires only two FFT (2n) is clearly effective compared to predecessors. A numerical accuracy and error analysis is also considered. Numerical examples are given to illustrate the effectiveness of our method. In addition, we introduced a fast algorithm for solving a linear banded Toeplitz system. This new approach is based on extending the given matrix with several rows on the top and several columns on the right and to assign zeros and some nonzero constants in each of these rows and columns in such a way that the augmented matrix has a lower triangular Toeplitz structure. Stability of the algorithm is discussed and its performance is showed by numerical experiments. This is essential to connect our algorithms to applications such as image restoration applications, a key area in applied mathematics.