# Methods to improve the efficiency and the memory requirement of LU factorization for complex symmetric system matrix

I want to solve a linear set of equations (Ax=b) using LU decomposition. My "A" matrix is a complex matrix which is symmetrical. The code that I work has two parts. In the first part I do all the initialization where I compute L and U factors of the matrix. The second part of the code runs in every time-step (which is specified in the beginning). In this section, I solve the equations Ld=b and Ux=d to find the solution vector x. And the computer which runs this part has limited memory. Also, I want this part to be as efficient as possible.

So my questions are:

1. Is there a way to save some memory for storing the L and U matrices of a symmetric complex matrix? If I deal with just an inverse instead of L and U I can just store the half of the matrix. Is there a similar way to save some storage for L and U matrices.

2. What are the methods that I can use to improve the efficiency of a LU decomposition for a complex symmetric matrix?

• Is using Cholesky factorization an option? (meaning a bit more is required from the matrix itself) Oct 13 '20 at 21:45
• Is the matrix sparse or dense? Oct 14 '20 at 1:15
• @Anton Menshov: I was under the impression that I cant use Cholesky for complex symmetric matrices @ Wolfgang Bangerth: The matrix is a sparse matrix
– HKK
Oct 14 '20 at 15:32
• There is a Cholesky-like algorithm that can be applied to a complex-symmetric matrix (really, it's identically the Cholesky algorithm, but using the complex-valued version of square root). It is not very stable so I wouldn't recommend it as a black box solver for an arbitrary complex-symmetric matrix. But the Bunch-Kauffman algorithm [used by all these LDL' factorizations] is about as stable as partial-pivoted LU (ie good enough for general purpose use). Oct 14 '20 at 15:59

You probably want a factorization of the form $$\mathbf A = \mathbf L \mathbf D \mathbf L^T$$, it can certainly be applied to a complex symmetric $$\mathbf A$$. LAPACK implements this factorization within [zsytrf] and provides a corresponding backsolution routine within [zsytrs]. There are sparse-direct versions of this algorithm too. PARDISO, TAUCS, and MyraMath all implement it (disclaimer: I wrote that last one).
EDIT1: Regarding the efficiency of backsolution, it's probably not great. Unlike LU [zgetrf] and Cholesky [zpotrf], the algorithm used by [zsytrf] technically does not deliver a triangular factor that is layout compatible with the triangular routines of the BLAS (eg [ztrsm]). Instead, it stores $$\mathbf L$$ and $$\mathbf D$$ interleaved as a bunch of 1x1 and 2x2 blocks (kinda arbitrarily, based on the pivoting decisions), which means the backsolution requires a similar sequence of rank-1 and rank-2 steps (this fussiness of the backsolution process is why LAPACK provides [zsytrs] to begin with). Unfortunately this is all BLAS1/BLAS2 level of performance. The algorithm to untangle $$\mathbf L$$ into a BLAS3-compatible triangular matrix is tedious.