diff --git a/numerics/fixed_arrays.hpp b/numerics/fixed_arrays.hpp index 4eed4af7e4..2072557359 100644 --- a/numerics/fixed_arrays.hpp +++ b/numerics/fixed_arrays.hpp @@ -39,8 +39,6 @@ class FixedVector final { constexpr FixedVector( std::array&& data); // NOLINT(runtime/explicit) - TransposedView Transpose() const; - Scalar Norm() const; Square Norm²() const; @@ -104,11 +102,14 @@ class FixedMatrix final { Scalar const* row() const; FixedMatrix Transpose() const; + Scalar FrobeniusNorm() const; bool operator==(FixedMatrix const& right) const; bool operator!=(FixedMatrix const& right) const; + // Applies the matrix as a bilinear form. Present for compatibility with + // |SymmetricBilinearForm|. Prefer to use |TransposedView| and |operator*|. template Product> operator()(FixedVector const& left, @@ -214,6 +215,7 @@ class FixedUpperTriangularMatrix final { std::array data_; }; +// Prefer using the operator* that takes a TransposedView. template constexpr Product InnerProduct( FixedVector const& left, @@ -363,7 +365,8 @@ constexpr FixedVector, rows> operator*( FixedMatrix const& left, FixedVector const& right); -// Use this operator to multiply a row vector with a matrix. +// Use this operator to multiply a row vector with a matrix. We don't have an +// operator returning a TransposedView as that would cause dangling references. template constexpr FixedVector, columns> operator*( TransposedView> const& left, @@ -379,6 +382,11 @@ template std::ostream& operator<<(std::ostream& out, FixedMatrix const& matrix); +template +std::ostream& operator<<( + std::ostream& out, + FixedStrictlyLowerTriangularMatrix const& matrix); + template std::ostream& operator<<( std::ostream& out, diff --git a/numerics/fixed_arrays_body.hpp b/numerics/fixed_arrays_body.hpp index f519c2a574..adeaede9a9 100644 --- a/numerics/fixed_arrays_body.hpp +++ b/numerics/fixed_arrays_body.hpp @@ -71,12 +71,6 @@ constexpr FixedVector::FixedVector( std::array&& data) : data_(std::move(data)) {} -template -TransposedView> -FixedVector::Transpose() const { - return {.transpose = *this}; -} - template Scalar FixedVector::Norm() const { return Sqrt(Norm²()); @@ -167,15 +161,6 @@ Scalar const* FixedMatrix::row() const { return &data_[r * columns()]; } -template -template -Product> -FixedMatrix::operator()( - FixedVector const& left, - FixedVector const& right) const { - return left.Transpose() * (*this * right); -} - template FixedMatrix FixedMatrix::Transpose() const { @@ -211,6 +196,15 @@ bool FixedMatrix::operator!=( return data_ != right.data_; } +template +template +Product> +FixedMatrix::operator()( + FixedVector const& left, + FixedVector const& right) const { + return TransposedView{left} * (*this * right); // NOLINT +} + template FixedMatrix FixedMatrix::Identity() { @@ -435,7 +429,7 @@ constexpr FixedVector operator-( for (int i = 0; i < size; ++i) { result[i] = -right[i]; } - return FixedVector, size>(std::move(result)); + return FixedVector(std::move(result)); } template @@ -733,6 +727,24 @@ std::ostream& operator<<(std::ostream& out, return out; } +template +std::ostream& operator<<( + std::ostream& out, + FixedStrictlyLowerTriangularMatrix const& matrix) { + out << "rows: " << matrix.rows() << "\n"; + for (int i = 0; i < matrix.rows(); ++i) { + out << "{"; + for (int j = 0; j < i; ++j) { + out << matrix(i, j); + if (j < i - 1) { + out << ", "; + } + } + out << "}\n"; + } + return out; +} + template std::ostream& operator<<( std::ostream& out, diff --git a/numerics/fixed_arrays_test.cpp b/numerics/fixed_arrays_test.cpp index 92dc97a2d7..c60d26cbec 100644 --- a/numerics/fixed_arrays_test.cpp +++ b/numerics/fixed_arrays_test.cpp @@ -86,7 +86,7 @@ TEST_F(FixedArraysTest, Assignment) { } TEST_F(FixedArraysTest, Norm) { - EXPECT_EQ(35, v4_.Transpose() * v4_); + EXPECT_EQ(35, TransposedView{v4_} * v4_); // NOLINT EXPECT_EQ(Sqrt(35.0), v4_.Norm()); EXPECT_EQ(35, v4_.Norm²()); EXPECT_EQ(Sqrt(517.0), m34_.FrobeniusNorm()); @@ -121,11 +121,11 @@ TEST_F(FixedArraysTest, VectorSpaces) { } TEST_F(FixedArraysTest, Algebra) { - EXPECT_EQ(-535, u3_.Transpose() * v3_); + EXPECT_EQ(-535, TransposedView{u3_} * v3_); // NOLINT EXPECT_EQ((FixedMatrix({-30, -30, 10, 40, -93, -93, 31, 124, 141, 141, -47, -188})), - v3_ * v4_.Transpose()); + v3_ * TransposedView{v4_}); EXPECT_EQ((FixedMatrix({ 0, 14, -22, 3, 14, -63, 5, -92})), m23_ * m34_); diff --git a/numerics/matrix_computations_body.hpp b/numerics/matrix_computations_body.hpp index 3f26468591..4692a44027 100644 --- a/numerics/matrix_computations_body.hpp +++ b/numerics/matrix_computations_body.hpp @@ -924,7 +924,7 @@ template typename RayleighQuotientGenerator::Result RayleighQuotient(Matrix const& A, Vector const& x) { // [GV13], section 8.2.3. - return x.Transpose() * (A * x) / (x.Transpose() * x); + return TransposedView{x} * (A * x) / (TransposedView{x} * x); // NOLINT } template diff --git a/numerics/unbounded_arrays.hpp b/numerics/unbounded_arrays.hpp index 5f7f297492..64cc6bd81e 100644 --- a/numerics/unbounded_arrays.hpp +++ b/numerics/unbounded_arrays.hpp @@ -3,12 +3,14 @@ #include #include #include +#include #include #include "base/tags.hpp" #include "numerics/fixed_arrays.hpp" #include "numerics/transposed_view.hpp" #include "quantities/named_quantities.hpp" +#include "quantities/si.hpp" namespace principia { namespace numerics { @@ -21,6 +23,7 @@ using namespace principia::base::_tags; using namespace principia::numerics::_fixed_arrays; using namespace principia::numerics::_transposed_view; using namespace principia::quantities::_named_quantities; +using namespace principia::quantities::_si; // An allocator that does not initialize the allocated objects. template @@ -47,8 +50,6 @@ class UnboundedVector final { template explicit UnboundedVector(FixedVector const& data); - TransposedView Transpose() const; - void Extend(int extra_size); void Extend(int extra_size, uninitialized_t); void Extend(std::initializer_list data); @@ -60,12 +61,23 @@ class UnboundedVector final { int size() const; + typename std::vector::const_iterator begin() const; + typename std::vector::const_iterator end() const; + Scalar& operator[](int index); Scalar const& operator[](int index) const; bool operator==(UnboundedVector const& right) const; bool operator!=(UnboundedVector const& right) const; + template + friend H AbslHashValue(H h, UnboundedVector const& vector) { + for (auto const& scalar : vector) { + h = H::combine(std::move(h), scalar / si::Unit); + } + return h; + } + private: std::vector> data_; }; @@ -76,11 +88,13 @@ class UnboundedMatrix final { UnboundedMatrix(int rows, int columns); UnboundedMatrix(int rows, int columns, uninitialized_t); - // The |data| must be in row-major format. + // The |data| must be in row-major format and must be for a square matrix. UnboundedMatrix(std::initializer_list data); - int columns() const; + UnboundedMatrix(int rows, int columns, std::initializer_list data); + int rows() const; + int columns() const; int size() const; // For 0 ≤ i < rows and 0 ≤ j < columns, the entry a_ij is accessed as @@ -90,11 +104,19 @@ class UnboundedMatrix final { Scalar const& operator()(int row, int column) const; UnboundedMatrix Transpose() const; + Scalar FrobeniusNorm() const; bool operator==(UnboundedMatrix const& right) const; bool operator!=(UnboundedMatrix const& right) const; + // Applies the matrix as a bilinear form. Present for compatibility with + // |SymmetricBilinearForm|. Prefer to use |TransposedView| and |operator*|. + template + Product> + operator()(UnboundedVector const& left, + UnboundedVector const& right) const; + static UnboundedMatrix Identity(int rows, int columns); private: @@ -124,8 +146,8 @@ class UnboundedLowerTriangularMatrix final { void EraseToEnd(int begin_row_index); - int columns() const; int rows() const; + int columns() const; int size() const; // For 0 ≤ j ≤ i < rows, the entry a_ij is accessed as |a(i, j)|. @@ -166,8 +188,8 @@ class UnboundedUpperTriangularMatrix final { void EraseToEnd(int begin_column_index); - int columns() const; int rows() const; + int columns() const; int size() const; // For 0 ≤ i ≤ j < columns, the entry a_ij is accessed as |a(i, j)|. @@ -203,28 +225,71 @@ class UnboundedUpperTriangularMatrix final { friend class Row; }; +// Prefer using the operator* that takes a TransposedView. +template +Product InnerProduct( + UnboundedVector const& left, + UnboundedVector const& right); + template UnboundedVector Normalize(UnboundedVector const& vector); +template +UnboundedMatrix> SymmetricProduct( + UnboundedVector const& left, + UnboundedVector const& right); + +template +UnboundedMatrix> SymmetricSquare( + UnboundedVector const& vector); + // Additive groups. template -constexpr UnboundedVector& operator+=( +UnboundedVector operator-( + UnboundedVector const& right); + +template +UnboundedMatrix operator-( + UnboundedMatrix const& right); + +template +UnboundedVector> operator+( + UnboundedVector const& left, + UnboundedVector const& right); + +template +UnboundedMatrix> operator+( + UnboundedMatrix const& left, + UnboundedMatrix const& right); + +template +UnboundedVector> operator-( + UnboundedVector const& left, + UnboundedVector const& right); + +template +UnboundedMatrix> operator-( + UnboundedMatrix const& left, + UnboundedMatrix const& right); + +template +UnboundedVector& operator+=( UnboundedVector& left, UnboundedVector const& right); template -constexpr UnboundedMatrix& operator+=( +UnboundedMatrix& operator+=( UnboundedMatrix& left, UnboundedMatrix const& right); template -constexpr UnboundedVector& operator-=( +UnboundedVector& operator-=( UnboundedVector& left, UnboundedVector const& right); template -constexpr UnboundedMatrix& operator-=( +UnboundedMatrix& operator-=( UnboundedMatrix& left, UnboundedMatrix const& right); @@ -256,32 +321,34 @@ UnboundedVector> operator/( RScalar const& right); template -constexpr UnboundedMatrix> +UnboundedMatrix> operator/(UnboundedMatrix const& left, RScalar const& right); template -constexpr UnboundedVector& operator*=( +UnboundedVector& operator*=( UnboundedVector& left, double right); template -constexpr UnboundedMatrix& operator*=( +UnboundedMatrix& operator*=( UnboundedMatrix& left, double right); template -constexpr UnboundedVector& operator/=( +UnboundedVector& operator/=( UnboundedVector& left, double right); template -constexpr UnboundedMatrix& operator/=( +UnboundedMatrix& operator/=( UnboundedMatrix& left, double right); // Hilbert space and algebra. +// TODO(phl): fixed_arrays.hpp has an operator* that takes a row. + template Product operator*( TransposedView> const& left, @@ -302,7 +369,8 @@ UnboundedVector> operator*( UnboundedMatrix const& left, UnboundedVector const& right); -// Use this operator to multiply a row vector with a matrix. +// Use this operator to multiply a row vector with a matrix. We don't have an +// operator returning a TransposedView as that would cause dangling references. template UnboundedVector> operator*( TransposedView> const& left, diff --git a/numerics/unbounded_arrays_body.hpp b/numerics/unbounded_arrays_body.hpp index e369ec864e..444ad3fc25 100644 --- a/numerics/unbounded_arrays_body.hpp +++ b/numerics/unbounded_arrays_body.hpp @@ -38,12 +38,6 @@ template UnboundedVector::UnboundedVector(FixedVector const& data) : data_(data.begin(), data.end()) {} -template -TransposedView> -UnboundedVector::Transpose() const { - return {.transpose = *this}; -} - template void UnboundedVector::Extend(int const extra_size) { DCHECK_LE(0, extra_size); @@ -85,6 +79,18 @@ int UnboundedVector::size() const { return data_.size(); } +template +typename std::vector::const_iterator UnboundedVector::begin() + const { + return data_.cbegin(); +} + +template +typename std::vector::const_iterator UnboundedVector::end() + const { + return data_.cend(); +} + template Scalar& UnboundedVector::operator[](int const index) { DCHECK_LE(0, index); @@ -131,8 +137,12 @@ UnboundedMatrix::UnboundedMatrix(std::initializer_list data) } template -int UnboundedMatrix::columns() const { - return columns_; +UnboundedMatrix::UnboundedMatrix(int const rows, int const columns, + std::initializer_list data) + : rows_(rows), + columns_(columns), + data_(std::move(data)) { + CHECK_EQ(data.size(), rows_ * columns_); } template @@ -140,6 +150,11 @@ int UnboundedMatrix::rows() const { return rows_; } +template +int UnboundedMatrix::columns() const { + return columns_; +} + template int UnboundedMatrix::size() const { return rows_ * columns_; @@ -197,6 +212,15 @@ bool UnboundedMatrix::operator!=(UnboundedMatrix const& right) const { return !(*this == right); } +template +template +Product> +UnboundedMatrix::operator()( + UnboundedVector const& left, + UnboundedVector const& right) const { + return TransposedView{left} * (*this * right); // NOLINT +} + template UnboundedMatrix UnboundedMatrix::Identity(int const rows, int const columns) { @@ -258,12 +282,12 @@ void UnboundedLowerTriangularMatrix::EraseToEnd( } template -int UnboundedLowerTriangularMatrix::columns() const { +int UnboundedLowerTriangularMatrix::rows() const { return rows_; } template -int UnboundedLowerTriangularMatrix::rows() const { +int UnboundedLowerTriangularMatrix::columns() const { return rows_; } @@ -375,12 +399,12 @@ void UnboundedUpperTriangularMatrix::EraseToEnd( } template -int UnboundedUpperTriangularMatrix::columns() const { +int UnboundedUpperTriangularMatrix::rows() const { return columns_; } template -int UnboundedUpperTriangularMatrix::rows() const { +int UnboundedUpperTriangularMatrix::columns() const { return columns_; } @@ -478,34 +502,153 @@ UnboundedUpperTriangularMatrix::Transpose( return result; } +template +Product InnerProduct(UnboundedVector const& left, + UnboundedVector const& right) { + return TransposedView{left} * right; // NOLINT +} + template UnboundedVector Normalize(UnboundedVector const& vector) { return vector / vector.Norm(); } +template +UnboundedMatrix> SymmetricProduct( + UnboundedVector const& left, + UnboundedVector const& right) { + CHECK_EQ(left.size(), right.size()); + UnboundedMatrix> result( + left.size(), right.size(), uninitialized); + for (int i = 0; i < left.size(); ++i) { + for (int j = 0; j < i; ++j) { + auto const r = 0.5 * (left[i] * right[j] + left[j] * right[i]); + result(i, j) = r; + result(j, i) = r; + } + result(i, i) = left[i] * right[i]; + } + return result; +} + +template +UnboundedMatrix> SymmetricSquare( + UnboundedVector const& vector) { + UnboundedMatrix> result( + vector.size(), vector.size(), uninitialized); + for (int i = 0; i < vector.size(); ++i) { + for (int j = 0; j < i; ++j) { + auto const r = vector[i] * vector[j]; + result(i, j) = r; + result(j, i) = r; + } + result(i, i) = Pow<2>(vector[i]); + } + return result; +} + +template +UnboundedVector operator-(UnboundedVector const& right) { + std::vector result; + result.reserve(right.size()); + for (int i = 0; i < right.size(); ++i) { + result[i] = -right[i]; + } + return UnboundedVector(std::move(result)); +} + +template +UnboundedMatrix operator-(UnboundedMatrix const& right) { + UnboundedMatrix result(right.rows(), right.columns(), uninitialized); + for (int i = 0; i < right.rows(); ++i) { + for (int j = 0; j < right.columns(); ++j) { + result(i, j) = -right(i, j); + } + } + return result; +} + +template +UnboundedVector> operator+( + UnboundedVector const& left, + UnboundedVector const& right) { + CHECK_EQ(left.size(), right.size()); + std::vector> result; + result.resize(right.size()); + for (int i = 0; i < right.size(); ++i) { + result[i] = left[i] + right[i]; + } + return UnboundedVector>(std::move(result)); +} + +template +UnboundedMatrix> operator+( + UnboundedMatrix const& left, + UnboundedMatrix const& right) { + CHECK_EQ(left.rows(), right.rows()); + CHECK_EQ(left.columns(), right.columns()); + UnboundedMatrix> result( + right.rows(), right.columns(), uninitialized); + for (int i = 0; i < right.rows(); ++i) { + for (int j = 0; j < right.columns(); ++j) { + result(i, j) = left(i, j) + right(i, j); + } + } + return result; +} + +template +UnboundedVector> operator-( + UnboundedVector const& left, + UnboundedVector const& right) { + CHECK_EQ(left.size(), right.size()); + std::vector> result; + result.resize(right.size()); + for (int i = 0; i < right.size(); ++i) { + result[i] = left[i] - right[i]; + } + return UnboundedVector>(std::move(result)); +} + +template +UnboundedMatrix> operator-( + UnboundedMatrix const& left, + UnboundedMatrix const& right) { + CHECK_EQ(left.rows(), right.rows()); + CHECK_EQ(left.columns(), right.columns()); + UnboundedMatrix> result( + right.rows(), right.columns(), uninitialized); + for (int i = 0; i < right.rows(); ++i) { + for (int j = 0; j < right.columns(); ++j) { + result(i, j) = left(i, j) + right(i, j); + } + } + return result; +} + template -constexpr UnboundedVector& operator+=( +UnboundedVector& operator+=( UnboundedVector& left, UnboundedVector const& right) { return left = left + right; } template -constexpr UnboundedMatrix& operator+=( +UnboundedMatrix& operator+=( UnboundedMatrix& left, UnboundedMatrix const& right) { return left = left + right; } template -constexpr UnboundedVector& operator-=( +UnboundedVector& operator-=( UnboundedVector& left, UnboundedVector const& right) { return left = left - right; } template -constexpr UnboundedMatrix& operator-=( +UnboundedMatrix& operator-=( UnboundedMatrix& left, UnboundedMatrix const& right) { return left = left - right; @@ -593,26 +736,26 @@ UnboundedMatrix> operator/( } template -constexpr UnboundedVector& operator*=(UnboundedVector& left, - double const right) { +UnboundedVector& operator*=(UnboundedVector& left, + double const right) { return left = left * right; } template -constexpr UnboundedMatrix& operator*=(UnboundedMatrix& left, - double const right) { +UnboundedMatrix& operator*=(UnboundedMatrix& left, + double const right) { return left = left * right; } template -constexpr UnboundedVector& operator/=(UnboundedVector& left, - double const right) { +UnboundedVector& operator/=(UnboundedVector& left, + double const right) { return left = left / right; } template -constexpr UnboundedMatrix& operator/=(UnboundedMatrix& left, - double const right) { +UnboundedMatrix& operator/=(UnboundedMatrix& left, + double const right) { return left = left / right; } diff --git a/numerics/unbounded_arrays_test.cpp b/numerics/unbounded_arrays_test.cpp index 6811ed6c1d..91616dc597 100644 --- a/numerics/unbounded_arrays_test.cpp +++ b/numerics/unbounded_arrays_test.cpp @@ -2,12 +2,14 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include "numerics/transposed_view.hpp" #include "quantities/elementary_functions.hpp" #include "testing_utilities/almost_equals.hpp" namespace principia { namespace numerics { +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; @@ -93,7 +95,7 @@ TEST_F(UnboundedArraysTest, Assignment) { } TEST_F(UnboundedArraysTest, Norm) { - EXPECT_EQ(35, v4_.Transpose() * v4_); + EXPECT_EQ(35, TransposedView{v4_} * v4_); // NOLINT EXPECT_EQ(Sqrt(35.0), v4_.Norm()); EXPECT_EQ(Sqrt(4'126'647.0), m4_.FrobeniusNorm()); } @@ -109,7 +111,7 @@ TEST_F(UnboundedArraysTest, MultiplicationDivision) { } TEST_F(UnboundedArraysTest, Algebra) { - EXPECT_EQ(3270, v3_.Transpose() * v3_); + EXPECT_EQ(3270, TransposedView{v3_} * v3_); // NOLINT } TEST_F(UnboundedArraysTest, VectorIndexing) { @@ -230,9 +232,9 @@ TEST_F(UnboundedArraysTest, Erase) { TEST_F(UnboundedArraysTest, Transpose) { EXPECT_EQ( UnboundedMatrix({1, 8, 55, 377, - 2, 13, 89, 610, - 3, 21, 144, 987, - 5, 34, 233, 1597}), m4_.Transpose()); + 2, 13, 89, 610, + 3, 21, 144, 987, + 5, 34, 233, 1597}), m4_.Transpose()); EXPECT_EQ( UnboundedUpperTriangularMatrix({1, 2, 5, 21, 3, 8, 34,