aligator  0.15.0
A primal-dual augmented Lagrangian-type solver for nonlinear trajectory optimization.
Loading...
Searching...
No Matches
block-tridiagonal.hpp
Go to the documentation of this file.
1
3#pragma once
4
7#include "aligator/tracy.hpp"
8
9namespace aligator {
10namespace gar {
11
12template <typename MatrixType>
13typename MatrixType::PlainObject
14blockTridiagToDenseMatrix(const std::vector<MatrixType> &subdiagonal,
15 const std::vector<MatrixType> &diagonal,
16 const std::vector<MatrixType> &superdiagonal) {
17 if (subdiagonal.size() != superdiagonal.size() ||
18 diagonal.size() != superdiagonal.size() + 1) {
19 ALIGATOR_DOMAIN_ERROR("Wrong lengths");
20 }
21
22 using PlainObjectType = typename MatrixType::PlainObject;
23
24 const size_t N = subdiagonal.size();
25 Eigen::Index dim = 0;
26 for (size_t i = 0; i <= N; i++) {
27 dim += diagonal[i].cols();
28 }
29
30 PlainObjectType out(dim, dim);
31 out.setZero();
32 Eigen::Index i0 = 0;
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];
36
37 if (i != N) {
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];
41 }
42
43 i0 += d;
44 }
45 return out;
46}
47
50template <typename MatrixType, typename InputType, typename OutType,
51 typename Scalar = typename MatrixType::Scalar>
52bool blockTridiagMatMul(const std::vector<MatrixType> &Asub,
53 const std::vector<MatrixType> &Adiag,
54 const std::vector<MatrixType> &Asuper,
56 BlkMatrix<OutType, -1, 1> &c, const Scalar beta) {
58 const size_t N = Asuper.size();
59
60 c.matrix() *= beta;
61
62 c[0].noalias() += Adiag[0] * b[0];
63 c[0].noalias() += Asuper[0] * b[1];
64
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];
69 }
70
71 c[N].noalias() += Asub[N - 1] * b[N - 1];
72 c[N].noalias() += Adiag[N] * b[N];
73
74 return true;
75}
76
80template <typename MatrixType, typename RhsType, typename DecType>
81bool symmetricBlockTridiagSolve(std::vector<MatrixType> &subdiagonal,
82 std::vector<MatrixType> &diagonal,
83 const std::vector<MatrixType> &superdiagonal,
85 std::vector<DecType> &facs) {
86 ALIGATOR_TRACY_ZONE_SCOPED;
88
89 if (subdiagonal.size() != superdiagonal.size() ||
90 diagonal.size() != superdiagonal.size() + 1 ||
91 rhs.rowDims().size() != diagonal.size()) {
92 return false;
93 }
94
95 // size of problem
96 const size_t N = superdiagonal.size();
97
98 size_t i = N - 1;
99 while (true) {
100 DecType &ldl = facs[i + 1];
101 ldl.compute(diagonal[i + 1]);
102 if (ldl.info() != Eigen::Success)
103 return false;
104
105 Eigen::Ref<RhsType> r = rhs[i + 1];
106 ldl.solveInPlace(r);
107
108 // the math has index of B starting at 1, array starts at 0
109 const auto &Bip1 = superdiagonal[i];
110 auto &Cip1 = subdiagonal[i]; // should be Bi.transpose()
111
112 rhs[i].noalias() -= Bip1 * rhs[i + 1];
113 ldl.solveInPlace(Cip1); // contains U.T = D[i+1]^-1 * B[i+1].transpose()
114
115 diagonal[i].noalias() -= Bip1 * Cip1;
116
117 if (i == 0)
118 break;
119 i--;
120 }
121
122 {
123 DecType &ldl = facs[0];
124 ldl.compute(diagonal[0]);
125 if (ldl.info() != Eigen::Success)
126 return false;
127 Eigen::Ref<RhsType> r = rhs[0];
128 ldl.solveInPlace(r);
129 }
130
131 for (size_t i = 0; i < N; i++) {
132 const auto &Uip1t = subdiagonal[i];
133 rhs[i + 1].noalias() -= Uip1t * rhs[i];
134 }
135
136 return true;
137}
138
144template <typename MatrixType, typename RhsType, typename DecType>
145bool blockTridiagRefinementStep(const std::vector<MatrixType> &transposedUfacs,
146 const std::vector<MatrixType> &superdiagonal,
147 const std::vector<DecType> &diagonalFacs,
150 // size of problem
151 const size_t N = superdiagonal.size();
152
153 size_t i = N - 1;
154 while (true) {
155 const DecType &ldl = diagonalFacs[i + 1];
156 Eigen::Ref<RhsType> r = rhs[i + 1];
157 ldl.solveInPlace(r);
158
159 // the math has index of B starting at 1, array starts at 0
160 const auto &Bip1 = superdiagonal[i];
161
162 rhs[i].noalias() -= Bip1 * rhs[i + 1];
163
164 if (i == 0)
165 break;
166 i--;
167 }
168
169 const DecType &ldl = diagonalFacs[0];
170 Eigen::Ref<RhsType> r = rhs[0];
171 ldl.solveInPlace(r);
172
173 for (size_t i = 0; i < N; i++) {
174 const auto &Uip1t = transposedUfacs[i];
175 rhs[i + 1].noalias() -= Uip1t * rhs[i];
176 }
177
178 return true;
179}
180
183template <typename MatrixType, typename RhsType, typename DecType>
185 std::vector<MatrixType> &subdiagonal, std::vector<MatrixType> &diagonal,
186 std::vector<MatrixType> &superdiagonal, BlkMatrix<RhsType, -1, 1> &rhs,
187 std::vector<DecType> &facs) {
188 ALIGATOR_TRACY_ZONE_SCOPED;
190
191 if (subdiagonal.size() != superdiagonal.size() ||
192 diagonal.size() != superdiagonal.size() + 1 ||
193 rhs.rowDims().size() != diagonal.size()) {
194 return false;
195 }
196
197 // size of problem
198 const size_t N = superdiagonal.size();
199
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)
204 return false;
205
206 Eigen::Ref<RhsType> r = rhs[i];
207 ldl.solveInPlace(r);
208
209 // the math has index of B starting at 1, array starts at 0
210 auto &Bip1 = superdiagonal[i];
211 const auto &Cip1 = subdiagonal[i]; // should be B[i+1].transpose()
212
213 rhs[i + 1].noalias() -= Cip1 * rhs[i];
214 ldl.solveInPlace(Bip1); // contains L.T = D[i]^-1 * B[i+1]
215
216 // substract B[i+1].T * L.T
217 diagonal[i + 1].noalias() -= Cip1 * Bip1;
218 }
219
220 {
221 DecType &ldl = facs[N];
222 ldl.compute(diagonal[N]);
223 if (ldl.info() != Eigen::Success)
224 return false;
225 Eigen::Ref<RhsType> r = rhs[N];
226 ldl.solveInPlace(r);
227 }
228
229 size_t i = N - 1;
230 while (true) {
231 const auto &Lim1t = superdiagonal[i];
232 rhs[i].noalias() -= Lim1t * rhs[i + 1];
233
234 if (i == 0)
235 break;
236 i--;
237 }
238 return true;
239}
240
241} // namespace gar
242} // namespace aligator
Block matrix class, with a fixed-size number of row and column blocks.
const row_dim_t & rowDims() const
MatrixType & matrix()
#define ALIGATOR_DOMAIN_ERROR(...)
#define ALIGATOR_NOMALLOC_SCOPED
Definition math.hpp:44
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.
Main package namespace.