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?

  • $\begingroup$ Is using Cholesky factorization an option? (meaning a bit more is required from the matrix itself) $\endgroup$ – Anton Menshov Oct 13 at 21:45
  • $\begingroup$ Is the matrix sparse or dense? $\endgroup$ – Wolfgang Bangerth Oct 14 at 1:15
  • $\begingroup$ @Anton Menshov: I was under the impression that I cant use Cholesky for complex symmetric matrices @ Wolfgang Bangerth: The matrix is a sparse matrix $\endgroup$ – HKK Oct 14 at 15:32
  • 1
    $\begingroup$ 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). $\endgroup$ – rchilton1980 Oct 14 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.

EDIT2: If your input is sparse, I'd just use a package that handles all this. Start with PARDISO, it's already present in MKL. It's probably not worth digging into any of the details.

| cite | improve this answer | |
  • $\begingroup$ LDLT would surely save some memory for me. Can you comment on how effective the forward and backward substitutions of LDLT regard to efficiency. $\endgroup$ – HKK Oct 14 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.