diff --git a/numerics/concepts.hpp b/numerics/concepts.hpp index 14ae4a180d..2c2a16f693 100644 --- a/numerics/concepts.hpp +++ b/numerics/concepts.hpp @@ -16,6 +16,8 @@ concept one_dimensional = requires(T& t, int index) { template concept two_dimensional = requires(T& t, int row, int column) { { t(row, column) } -> std::same_as; +} || requires(T const& t, int row, int column) { + { t(row, column) } -> std::same_as; }; template diff --git a/numerics/fixed_arrays.hpp b/numerics/fixed_arrays.hpp index 6a3cb5149d..813e280873 100644 --- a/numerics/fixed_arrays.hpp +++ b/numerics/fixed_arrays.hpp @@ -5,6 +5,7 @@ #include #include "base/tags.hpp" +#include "numerics/matrix_views.hpp" #include "numerics/transposed_view.hpp" #include "quantities/named_quantities.hpp" #include "quantities/si.hpp" @@ -15,6 +16,7 @@ namespace _fixed_arrays { namespace internal { using namespace principia::base::_tags; +using namespace principia::numerics::_matrix_views; using namespace principia::numerics::_transposed_view; using namespace principia::quantities::_named_quantities; using namespace principia::quantities::_si; @@ -38,6 +40,12 @@ class FixedVector final { constexpr FixedVector( std::array&& data); // NOLINT(runtime/explicit) + // Constructs a fixed vector by copying data from the view. Note that the + // result is fixed even if the matrix being viewed is an UnboundedMatrix. + template + requires std::same_as + constexpr explicit FixedVector(ColumnView const& view); + // Convertible to an array. explicit constexpr operator std::array() const; @@ -173,6 +181,8 @@ class FixedStrictlyLowerTriangularMatrix final { constexpr FixedStrictlyLowerTriangularMatrix( std::array const& data); + explicit constexpr operator FixedMatrix() const; + friend bool operator==(FixedStrictlyLowerTriangularMatrix const& left, FixedStrictlyLowerTriangularMatrix const& right) = default; @@ -215,6 +225,8 @@ class FixedLowerTriangularMatrix final { explicit FixedLowerTriangularMatrix( TransposedView> const& view); + explicit constexpr operator FixedMatrix() const; + friend bool operator==(FixedLowerTriangularMatrix const& left, FixedLowerTriangularMatrix const& right) = default; friend bool operator!=(FixedLowerTriangularMatrix const& left, @@ -252,6 +264,8 @@ class FixedUpperTriangularMatrix final { explicit FixedUpperTriangularMatrix( TransposedView> const& view); + explicit constexpr operator FixedMatrix() const; + friend bool operator==(FixedUpperTriangularMatrix const& left, FixedUpperTriangularMatrix const& right) = default; friend bool operator!=(FixedUpperTriangularMatrix const& left, diff --git a/numerics/fixed_arrays_body.hpp b/numerics/fixed_arrays_body.hpp index 0a1077b215..c873b64107 100644 --- a/numerics/fixed_arrays_body.hpp +++ b/numerics/fixed_arrays_body.hpp @@ -60,6 +60,17 @@ constexpr FixedVector::FixedVector( std::array&& data) : data_(std::move(data)) {} +template +template + requires std::same_as +constexpr FixedVector::FixedVector( + ColumnView const& view) : FixedVector(uninitialized) { + CONSTEXPR_DCHECK(view.size() == size_); + for (int i = 0; i < size_; ++i) { + (*this)[i] = view[i]; + } +} + template constexpr FixedVector::operator std::array() const { @@ -308,6 +319,18 @@ constexpr FixedStrictlyLowerTriangularMatrix:: FixedStrictlyLowerTriangularMatrix(std::array const& data) : data_(data) {} +template +constexpr FixedStrictlyLowerTriangularMatrix:: +operator FixedMatrix() const { + FixedMatrix result; // Initialized. + for (int i = 0; i < rows_; ++i) { + for (int j = 0; j < i; ++j) { + result(i, j) = (*this)(i, j); + } + } + return result; +} + template constexpr Scalar_& FixedStrictlyLowerTriangularMatrix:: operator()(int const row, int const column) { @@ -366,6 +389,18 @@ FixedLowerTriangularMatrix::FixedLowerTriangularMatrix( } } +template +constexpr FixedLowerTriangularMatrix:: +operator FixedMatrix() const { + FixedMatrix result; // Initialized. + for (int i = 0; i < rows_; ++i) { + for (int j = 0; j <= i; ++j) { + result(i, j) = (*this)(i, j); + } + } + return result; +} + template constexpr Scalar_& FixedLowerTriangularMatrix:: operator()(int const row, int const column) { @@ -417,6 +452,18 @@ FixedUpperTriangularMatrix::FixedUpperTriangularMatrix( } } +template +constexpr FixedUpperTriangularMatrix:: +operator FixedMatrix() const { + FixedMatrix result; // Initialized. + for (int j = 0; j < columns_; ++j) { + for (int i = 0; i <= j; ++i) { + result(i, j) = (*this)(i, j); + } + } + return result; +} + template constexpr Scalar_& FixedUpperTriangularMatrix:: operator()(int const row, int const column) { diff --git a/numerics/fixed_arrays_test.cpp b/numerics/fixed_arrays_test.cpp index ca4fb5dfe6..336e4e120b 100644 --- a/numerics/fixed_arrays_test.cpp +++ b/numerics/fixed_arrays_test.cpp @@ -261,6 +261,25 @@ TEST_F(FixedArraysTest, Row) { EXPECT_EQ(-24, r1 * v); } +TEST_F(FixedArraysTest, Conversions) { + using M = FixedMatrix; + EXPECT_EQ(M({0, 0, 0, 0, + 1, 0, 0, 0, + 2, 3, 0, 0, + 5, 8, 13, 0}), + M(sl4_)); + EXPECT_EQ(M({1, 0, 0, 0, + 2, 3, 0, 0, + 5, 8, 13, 0, + 21, 34, 55, 89}), + M(l4_)); + EXPECT_EQ(M({1, 2, 3, 5, + 0, 8, 13, 21, + 0, 0, 34, 55, + 0, 0, 0, 89}), + M(u4_)); +} + TEST_F(FixedArraysTest, Transpose) { EXPECT_EQ( (FixedMatrix({-8, -4, 6, diff --git a/numerics/matrix_computations.hpp b/numerics/matrix_computations.hpp index 440bb45e61..a56af77140 100644 --- a/numerics/matrix_computations.hpp +++ b/numerics/matrix_computations.hpp @@ -29,6 +29,14 @@ struct ᵗRDRDecompositionGenerator; template struct SubstitutionGenerator; +// Declares: +// struct Result { +// ⟨upper triangular matrix⟩ R; +// ⟨matrix⟩ Q; +// }; +template +struct GramSchmidtGenerator; + // Declares: // struct Result { // (matrix) H; @@ -99,6 +107,10 @@ typename SubstitutionGenerator::Result ForwardSubstitution(LowerTriangularMatrix const& L, Vector const& b); +template +typename GramSchmidtGenerator::Result +ClassicalGramSchmidt(Matrix const& L); + // If A is a square matrix, returns U and H so that A = ᵗU H U, where H is an // upper Hessenberg matrix. // TODO(phl): Add support for returning U. @@ -141,6 +153,7 @@ Solve(Matrix A, Vector b); using internal::BackSubstitution; using internal::CholeskyDecomposition; +using internal::ClassicalGramSchmidt; using internal::ClassicalJacobi; using internal::ForwardSubstitution; using internal::HessenbergDecomposition; diff --git a/numerics/matrix_computations_body.hpp b/numerics/matrix_computations_body.hpp index 826e59fa57..6f6eb4124d 100644 --- a/numerics/matrix_computations_body.hpp +++ b/numerics/matrix_computations_body.hpp @@ -285,6 +285,29 @@ struct SubstitutionGenerator, static Result Uninitialized(TriangularMatrix const& m); }; +template +struct GramSchmidtGenerator> { + struct Result { + UnboundedMatrix Q; + UnboundedUpperTriangularMatrix R; + }; + using AVector = UnboundedVector; + using QVector = UnboundedVector; + static Result Uninitialized(UnboundedMatrix const& m); +}; + +template +struct GramSchmidtGenerator> { + struct Result { + FixedMatrix Q; + FixedUpperTriangularMatrix R; + }; + using AVector = FixedVector; + using QVector = FixedVector; + static Result Uninitialized( + FixedMatrix const& m); +}; + template struct HessenbergDecompositionGenerator> { struct Result { @@ -432,6 +455,22 @@ Uninitialized(TriangularMatrix const& m) -> Result { return Result(m.columns(), uninitialized); } +template +auto GramSchmidtGenerator>:: +Uninitialized(UnboundedMatrix const& m) -> Result { + return Result{ + .Q = UnboundedMatrix(m.rows(), m.columns(), uninitialized), + .R = UnboundedUpperTriangularMatrix(m.columns(), uninitialized)}; +} + +template +auto GramSchmidtGenerator>:: +Uninitialized(FixedMatrix const& m) -> Result { + return Result{ + .Q = FixedMatrix(uninitialized), + .R = FixedUpperTriangularMatrix(uninitialized)}; +} + template auto ClassicalJacobiGenerator>::Identity( UnboundedMatrix const& m) -> Rotation { @@ -603,6 +642,44 @@ ForwardSubstitution(LowerTriangularMatrix const& L, return x; } +template +typename GramSchmidtGenerator::Result ClassicalGramSchmidt( + Matrix const& A) { + using G = GramSchmidtGenerator; + auto result = G::Uninitialized(A); + auto& Q = result.Q; + auto& R = result.R; + int const n = A.rows(); + + // [Hig02], Algorithm 19.11. + for (int j = 0; j < n; ++j) { + auto const aⱼ = ColumnView{.matrix = A, + .first_row = 0, + .last_row = n - 1, + .column = j}; + for (int i = 0; i < j; ++i) { + auto const qᵢ = ColumnView{.matrix = Q, + .first_row = 0, + .last_row = n - 1, + .column = i}; + R(i, j) = TransposedView{.transpose = qᵢ} * aⱼ; // NOLINT + } + typename G::AVector qʹⱼ(aⱼ); + for (int k = 0; k < j; ++k) { + qʹⱼ -= R(k, j) * typename G::QVector(ColumnView{.matrix = Q, + .first_row = 0, + .last_row = n - 1, + .column = k}); + } + R(j, j) = qʹⱼ.Norm(); + ColumnView{.matrix = Q, // NOLINT + .first_row = 0, + .last_row = n - 1, + .column = j} = qʹⱼ / R(j, j); + } + return result; +} + template typename HessenbergDecompositionGenerator::Result HessenbergDecomposition(Matrix const& A) { diff --git a/numerics/matrix_computations_test.cpp b/numerics/matrix_computations_test.cpp index 3da88137f9..6229c3fb45 100644 --- a/numerics/matrix_computations_test.cpp +++ b/numerics/matrix_computations_test.cpp @@ -5,6 +5,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" #include "numerics/fixed_arrays.hpp" +#include "numerics/transposed_view.hpp" #include "numerics/unbounded_arrays.hpp" #include "quantities/elementary_functions.hpp" #include "testing_utilities/almost_equals.hpp" @@ -18,6 +19,7 @@ namespace numerics { using namespace principia::numerics::_fixed_arrays; using namespace principia::numerics::_matrix_computations; +using namespace principia::numerics::_transposed_view; using namespace principia::numerics::_unbounded_arrays; using namespace principia::quantities::_elementary_functions; using namespace principia::testing_utilities::_almost_equals; @@ -110,6 +112,31 @@ TYPED_TEST(MatrixComputationsTest, ForwardSubstitution) { EXPECT_THAT(x4_actual, AlmostEquals(x4_expected, 0)); } +TYPED_TEST(MatrixComputationsTest, ClassicalGramSchmidt) { + using Matrix = typename std::tuple_element<3, TypeParam>::type; + Matrix const m4({1, 2, 3, -4, + 5, 6, 7, 8, + 9, 8, -7, 6, + 5, 4, 3, 2}); + auto const qr = ClassicalGramSchmidt(m4); + + // Check that the decomposition is correct. + auto const near_m4 = qr.Q * Matrix(qr.R); + EXPECT_THAT(near_m4, AlmostEquals(m4, 1)); + + // Check that Q is nearly orthogonal. + auto const near_identity = Matrix(TransposedView{.transpose = qr.Q}) * qr.Q; + for (int i = 0; i < near_identity.rows(); ++i) { + for (int j = 0; j < near_identity.columns(); ++j) { + if (i == j) { + EXPECT_THAT(near_identity(i, j), AlmostEquals(1, 0, 2)); + } else { + EXPECT_THAT(near_identity(i, j), VanishesBefore(1, 0, 3)); + } + } + } +} + TYPED_TEST(MatrixComputationsTest, HessenbergForm) { using Matrix = typename std::tuple_element<3, TypeParam>::type; Matrix const m4({1, 2, 3, -4, diff --git a/numerics/matrix_views.hpp b/numerics/matrix_views.hpp index ca9229e3b5..a1fbec1e41 100644 --- a/numerics/matrix_views.hpp +++ b/numerics/matrix_views.hpp @@ -31,6 +31,13 @@ struct BlockView { constexpr Scalar& operator()(int row, int column); constexpr Scalar const& operator()(int row, int column) const; + template + requires two_dimensional && same_elements_as + BlockView& operator=(T const& right); + template + requires two_dimensional && same_elements_as + BlockView& operator=(T&& right); + template requires two_dimensional && same_elements_as BlockView& operator+=(T const& right); @@ -59,6 +66,13 @@ struct ColumnView { constexpr Scalar& operator[](int index); constexpr Scalar const& operator[](int index) const; + template + requires one_dimensional && same_elements_as + ColumnView& operator=(T const& right); + template + requires one_dimensional && same_elements_as + ColumnView& operator=(T&& right); + template requires one_dimensional && same_elements_as ColumnView& operator+=(T const& right); @@ -74,6 +88,11 @@ struct ColumnView { constexpr int size() const; }; +template +Product operator*( + TransposedView> const& left, + ColumnView const& right); + // The declarations below also declare the symmetric operator== and // operator!=. diff --git a/numerics/matrix_views_body.hpp b/numerics/matrix_views_body.hpp index 1290b94ada..613d9936dc 100644 --- a/numerics/matrix_views_body.hpp +++ b/numerics/matrix_views_body.hpp @@ -2,6 +2,8 @@ #include "numerics/matrix_views.hpp" +#include + #include "base/tags.hpp" #include "quantities/elementary_functions.hpp" @@ -32,6 +34,36 @@ constexpr auto BlockView::operator()( return matrix(first_row + row, first_column + column); } +template + requires two_dimensional +template + requires two_dimensional && same_elements_as +BlockView& BlockView::operator=(T const& right) { + DCHECK_EQ(rows(), right.rows()); + DCHECK_EQ(columns(), right.columns()); + for (int i = 0; i < rows(); ++i) { + for (int j = 0; j < columns(); ++j) { + matrix(first_row + i, first_column + j) = right(i, j); + } + } + return *this; +} + +template + requires two_dimensional +template + requires two_dimensional && same_elements_as +BlockView& BlockView::operator=(T&& right) { + DCHECK_EQ(rows(), right.rows()); + DCHECK_EQ(columns(), right.columns()); + for (int i = 0; i < rows(); ++i) { + for (int j = 0; j < columns(); ++j) { + matrix(first_row + i, first_column + j) = std::move(right(i, j)); + } + } + return *this; +} + template requires two_dimensional template @@ -112,6 +144,30 @@ constexpr auto ColumnView::operator[]( return matrix(first_row + index, column); } +template + requires two_dimensional +template + requires one_dimensional && same_elements_as +ColumnView& ColumnView::operator=(T const& right) { + DCHECK_EQ(size(), right.size()); + for (int i = 0; i < size(); ++i) { + matrix(first_row + i, column) = right[i]; + } + return *this; +} + +template + requires two_dimensional +template + requires one_dimensional && same_elements_as +ColumnView& ColumnView::operator=(T&& right) { + DCHECK_EQ(size(), right.size()); + for (int i = 0; i < size(); ++i) { + matrix(first_row + i, column) = std::move(right[i]); + } + return *this; +} + template requires two_dimensional template @@ -176,6 +232,18 @@ constexpr auto ColumnView::size() const -> int { return last_row - first_row + 1; } +template +Product operator*( + TransposedView> const& left, + ColumnView const& right) { + DCHECK_EQ(left.size(), right.size()); + Product result{}; + for (int i = 0; i < left.size(); ++i) { + result += left[i] * right[i]; + } + return result; +} + template requires two_dimensional && same_elements_as bool operator==(BlockView const& left, T const& right) { diff --git a/numerics/matrix_views_test.cpp b/numerics/matrix_views_test.cpp index feb8f3d22b..d2d3fafdae 100644 --- a/numerics/matrix_views_test.cpp +++ b/numerics/matrix_views_test.cpp @@ -2,6 +2,7 @@ #include "gtest/gtest.h" #include "numerics/fixed_arrays.hpp" +#include "numerics/transposed_view.hpp" #include "numerics/unbounded_arrays.hpp" #include "quantities/elementary_functions.hpp" #include "testing_utilities/almost_equals.hpp" @@ -11,6 +12,7 @@ namespace numerics { using namespace principia::numerics::_fixed_arrays; using namespace principia::numerics::_matrix_views; +using namespace principia::numerics::_transposed_view; using namespace principia::numerics::_unbounded_arrays; using namespace principia::quantities::_elementary_functions; using namespace principia::testing_utilities::_almost_equals; @@ -188,5 +190,18 @@ TEST_F(MatrixViewsTest, ColumnView_Norm) { EXPECT_THAT(cum34.Norm(), AlmostEquals(Sqrt(106), 0)); } +TEST_F(MatrixViewsTest, ColumnView_DotProduct) { + EXPECT_EQ( + -60, + (TransposedView{.transpose = ColumnView{.matrix = fm34_, // NOLINT + .first_row = 0, + .last_row = 2, + .column = 1}} * + ColumnView{.matrix = um34_, + .first_row = 0, + .last_row = 2, + .column = 2})); +} + } // namespace numerics } // namespace principia diff --git a/numerics/transposed_view_body.hpp b/numerics/transposed_view_body.hpp index 659adf1b4c..46ffb1a1fa 100644 --- a/numerics/transposed_view_body.hpp +++ b/numerics/transposed_view_body.hpp @@ -28,20 +28,21 @@ constexpr int TransposedView::size() const } template -constexpr typename T::Scalar& TransposedView::operator[](int index) +constexpr typename T::Scalar& TransposedView::operator[](int const index) requires one_dimensional { return transpose[index]; } template constexpr typename T::Scalar const& TransposedView::operator[]( - int index) const + int const index) const requires one_dimensional { return transpose[index]; } template -constexpr typename T::Scalar& TransposedView::operator()(int row, int column) +constexpr typename T::Scalar& TransposedView::operator()(int const row, + int const column) requires two_dimensional { // Note the transposition. return transpose(column, row); @@ -49,8 +50,8 @@ constexpr typename T::Scalar& TransposedView::operator()(int row, int column) template constexpr typename T::Scalar const& TransposedView::operator()( - int row, - int column) const + int const row, + int const column) const requires two_dimensional { // Note the transposition. return transpose(column, row); diff --git a/numerics/unbounded_arrays.hpp b/numerics/unbounded_arrays.hpp index e8c9b1d3ae..ab6addbde9 100644 --- a/numerics/unbounded_arrays.hpp +++ b/numerics/unbounded_arrays.hpp @@ -178,6 +178,8 @@ class UnboundedLowerTriangularMatrix final { explicit UnboundedLowerTriangularMatrix( TransposedView> const& view); + explicit operator UnboundedMatrix() const; + friend bool operator==(UnboundedLowerTriangularMatrix const& left, UnboundedLowerTriangularMatrix const& right) = default; friend bool operator!=(UnboundedLowerTriangularMatrix const& left, @@ -227,6 +229,8 @@ class UnboundedUpperTriangularMatrix final { explicit UnboundedUpperTriangularMatrix( TransposedView> const& view); + explicit operator UnboundedMatrix() const; + friend bool operator==(UnboundedUpperTriangularMatrix const& left, UnboundedUpperTriangularMatrix const& right) = default; friend bool operator!=(UnboundedUpperTriangularMatrix const& left, diff --git a/numerics/unbounded_arrays_body.hpp b/numerics/unbounded_arrays_body.hpp index f8187de7d2..6b668cff91 100644 --- a/numerics/unbounded_arrays_body.hpp +++ b/numerics/unbounded_arrays_body.hpp @@ -355,6 +355,18 @@ UnboundedLowerTriangularMatrix::UnboundedLowerTriangularMatrix( } } +template +UnboundedLowerTriangularMatrix::operator UnboundedMatrix< + Scalar_>() const { + UnboundedMatrix result(rows_, rows_); // Initialized. + for (int i = 0; i < rows(); ++i) { + for (int j = 0; j <= i; ++j) { + result(i, j) = (*this)(i, j); + } + } + return result; +} + template Scalar_& UnboundedLowerTriangularMatrix::operator()( int const row, int const column) { @@ -456,6 +468,18 @@ UnboundedUpperTriangularMatrix::UnboundedUpperTriangularMatrix( } } +template +UnboundedUpperTriangularMatrix:: +operator UnboundedMatrix() const { + UnboundedMatrix result(columns_, columns_); // Initialized. + for (int j = 0; j < columns_; ++j) { + for (int i = 0; i <= j; ++i) { + result(i, j) = (*this)(i, j); + } + } + return result; +} + template Scalar_& UnboundedUpperTriangularMatrix::operator()( int const row, int const column) { diff --git a/numerics/unbounded_arrays_test.cpp b/numerics/unbounded_arrays_test.cpp index efe4e35f1b..3534adb146 100644 --- a/numerics/unbounded_arrays_test.cpp +++ b/numerics/unbounded_arrays_test.cpp @@ -228,6 +228,21 @@ TEST_F(UnboundedArraysTest, UpperTriangularMatrixIndexing) { EXPECT_EQ(1, u4(0, 0)); } +TEST_F(UnboundedArraysTest, Conversions) { + EXPECT_EQ(UnboundedMatrix(4, 4, + {1, 0, 0, 0, + 2, 3, 0, 0, + 5, 8, 13, 0, + 21, 34, 55, 89}), + UnboundedMatrix(l4_)); + EXPECT_EQ(UnboundedMatrix(4, 4, + {1, 2, 3, 5, + 0, 8, 13, 21, + 0, 0, 34, 55, + 0, 0, 0, 89}), + UnboundedMatrix(u4_)); +} + TEST_F(UnboundedArraysTest, Extend) { { UnboundedVector u2({1, 2}); diff --git a/testing_utilities/approximate_quantity_body.hpp b/testing_utilities/approximate_quantity_body.hpp index c814ae2823..56038c0f03 100644 --- a/testing_utilities/approximate_quantity_body.hpp +++ b/testing_utilities/approximate_quantity_body.hpp @@ -101,10 +101,10 @@ inline ApproximateQuantity ApproximateQuantity::Parse( } } if (ulp <= 9) { - error_representation[*last_digit_index] = '0' + ulp; + error_representation[*last_digit_index] = static_cast('0' + ulp); } else { CHECK(is_hexadecimal); - error_representation[*last_digit_index] = 'A' + ulp - 10; + error_representation[*last_digit_index] = static_cast('A' + ulp - 10); } // Apparently the stupid language doesn't know how to parse literals with