# Numerical Linear Algebra Problems in Structural Analysis

Kannan, Ramaseshan (2014) Numerical Linear Algebra Problems in Structural Analysis. Doctoral thesis, Manchester Institute for Mathematical Sciences, The University of Manchester.

A range of numerical linear algebra problems that arise in finite element-based structural analysis are considered. These problems were encountered when implementing the finite element method in the software package Oasys GSA. We present novel solutions to these problems in the form of a new method for error detection, algorithms with superior numerical efficiency and algorithms with scalable performance on parallel computers. The solutions and their corresponding software implementations have been integrated into GSA's program code and we present results that demonstrate the use of these implementations by engineers to solve real-world structural analysis problems. We start by introducing a new method for detecting the sources of ill conditioning in finite element stiffness matrices. The stiffness matrix can become ill-conditioned due to a variety of errors including user errors, i.e., errors in the definition of the finite element model when using a software package. We identify two classes of errors and develop a new method called model stability analysis for detecting where in the structural model these errors arise. Our method for debugging models uses the properties of the eigenvectors associated with the smallest and largest eigenvalues of the stiffness matrix to identify parts of the structural model that are badly defined. Through examples, we demonstrate the use of model stability analysis on real-life models and show that upon fixing the errors identified by the method, the condition number of the matrices typically drops from $O(10^{16})$ to $O(10^8)$. In the second part we introduce the symmetric definite generalized eigenproblem and the symmetric eigenproblem that are encountered in structural analysis. The symmetric definite generalized eigenvalue problem is to solve for the smallest eigenvalues and associated eigenvectors of $Kx = \lambda Mx$ for a symmetric positive-definite, sparse stiffness matrix $K$ and a symmetric semidefinite mass matrix $M$. %We seek the eigenvalues and associated eigenvectors closest to 0 for this problem. We analyse the existing solution algorithm, which uses the subspace iteration method. We improve the numerical efficiency of the algorithm by accelerating convergence via novel shift strategies and by reducing floating point operations using selective re-orthogonalization and locking of converged eigenvectors. The software implementation also benefits from improvements such as the use of faster and more robust libraries for linear algebra operations encountered during the iteration. We then develop an implementation of the subspace iteration method for solving the symmetric eigenvalue problem $K x = \lambda x$ for the smallest and largest eigenvalues and their eigenvectors; this problem arises in model stability analysis. In the next part we tackle the main bottleneck of sparse matrix computations in our eigensolvers: the sparse matrix-multiple vector multiplication kernel. We seek to obtain an algorithm that has high computational throughput for the operation $Y = AX$ for a square sparse $A$ and a conforming skinny, dense $X$ which, in turn, depends on the underlying storage format of $A$. Current sparse matrix formats and algorithms have high bandwidth requirements and poor reuse of cache and register loaded entries, which restricts their performance. We propose the mapped blocked row format: a bitmapped sparse matrix format that stores entries as blocks without a fill overhead, thereby offering blocking without additional storage and bandwidth overheads. An efficient algorithm decodes bitmaps using de Bruijn sequences and minimizes the number of conditionals evaluated. Performance is compared with that of popular formats, including vendor implementations of sparse BLAS. Our sparse matrix multiple-vector multiplication algorithm achieves high throughput on all platforms and is implemented using platform-neutral optimizations. The last part of the thesis deals with the $B$-orthogonalization problem, i.e., for a sparse symmetric positive definite $B \in \R^{\nbyn}$ and a tall skinny matrix $X \in \R^{\nbyl}$, we wish to rotate the columns of $X$ such that $X^TBX = I$. This problem arises when the engineer wishes to orthonormalize the eigenvectors of the generalized eigenproblem such that they are orthogonal with respect to the stiffness matrix. Conventionally in the literature Gram-Schmidt like methods are employed for \bo ization but they have poor cache efficiency. We recall a method that uses the Cholesky factorization of the inner product matrix $X^TBX$ and derive a stable algorithm that increases the parallel scalability through cache-reuse and is also numerically stable. Our experiments demonstrate orders of magnitude improvement in performance compared with Gram-Schmidt family of methods.