emg3d.core.solve(amat, bvec)[source]

Solve A x = b using a non-standard Cholesky factorisation.

Solve the system A x = b using a non-standard Cholesky factorisation without pivoting for a symmetric, complex matrix A tailored to the problem of the multigrid solver. The matrix A (amat) is an array of length 6*n, containing the main diagonal and the first five lower off-diagonals (ordered so that the first element of the main diagonal is followed by the first elements of the off diagonals, then the second elements and so on). The vector bvec has length b.

The solution is placed in b (bvec), and A (amat) is replaced by its decomposition.

  1. Non-standard Cholesky factorisation.

    From [Muld07]: «We use a non-standard Cholesky factorisation. The standard factorisation factors a Hermitian matrix A into L L^H, where L is a lower triangular matrix and L^H its complex conjugate transpose. In our case, the discretisation is based on the Finite Integration Technique ([Weil77]) and provides a matrix A that is complex-valued and symmetric: A = A^T, where the superscript T denotes the transpose. The line relaxation scheme takes a matrix B that is a subset of A along the line. B is a complex symmetric band matrix with eleven diagonals. The non-standard Cholesky factorisation factors the matrix B into L L^T. Because of the symmetry, only the main diagonal and five lower diagonal elements of B need to be computed. The Cholesky factorisation replaces this matrix by L, containing six diagonals, after which the line relaxation can be carried out by simple back-substitution.»

    \(A = L D L^T\) factorisation without pivoting:

    \[\begin{split}D(j) &= A(j,j)-\sum_{k=1}^{j-1} L(j,k)^2 D(k),\ j=1,..,n ;\\ L(i,j) &= \frac{1}{D(j)} \left[A(i,j)-\sum_{k=1}^{j-1} L(i,k)L(j,k)D(k)\right], \ i=j+1,..,n .\end{split}\]

    A and L are in this case arrays, where \(A(i, j) \rightarrow A(i+5j)\).

  2. Solve A x = b.

    Solve A x = b, given L which is the result from the factorisation in the first step (and stored in A), hence, solve L x = b, where x is stored in b:

    \[b(j) = b(j) - \sum_{k=1}^{j-1} L(j,k) x(k), j = 2,..,n .\]

The result is equivalent with simply using numpy.linalg.solve, but faster for the particular use-case of this code.

Note that in this custom solver there is no pivoting, and the diagonals of the matrix cannot be zero.


Banded matrix A provided as a vector of length 6*n, containing main diagonal plus first five lower diagonals.


Right-hand-side vector b of length n.