7#include "aligator/tracy.hpp"
12template <
typename MatrixType,
class A>
13typename MatrixType::PlainObject
15 const std::vector<MatrixType, A> &diagonal,
16 const std::vector<MatrixType, A> &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,
class A,
typename InputType,
typename OutType,
51 typename Scalar =
typename MatrixType::Scalar>
53 const std::vector<MatrixType, A> &Adiag,
54 const std::vector<MatrixType, A> &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,
class A,
class A2,
typename RhsType,
83 std::vector<MatrixType, A> &diagonal,
84 const std::vector<MatrixType, A> &superdiagonal,
86 std::vector<DecType, A2> &facs) {
87 ALIGATOR_TRACY_ZONE_SCOPED;
90 if (subdiagonal.size() != superdiagonal.size() ||
91 diagonal.size() != superdiagonal.size() + 1 ||
92 rhs.
rowDims().size() != diagonal.size()) {
97 const size_t N = superdiagonal.size();
101 DecType &ldl = facs[i + 1];
102 ldl.compute(diagonal[i + 1]);
103 if (ldl.info() != Eigen::Success)
106 Eigen::Ref<RhsType> r = rhs[i + 1];
110 const auto &Bip1 = superdiagonal[i];
111 auto &Cip1 = subdiagonal[i];
113 rhs[i].noalias() -= Bip1 * rhs[i + 1];
114 ldl.solveInPlace(Cip1);
116 diagonal[i].noalias() -= Bip1 * Cip1;
124 DecType &ldl = facs[0];
125 ldl.compute(diagonal[0]);
126 if (ldl.info() != Eigen::Success)
128 Eigen::Ref<RhsType> r = rhs[0];
132 for (
size_t i = 0; i < N; i++) {
133 const auto &Uip1t = subdiagonal[i];
134 rhs[i + 1].noalias() -= Uip1t * rhs[i];
145template <
typename MatrixType,
class A,
class A2,
typename RhsType,
148 const std::vector<MatrixType, A> &transposedUfacs,
149 const std::vector<MatrixType, A> &superdiagonal,
150 const std::vector<DecType, A2> &diagonalFacs,
154 const size_t N = superdiagonal.size();
158 const DecType &ldl = diagonalFacs[i + 1];
159 Eigen::Ref<RhsType> r = rhs[i + 1];
163 const auto &Bip1 = superdiagonal[i];
165 rhs[i].noalias() -= Bip1 * rhs[i + 1];
172 const DecType &ldl = diagonalFacs[0];
173 Eigen::Ref<RhsType> r = rhs[0];
176 for (
size_t i = 0; i < N; i++) {
177 const auto &Uip1t = transposedUfacs[i];
178 rhs[i + 1].noalias() -= Uip1t * rhs[i];
186template <
typename MatrixType,
class A,
class A2,
typename RhsType,
189 std::vector<MatrixType, A> &subdiagonal,
190 std::vector<MatrixType, A> &diagonal,
192 std::vector<DecType, A2> &facs) {
193 ALIGATOR_TRACY_ZONE_SCOPED;
196 if (subdiagonal.size() != superdiagonal.size() ||
197 diagonal.size() != superdiagonal.size() + 1 ||
198 rhs.
rowDims().size() != diagonal.size()) {
203 const size_t N = superdiagonal.size();
205 for (
size_t i = 0; i < N; i++) {
206 DecType &ldl = facs[i];
207 ldl.compute(diagonal[i]);
208 if (ldl.info() != Eigen::Success)
211 Eigen::Ref<RhsType> r = rhs[i];
215 auto &Bip1 = superdiagonal[i];
216 const auto &Cip1 = subdiagonal[i];
218 rhs[i + 1].noalias() -= Cip1 * rhs[i];
219 ldl.solveInPlace(Bip1);
222 diagonal[i + 1].noalias() -= Cip1 * Bip1;
226 DecType &ldl = facs[N];
227 ldl.compute(diagonal[N]);
228 if (ldl.info() != Eigen::Success)
230 Eigen::Ref<RhsType> r = rhs[N];
236 const auto &Lim1t = superdiagonal[i];
237 rhs[i].noalias() -= Lim1t * rhs[i + 1];
Block matrix class, with a fixed or dynamic-size number of row and column blocks.
const RowDimsType & rowDims() const
#define ALIGATOR_DOMAIN_ERROR(...)
#define ALIGATOR_NOMALLOC_SCOPED
bool symmetricBlockTridiagSolve(std::vector< MatrixType, A > &subdiagonal, std::vector< MatrixType, A > &diagonal, const std::vector< MatrixType, A > &superdiagonal, BlkMatrix< RhsType, -1, 1 > &rhs, std::vector< DecType, A2 > &facs)
Solve a symmetric block-tridiagonal problem by in-place factorization. The subdiagonal will be used ...
bool symmetricBlockTridiagSolveDownLooking(std::vector< MatrixType, A > &subdiagonal, std::vector< MatrixType, A > &diagonal, std::vector< MatrixType, A > &superdiagonal, BlkMatrix< RhsType, -1, 1 > &rhs, std::vector< DecType, A2 > &facs)
Solve a symmetric block-tridiagonal problem by in-place factorization. The subdiagonal will be used ...
bool blockTridiagMatMul(const std::vector< MatrixType, A > &Asub, const std::vector< MatrixType, A > &Adiag, const std::vector< MatrixType, A > &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 ...
MatrixType::PlainObject blockTridiagToDenseMatrix(const std::vector< MatrixType, A > &subdiagonal, const std::vector< MatrixType, A > &diagonal, const std::vector< MatrixType, A > &superdiagonal)
bool blockTridiagRefinementStep(const std::vector< MatrixType, A > &transposedUfacs, const std::vector< MatrixType, A > &superdiagonal, const std::vector< DecType, A2 > &diagonalFacs, BlkMatrix< RhsType, -1, 1 > &rhs)
Given the decomposed matrix data, just perform the backward-forward step of the algorithm.