2
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • $\begingroup$ Is using Cholesky factorization an option? (meaning a bit more is required from the matrix itself) $\endgroup$
    – Anton Menshov ♦
    Oct 13, 2020 at 21:45
  • $\begingroup$ Is the matrix sparse or dense? $\endgroup$
    – Wolfgang Bangerth
    Oct 14, 2020 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, 2020 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, 2020 at 15:59

1 Answer 1

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\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, 2020 at 15:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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