7#include "aligator/tracy.hpp"
12template <
typename MatrixType>
13typename MatrixType::PlainObject
15 const std::vector<MatrixType> &diagonal,
16 const std::vector<MatrixType> &superdiagonal) {
17 if (subdiagonal.size() != superdiagonal.size() ||
18 diagonal.size() != superdiagonal.size() + 1) {
22 using PlainObjectType =
typename MatrixType::PlainObject;
24 const size_t N = subdiagonal.size();
26 for (
size_t i = 0; i <= N; i++) {
27 dim += diagonal[i].cols();
30 PlainObjectType out(dim, dim);
33 for (
size_t i = 0; i <= N; i++) {
34 Eigen::Index d = diagonal[i].cols();
35 out.block(i0, i0, d, d) = diagonal[i];
38 Eigen::Index dn = superdiagonal[i].cols();
39 out.block(i0, i0 + d, d, dn) = superdiagonal[i];
40 out.block(i0 + d, i0, dn, d) = subdiagonal[i];
50template <
typename MatrixType,
typename InputType,
typename OutType,
51 typename Scalar =
typename MatrixType::Scalar>
53 const std::vector<MatrixType> &Adiag,
54 const std::vector<MatrixType> &Asuper,
58 const size_t N = Asuper.size();
62 c[0].noalias() += Adiag[0] * b[0];
63 c[0].noalias() += Asuper[0] * b[1];
65 for (
size_t i = 1; i < N; i++) {
66 c[i].noalias() += Asub[i - 1] * b[i - 1];
67 c[i].noalias() += Adiag[i] * b[i];
68 c[i].noalias() += Asuper[i] * b[i + 1];
71 c[N].noalias() += Asub[N - 1] * b[N - 1];
72 c[N].noalias() += Adiag[N] * b[N];
80template <
typename MatrixType,
typename RhsType,
typename DecType>
82 std::vector<MatrixType> &diagonal,
83 const std::vector<MatrixType> &superdiagonal,
85 std::vector<DecType> &facs) {
86 ALIGATOR_TRACY_ZONE_SCOPED;
89 if (subdiagonal.size() != superdiagonal.size() ||
90 diagonal.size() != superdiagonal.size() + 1 ||
91 rhs.
rowDims().size() != diagonal.size()) {
96 const size_t N = superdiagonal.size();
100 DecType &ldl = facs[i + 1];
101 ldl.compute(diagonal[i + 1]);
102 if (ldl.info() != Eigen::Success)
105 Eigen::Ref<RhsType> r = rhs[i + 1];
109 const auto &Bip1 = superdiagonal[i];
110 auto &Cip1 = subdiagonal[i];
112 rhs[i].noalias() -= Bip1 * rhs[i + 1];
113 ldl.solveInPlace(Cip1);
115 diagonal[i].noalias() -= Bip1 * Cip1;
123 DecType &ldl = facs[0];
124 ldl.compute(diagonal[0]);
125 if (ldl.info() != Eigen::Success)
127 Eigen::Ref<RhsType> r = rhs[0];
131 for (
size_t i = 0; i < N; i++) {
132 const auto &Uip1t = subdiagonal[i];
133 rhs[i + 1].noalias() -= Uip1t * rhs[i];
144template <
typename MatrixType,
typename RhsType,
typename DecType>
146 const std::vector<MatrixType> &superdiagonal,
147 const std::vector<DecType> &diagonalFacs,
151 const size_t N = superdiagonal.size();
155 const DecType &ldl = diagonalFacs[i + 1];
156 Eigen::Ref<RhsType> r = rhs[i + 1];
160 const auto &Bip1 = superdiagonal[i];
162 rhs[i].noalias() -= Bip1 * rhs[i + 1];
169 const DecType &ldl = diagonalFacs[0];
170 Eigen::Ref<RhsType> r = rhs[0];
173 for (
size_t i = 0; i < N; i++) {
174 const auto &Uip1t = transposedUfacs[i];
175 rhs[i + 1].noalias() -= Uip1t * rhs[i];
183template <
typename MatrixType,
typename RhsType,
typename DecType>
185 std::vector<MatrixType> &subdiagonal, std::vector<MatrixType> &diagonal,
187 std::vector<DecType> &facs) {
188 ALIGATOR_TRACY_ZONE_SCOPED;
191 if (subdiagonal.size() != superdiagonal.size() ||
192 diagonal.size() != superdiagonal.size() + 1 ||
193 rhs.
rowDims().size() != diagonal.size()) {
198 const size_t N = superdiagonal.size();
200 for (
size_t i = 0; i < N; i++) {
201 DecType &ldl = facs[i];
202 ldl.compute(diagonal[i]);
203 if (ldl.info() != Eigen::Success)
206 Eigen::Ref<RhsType> r = rhs[i];
210 auto &Bip1 = superdiagonal[i];
211 const auto &Cip1 = subdiagonal[i];
213 rhs[i + 1].noalias() -= Cip1 * rhs[i];
214 ldl.solveInPlace(Bip1);
217 diagonal[i + 1].noalias() -= Cip1 * Bip1;
221 DecType &ldl = facs[N];
222 ldl.compute(diagonal[N]);
223 if (ldl.info() != Eigen::Success)
225 Eigen::Ref<RhsType> r = rhs[N];
231 const auto &Lim1t = superdiagonal[i];
232 rhs[i].noalias() -= Lim1t * rhs[i + 1];
Block matrix class, with a fixed-size number of row and column blocks.
const row_dim_t & rowDims() const
#define ALIGATOR_DOMAIN_ERROR(...)
#define ALIGATOR_NOMALLOC_SCOPED
bool symmetricBlockTridiagSolveDownLooking(std::vector< MatrixType > &subdiagonal, std::vector< MatrixType > &diagonal, std::vector< MatrixType > &superdiagonal, BlkMatrix< RhsType, -1, 1 > &rhs, std::vector< DecType > &facs)
Solve a symmetric block-tridiagonal problem by in-place factorization. The subdiagonal will be used ...
MatrixType::PlainObject blockTridiagToDenseMatrix(const std::vector< MatrixType > &subdiagonal, const std::vector< MatrixType > &diagonal, const std::vector< MatrixType > &superdiagonal)
bool blockTridiagMatMul(const std::vector< MatrixType > &Asub, const std::vector< MatrixType > &Adiag, const std::vector< MatrixType > &Asuper, const BlkMatrix< InputType, -1, 1 > &b, BlkMatrix< OutType, -1, 1 > &c, const Scalar beta)
Evaluate c <- beta * c + A * b, where A is the block-tridiagonal matrix described by the first three ...
bool symmetricBlockTridiagSolve(std::vector< MatrixType > &subdiagonal, std::vector< MatrixType > &diagonal, const std::vector< MatrixType > &superdiagonal, BlkMatrix< RhsType, -1, 1 > &rhs, std::vector< DecType > &facs)
Solve a symmetric block-tridiagonal problem by in-place factorization. The subdiagonal will be used ...
bool blockTridiagRefinementStep(const std::vector< MatrixType > &transposedUfacs, const std::vector< MatrixType > &superdiagonal, const std::vector< DecType > &diagonalFacs, BlkMatrix< RhsType, -1, 1 > &rhs)
Given the decomposed matrix data, just perform the backward-forward step of the algorithm.