Line data Source code
1 : #ifndef BC7FD90F_FBB7_481C_89C4_89BEE41309C5
2 : #define BC7FD90F_FBB7_481C_89C4_89BEE41309C5
3 :
4 : #include "base/first_include.h" // IWYU pragma: keep
5 : #include "math/linalg/square_matrix.h" // IWYU pragma: keep
6 : #include <initializer_list>
7 :
8 : namespace tracking
9 : {
10 : namespace math
11 : {
12 :
13 : template <typename ValueType_, sint32 Rows_, sint32 Cols_, bool IsRowMajor_>
14 : class Matrix;
15 :
16 : template <typename ValueType_, sint32 Size_>
17 : class DiagonalMatrix;
18 :
19 : // TODO(matthias): add interface contract
20 : // TODO(matthias): use optimized menory storage for triangular matrixes
21 :
22 : /// \brief A triangular matrix specialization that stores only the upper or lower triangular part.
23 : ///
24 : /// This class represents triangular matrixes (upper or lower) and provides operations
25 : /// optimized for triangular structure. It inherits from SquareMatrix but restricts
26 : /// access to maintain triangular properties. Memory usage is not yet optimized and
27 : /// still stores the full square matrix internally.
28 : ///
29 : /// \tparam ValueType_ The atomic data type of internal elements
30 : /// \tparam Size_ The dimension of the triangular matrix (compile-time constant)
31 : /// \tparam IsLower_ Triangular type flag (true for lower triangular, false for upper triangular)
32 : /// \tparam IsRowMajor_ Storage layout flag (true for row-major, false for column-major)
33 : ///
34 : /// \note Currently uses full square matrix storage. Future optimization should use
35 : /// triangular storage to save memory (Size_*(Size_+1)/2 elements instead of Size_^2).
36 : ///
37 : /// \see SquareMatrix for general square matrix operations
38 : /// \see DiagonalMatrix for diagonal matrix operations
39 : template <typename ValueType_, sint32 Size_, bool IsLower_, bool IsRowMajor_ = true>
40 1390 : class TriangularMatrix TEST_REMOVE_FINAL: public SquareMatrix<ValueType_, Size_, IsRowMajor_>
41 : {
42 : public:
43 : using BaseSquareMatrix = SquareMatrix<ValueType_, Size_, IsRowMajor_>; ///< type of the parent class
44 :
45 : // unhide ctor of base class to allow implicit call in derived default ctors
46 155 : using BaseSquareMatrix::BaseSquareMatrix;
47 :
48 : /// \brief Type of the transposed matrix without changing the memory layout
49 : using transpose_type = TriangularMatrix<ValueType_, Size_, !IsLower_, !IsRowMajor_>;
50 :
51 : /// \brief Construct a new Triangular Matrix object
52 : /// \param[in] other
53 : explicit TriangularMatrix(const BaseSquareMatrix& other);
54 :
55 : /// \brief Move construct a new Triangular Matrix object
56 : /// \param[in] other
57 342 : explicit TriangularMatrix(BaseSquareMatrix&& other) noexcept
58 304 : : BaseSquareMatrix{std::move(other)} {}; // TODO(matthias): might be dangerous due to memory artifacts
59 :
60 :
61 : /// \brief Construct an identity triangular matrix.
62 : ///
63 : /// Creates a triangular matrix with ones on the diagonal and zeros elsewhere,
64 : /// maintaining the triangular structure constraint.
65 : ///
66 : /// \return TriangularMatrix An identity matrix respecting triangular constraints
67 : ///
68 : /// \note For unit triangular matrixes, this creates a matrix with ones on the diagonal
69 88 : [[nodiscard]] static auto Identity() -> TriangularMatrix { return TriangularMatrix{BaseSquareMatrix::Identity()}; }
70 :
71 : /// \brief Creates a TriangularMatrix from a nested initializer list
72 : ///
73 : /// This function constructs a TriangularMatrix from a nested initializer list.
74 : /// The triangular structure (upper or lower) is determined by the IsLower_ template parameter.
75 : ///
76 : /// \param[in] list Nested initializer list representing the triangular matrix
77 : /// \return TriangularMatrix instance initialized with the provided values
78 : /// \see TriangularFromSquare() for creating from square matrixes
79 : /// \see SquareFromList() for the underlying conversion
80 : [[nodiscard]] static auto FromList(const std::initializer_list<std::initializer_list<ValueType_>>& list) -> TriangularMatrix;
81 :
82 : /// \brief Set a triangular block within this triangular matrix.
83 : ///
84 : /// Copies a triangular block from the source matrix to the specified position
85 : /// in this matrix, maintaining triangular structure constraints.
86 : ///
87 : /// \tparam SrcSize Size of the source triangular block
88 : /// \tparam SrcCount Number of diagonal elements to copy
89 : /// \tparam SrcRowBeg Starting row index in the source block
90 : /// \tparam SrcColBeg Starting column index in the source block
91 : /// \tparam DstRowBeg Starting row index in destination (this matrix)
92 : /// \tparam DstColBeg Starting column index in destination (this matrix)
93 : /// \param[in] block The source triangular block matrix to copy from
94 : ///
95 : /// \note All indices are compile-time constants for efficiency
96 : template <sint32 SrcSize, sint32 SrcCount, sint32 SrcRowBeg, sint32 SrcColBeg, sint32 DstRowBeg, sint32 DstColBeg>
97 : void setBlock(const TriangularMatrix<ValueType_, SrcSize, IsLower_, IsRowMajor_>& block);
98 :
99 : // TODO(matthias): add setBlock with params defined at runtime
100 :
101 : /// \brief Multiply triangular matrix with a general matrix.
102 : ///
103 : /// Performs matrix multiplication T * M where T is this triangular matrix
104 : /// and M is a general matrix. Optimized for triangular structure.
105 : ///
106 : /// \tparam Cols_ Number of columns in the input matrix
107 : /// \tparam IsRowMajor2_ Storage layout of the input matrix
108 : /// \param[in] mat The matrix to multiply with (right-hand side)
109 : /// \return Matrix The result of the multiplication T * M
110 : ///
111 : /// \note Complexity is O(Size_² * Cols_), optimized for triangular operations
112 : template <sint32 Cols_, bool IsRowMajor2_>
113 : [[nodiscard]] auto operator*(const Matrix<ValueType_, Size_, Cols_, IsRowMajor2_>& mat) const
114 : -> Matrix<ValueType_, Size_, Cols_, IsRowMajor2_>;
115 :
116 : /// \brief Multiplication with triangular matrix: Tria * Matrix
117 : /// \param[in] mat A triangular matrix
118 : /// \return TriangularMatrix<value_type, Size, isLower>
119 : [[nodiscard]] auto operator*(const TriangularMatrix& mat) const -> TriangularMatrix;
120 :
121 : /// \brief Multiplication with triangular matrix: Tria * Matrix
122 : /// \param[in] mat A triangular matrix
123 : /// \return SquareMatrix<value_type, Size>
124 : [[nodiscard]] auto operator*(const TriangularMatrix<ValueType_, Size_, !IsLower_, IsRowMajor_>& mat) const -> BaseSquareMatrix;
125 :
126 : /// \brief Multiplication with diagonal matrix: Tria * Matrix
127 : /// \param[in] diag A diagonal matrix
128 : /// \return TriangularMatrix<value_type, Size, isLower>
129 : [[nodiscard]] auto operator*(const DiagonalMatrix<ValueType_, Size_>& diag) const -> TriangularMatrix;
130 :
131 : /// \brief Multiplication with scalar: Tria * scalar
132 : /// \param[in] scalar A scalar value
133 : /// \return TriangularMatrix<value_type, Size, isLower>
134 : [[nodiscard]] auto operator*(const ValueType_ scalar) const -> TriangularMatrix;
135 :
136 : /// \brief Inplace Multiplication with scalar: Tria * scalar
137 : /// \param[in] scalar A scalar value
138 : void operator*=(const ValueType_ scalar);
139 :
140 : /// \brief Element read-only access to a scalar triangular value
141 : /// \param[in] row Row index of the element
142 : /// \param[in] col Col index of the element
143 : /// \return value_type scalar triangular value
144 : [[nodiscard]] auto operator()(sint32 row, sint32 col) const -> tl::expected<ValueType_, Errors>;
145 :
146 : /// \brief Element access to a scalar triangular value
147 : /// \param[in] row Row index of the element
148 : /// \param[in] col Col index of the element
149 : /// \return value_type& Reference to the scalar triangular value
150 : [[nodiscard]] auto operator()(sint32 row, sint32 col) -> tl::expected<std::reference_wrapper<ValueType_>, Errors>;
151 :
152 : /// \brief Calculate the transposed matrix without changing the layout
153 : /// \return const transpose_type& const reference to same data as Self, but differently interpreted
154 : [[nodiscard]] auto transpose() const -> const transpose_type&;
155 :
156 : /// \brief Calculate the transposed matrix without changing the layout
157 : /// \return transpose_type& reference to same data as Self, but differently interpreted
158 : [[nodiscard]] auto transpose() -> transpose_type&;
159 :
160 : /// \brief Solve the triangular system T*x = b using forward/backward substitution.
161 : ///
162 : /// Solves the linear system T*x = b where T is this triangular matrix.
163 : /// Uses efficient forward substitution for lower triangular or backward
164 : /// substitution for upper triangular matrixes.
165 : ///
166 : /// \tparam Cols_ Number of columns in the right-hand side matrix b
167 : /// \tparam IsRowMajor2_ Storage layout of the right-hand side matrix
168 : /// \param[in] b Right-hand side matrix of the equation T*x = b
169 : /// \return Matrix Solution matrix x such that T*x ≈ b
170 : ///
171 : /// \note Complexity is O(Size_² * Cols_), much faster than general matrix solving
172 : /// \note Assumes the triangular matrix is non-singular
173 : template <sint32 Cols_, bool IsRowMajor2_>
174 : [[nodiscard]] auto solve(const Matrix<ValueType_, Size_, Cols_, IsRowMajor2_>& b) const
175 : -> Matrix<ValueType_, Size_, Cols_, IsRowMajor2_>;
176 :
177 : /// \brief Compute the inverse of the triangular matrix.
178 : ///
179 : /// Calculates the inverse using specialized triangular matrix inversion algorithms.
180 : /// The result maintains the triangular structure.
181 : ///
182 : /// \return TriangularMatrix The inverse matrix such that T * T^(-1) = I
183 : ///
184 : /// \note More efficient than general matrix inversion due to triangular structure
185 : /// \note Assumes the matrix is non-singular
186 : [[nodiscard]] auto inverse() const -> TriangularMatrix;
187 :
188 : /// \brief Calculate the determinant of the triangular matrix.
189 : ///
190 : /// Computes the determinant as the product of the diagonal elements.
191 : ///
192 : /// \return ValueType_ The determinant of the matrix
193 : /// \note Time complexity: O(n) where n is the matrix dimension
194 : /// \note For singular matrixes, the determinant will be zero or very close to zero
195 : [[nodiscard]] auto determinant() const -> ValueType_;
196 :
197 : // TODO(matthias): UnitUpper inplace inverse, Grewal Table 6.7 p.235
198 :
199 : /// \brief Checks for Unit Upper condition
200 : /// \return true
201 : [[nodiscard]] auto isUnitUpperTriangular() const -> bool;
202 :
203 : //////////////////////////////////////////////////
204 : // unsafe access operators --->
205 : /// \brief Element read-only access to a scalar triangular value
206 : /// \param[in] row Row index of the element
207 : /// \param[in] col Col index of the element
208 : /// \return value_type scalar triangular value
209 : [[nodiscard]] auto at_unsafe(sint32 row, sint32 col) const -> ValueType_;
210 :
211 : /// \brief Element access to a scalar triangular value
212 : /// \param[in] row Row index of the element
213 : /// \param[in] col Col index of the element
214 : /// \return value_type& Reference to the scalar triangular value
215 : [[nodiscard]] auto at_unsafe(sint32 row, sint32 col) -> ValueType_&;
216 : // <---
217 :
218 : private:
219 : /// \brief hide inherited transpose function
220 : using BaseSquareMatrix::transpose;
221 :
222 : /// \brief hide inherited determinant function
223 : using BaseSquareMatrix::determinant;
224 :
225 : /// \brief hide inherited operator() to prevent accessing off-triangular elements
226 : using BaseSquareMatrix::operator();
227 : };
228 :
229 : } // namespace math
230 : } // namespace tracking
231 :
232 : #endif // BC7FD90F_FBB7_481C_89C4_89BEE41309C5
|
