diff --git a/benchmarks/benchmark_matrix.cpp b/benchmarks/benchmark_matrix.cpp index 662a84a..b4e9c78 100644 --- a/benchmarks/benchmark_matrix.cpp +++ b/benchmarks/benchmark_matrix.cpp @@ -12,7 +12,7 @@ using namespace cramer; // Helper function to create a random matrix template Matrix createRandomMatrix(size_t rows, size_t cols) { - static std::random_device rd; + static std::random_device rand; static std::mt19937 gen(rd()); static std::uniform_real_distribution<> dis(-1.0, 1.0); @@ -28,9 +28,9 @@ Matrix createRandomMatrix(size_t rows, size_t cols) { // Benchmark Matrix Constructor static void BM_MatrixConstructor(benchmark::State& state) { const size_t size = state.range(0); - for (auto _ : state) { - Matrix m(size, size); - benchmark::DoNotOptimize(m); + for (auto iter : state) { + Matrix matrix(size, size); + benchmark::DoNotOptimize(matrix); } } BENCHMARK(BM_MatrixConstructor)->Range(8, 1024); @@ -38,10 +38,10 @@ BENCHMARK(BM_MatrixConstructor)->Range(8, 1024); // Benchmark Matrix Addition static void BM_MatrixAddition(benchmark::State& state) { const size_t size = state.range(0); - auto m1 = createRandomMatrix(size, size); - auto m2 = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m1 + m2; + auto matrix1 = createRandomMatrix(size, size); + auto matrix2 = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix1 + matrix2; benchmark::DoNotOptimize(result); } } @@ -50,10 +50,10 @@ BENCHMARK(BM_MatrixAddition)->Range(8, 1024); // Benchmark Matrix Subtraction static void BM_MatrixSubtraction(benchmark::State& state) { const size_t size = state.range(0); - auto m1 = createRandomMatrix(size, size); - auto m2 = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m1 - m2; + auto matrix1 = createRandomMatrix(size, size); + auto matrix2 = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix1 - matrix2; benchmark::DoNotOptimize(result); } } @@ -62,10 +62,10 @@ BENCHMARK(BM_MatrixSubtraction)->Range(8, 1024); // Benchmark Matrix Multiplication static void BM_MatrixMultiplication(benchmark::State& state) { const size_t size = state.range(0); - auto m1 = createRandomMatrix(size, size); - auto m2 = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m1 * m2; + auto matrix1 = createRandomMatrix(size, size); + auto matrix2 = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix1 * matrix2; benchmark::DoNotOptimize(result); } } @@ -74,10 +74,10 @@ BENCHMARK(BM_MatrixMultiplication)->Range(8, 256); // Benchmark Matrix Scalar Multiplication static void BM_MatrixScalarMultiplication(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); + auto matrix = createRandomMatrix(size, size); double scalar = 2.0; - for (auto _ : state) { - auto result = m * scalar; + for (auto iter : state) { + auto result = matrix * 2.0; benchmark::DoNotOptimize(result); } } @@ -86,9 +86,9 @@ BENCHMARK(BM_MatrixScalarMultiplication)->Range(8, 1024); // Benchmark Matrix Transpose static void BM_MatrixTranspose(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.transpose(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.transpose(); benchmark::DoNotOptimize(result); } } @@ -97,9 +97,9 @@ BENCHMARK(BM_MatrixTranspose)->Range(8, 1024); // Benchmark Matrix Determinant static void BM_MatrixDeterminant(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.det(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.det(); benchmark::DoNotOptimize(result); } } @@ -108,9 +108,9 @@ BENCHMARK(BM_MatrixDeterminant)->Range(2, 128); // Benchmark Matrix Inverse static void BM_MatrixInverse(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.inverse(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.inverse(); benchmark::DoNotOptimize(result); } } @@ -119,9 +119,9 @@ BENCHMARK(BM_MatrixInverse)->Range(2, 128); // Benchmark Matrix LU Decomposition static void BM_MatrixLU(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto [L, U] = m.lu(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto [L, U] = matrix.lu(); benchmark::DoNotOptimize(L); benchmark::DoNotOptimize(U); } @@ -131,9 +131,9 @@ BENCHMARK(BM_MatrixLU)->Range(2, 128); // Benchmark Matrix QR Decomposition static void BM_MatrixQR(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto [Q, R] = m.qr(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto [Q, R] = matrix.qr(); benchmark::DoNotOptimize(Q); benchmark::DoNotOptimize(R); } @@ -143,9 +143,9 @@ BENCHMARK(BM_MatrixQR)->Range(2, 128); // Benchmark Matrix SVD static void BM_MatrixSVD(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto [U, S, V] = m.svd(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto [U, S, V] = matrix.svd(); benchmark::DoNotOptimize(U); benchmark::DoNotOptimize(S); benchmark::DoNotOptimize(V); @@ -156,9 +156,9 @@ BENCHMARK(BM_MatrixSVD)->Range(2, 64); // Benchmark Matrix Eigenvalues static void BM_MatrixEigenvalues(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.eigenvalues(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.eigenvalues(); benchmark::DoNotOptimize(result); } } @@ -167,9 +167,9 @@ BENCHMARK(BM_MatrixEigenvalues)->Range(2, 64); // Benchmark Matrix Eigenvectors static void BM_MatrixEigenvectors(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.eigenvectors(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.eigenvectors(); benchmark::DoNotOptimize(result); } } @@ -178,9 +178,9 @@ BENCHMARK(BM_MatrixEigenvectors)->Range(2, 64); // Benchmark Matrix Rank static void BM_MatrixRank(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.rank(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.rank(); benchmark::DoNotOptimize(result); } } @@ -189,9 +189,9 @@ BENCHMARK(BM_MatrixRank)->Range(2, 128); // Benchmark Matrix Trace static void BM_MatrixTrace(benchmark::State& state) { const size_t size = state.range(0); - auto m = createRandomMatrix(size, size); - for (auto _ : state) { - auto result = m.trace(); + auto matrix = createRandomMatrix(size, size); + for (auto iter : state) { + auto result = matrix.trace(); benchmark::DoNotOptimize(result); } } diff --git a/include/matrix.hpp b/include/matrix.hpp index 19b2356..f10ac93 100644 --- a/include/matrix.hpp +++ b/include/matrix.hpp @@ -264,25 +264,6 @@ class Matrix { */ Matrix conjugate() const; - /** - * @brief Calculates the matrix exponential. - * @return The exponential of the matrix. - */ - Matrix exp() const; - - /** - * @brief Calculates the matrix power. - * @param n The power to raise the matrix to. - * @return The matrix raised to the power n. - */ - Matrix pow(int n) const; - - /** - * @brief Calculates the square root of the matrix. - * @return The square root of the matrix. - */ - Matrix sqrt() const; - /** * @brief Performs LU decomposition of the matrix. * @return A pair of matrices (L, U) where L is lower triangular and U is diff --git a/src/matrix.cpp b/src/matrix.cpp index ef04ee4..4f771a5 100644 --- a/src/matrix.cpp +++ b/src/matrix.cpp @@ -621,63 +621,6 @@ Matrix Matrix::exp() const { return F; } -template -Matrix Matrix::pow(int n) const { - if (!is_square()) { - throw std::invalid_argument( - "Matrix must be square to calculate its integer power."); - } - - if (n == 0) { - return identity(rows); - } - - Matrix result = identity(rows); - Matrix base = *this; - bool negative = (n < 0); - n = std::abs(n); - - while (n > 0) { - if (n % 2 == 1) { - result = result * base; // Use matrix multiplication instead of *= - } - base = base * base; // Use matrix multiplication instead of *= - n /= 2; - } - - return negative ? result.inverse() : result; -} - -template -Matrix Matrix::sqrt() const { - if (!is_square()) { - throw std::invalid_argument( - "Matrix must be square to calculate its square root."); - } - - // This is a placeholder implementation using Denman-Beavers iteration - // For a more robust implementation, consider using Schur decomposition or - // other methods - Matrix X = *this; - Matrix Y = identity(rows); - const int max_iterations = 20; - const T tolerance = std::numeric_limits::epsilon(); - - for (int i = 0; i < max_iterations; ++i) { - Matrix X_next = (X + Y.inverse()) * static_cast(0.5); - Matrix Y_next = (Y + X.inverse()) * static_cast(0.5); - - if (std::abs((X_next - X).max_norm()) < std::abs(tolerance)) { - return X_next; - } - - X = X_next; - Y = Y_next; - } - - throw std::runtime_error("Square root iteration did not converge."); -} - // Matrix decompositions template @@ -717,48 +660,6 @@ std::pair, Matrix> Matrix::lu() const { return std::make_pair(L, U); } -template -Matrix Matrix::exp() const { - if (!is_square()) { - throw std::invalid_argument("Matrix must be square to calculate the exponential."); - } - - const int q = 6; - const int p = 1; - T norm = max_norm(); - int s = std::max(0, static_cast(std::ceil(std::log2(norm)))); - - Matrix A = *this * (1.0 / std::pow(2.0, s)); - Matrix X = identity(rows) + A; - Matrix cX = identity(rows) - A; - Matrix A2 = A * A; - - for (int k = 2; k <= q; ++k) { - T c = 1.0; - for (int j = 1; j <= std::min(k, p); ++j) { - c *= static_cast(k - j + 1) / static_cast(j * (2 * k - j + 1)); - } - Matrix Y = c * A2; - X += Y; - if (k % 2 == 0) { - cX += Y; - } else { - cX -= Y; - } - if (k < q) { - A2 *= A; - } - } - - Matrix E = X * cX.inverse(); - - for (int k = 0; k < s; ++k) { - E = E * E; - } - - return E; -} - template std::tuple, Matrix, Matrix> Matrix::svd() const { const size_t m = rows; diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp index 1f1d43d..f9392ac 100644 --- a/tests/test_matrix.cpp +++ b/tests/test_matrix.cpp @@ -334,33 +334,6 @@ TEST_F(MatrixTest, Conjugate) { EXPECT_EQ(result(1, 1), std::complex(4, -4)); } -TEST_F(MatrixTest, Exp) { - Matrix mat(2, 2, {0, 1, -1, 0}); - Matrix result = mat.exp(); - EXPECT_NEAR(result(0, 0), std::cos(1), 1e-9); - EXPECT_NEAR(result(0, 1), std::sin(1), 1e-9); - EXPECT_NEAR(result(1, 0), -std::sin(1), 1e-9); - EXPECT_NEAR(result(1, 1), std::cos(1), 1e-9); -} - -TEST_F(MatrixTest, Pow) { - Matrix mat(2, 2, {1, 1, 1, 0}); - Matrix result = mat.pow(3); - EXPECT_EQ(result(0, 0), 3); - EXPECT_EQ(result(0, 1), 2); - EXPECT_EQ(result(1, 0), 2); - EXPECT_EQ(result(1, 1), 1); -} - -TEST_F(MatrixTest, Sqrt) { - Matrix mat(2, 2, {4, 0, 0, 9}); - Matrix result = mat.sqrt(); - EXPECT_NEAR(result(0, 0), 2, 1e-9); - EXPECT_NEAR(result(0, 1), 0, 1e-9); - EXPECT_NEAR(result(1, 0), 0, 1e-9); - EXPECT_NEAR(result(1, 1), 3, 1e-9); -} - TEST_F(MatrixTest, LU) { Matrix mat(3, 3, {1, 2, 3, 4, 5, 6, 7, 8, 9}); auto [lower, upper] = mat.lu();