tridiag_complex.hpp

Go to the documentation of this file.

/// @file factorization/tridiag_complex.hpp
/// @brief Precomputed LU Thomas solver for constant-coefficient complex
/// tridiagonal systems.
///
/// Solves the n x n system:
///
///   a*x[k-1] + b*x[k] + c*x[k+1] = d[k],  k = 0 .. n-1
///
/// with Dirichlet boundary conditions  x[-1] = x[n] = 0.
///
/// Because a, b, c are the same for every row the LU factorization depends only
/// on the coefficients, not on d.  Factor once, then call solve() for each
/// new right-hand side in O(n).
///
/// Typical use: Crank-Nicolson kinetic sweeps in the TDSE solver, where
/// a = c = -i*alpha  and  b = 1 + 2i*alpha  for some real alpha.
#pragma once
 
#include "core/types.hpp"
#include <complex>
#include <vector>
 
namespace num {
 
/// @brief Constant-coefficient complex tridiagonal solver (precomputed LU).
///
/// @code
/// num::ComplexTriDiag td;
/// std::complex<double> a(0.0, -alpha);
/// std::complex<double> b(1.0,  2.0*alpha);
/// std::complex<double> c(0.0, -alpha);
/// td.factor(n, a, b, c);
///
/// std::vector<std::complex<double>> rhs = ...;
/// td.solve(rhs);   // rhs is overwritten with the solution
/// @endcode
struct ComplexTriDiag {
  using cplx = std::complex<double>;
 
  std::vector<cplx> c_mod; ///< Modified super-diagonal (precomputed from LU)
  std::vector<cplx> inv_b; ///< Inverse of modified main diagonal (precomputed)
  int n = 0;
  cplx a_coeff = {}; ///< Sub-diagonal value (constant across all rows)
 
  /// @brief Factor the tridiagonal matrix.
  ///
  /// @param n_  System size
  /// @param a_  Sub-diagonal coefficient (row k depends on x[k-1])
  /// @param b_  Main-diagonal coefficient
  /// @param c_  Super-diagonal coefficient (row k depends on x[k+1])
  void factor(int n_, cplx a_, cplx b_, cplx c_);
 
  /// @brief In-place Thomas solve.
  ///
  /// On entry  d holds the right-hand side; on exit it holds the solution.
  /// The size of d must equal n (set by the last factor() call).
  void solve(std::vector<cplx> &d) const;
};
 
} // namespace num