diff --git a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md index 5f0139f..5b978e7 100644 --- a/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md +++ b/docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md @@ -85,3 +85,45 @@ - handoff: `SourceLocation`, `SourceEntityId`, `Diagnostic`, deterministic diagnostic sorting, `FailureCategory`, `Status`, and `Result` are available to Step 9 and later parser/model/solver tasks. + +## Step 9 — dense-math-adapters + +- task_id: `TASK-09` +- status: `completed` +- changed_files: `include/fesa/math/vector.hpp`, + `include/fesa/math/matrix.hpp`, `src/fesa/math/vector.cpp`, + `src/fesa/math/matrix.cpp`, `tests/unit/math/vector_test.cpp`, + `tests/unit/math/matrix_test.cpp`, `src/fesa/CMakeLists.txt`, + `tests/CMakeLists.txt`, + `docs/implementation-plans/linear-static-3d-euler-beam-implementation-report.md`, + `phases/linear-static-3d-euler-beam/index.json` +- requirement_ids: `FESA-REQ-LS3DEB-025`, `FESA-REQ-LS3DEB-034` +- test_ids: `T09-DENSE-001`, `T09-DENSE-002` + +| stage | exact command | exit_code | expected_or_observed_result | evidence_tail | +| --- | --- | ---: | --- | --- | +| RED | `cmake --build .harness/build --config Debug --target fesa_tests`; `ctest --test-dir .harness/build -C Debug -R DenseMath --output-on-failure` | 1 | Both planned tests are registered before production and the build fails for the missing dense-math API | MSVC C1083 for `fesa/math/matrix.hpp` and `fesa/math/vector.hpp`; DenseMath discovery is unavailable until the test executable builds | +| GREEN-build | `cmake --build .harness/build --config Debug --target fesa_tests` | 0 | Minimal dense adapters, tests, and runtime staging compile, link, and complete GoogleTest discovery | `matrix.cpp`, `vector.cpp`, both math tests, `fesa_solver.lib`, and `fesa_unit_tests.exe` built; no FESA warning under `/W4 /WX` | +| GREEN-test | `ctest --test-dir .harness/build -C Debug -R DenseMath --output-on-failure` | 0 | Storage/ownership/checking and row-major BLAS contracts pass | 2/2 DenseMath tests passed | +| VERIFY-configure | `cmake -S . -B .harness/build -A x64 -DFESA_GTEST_SOURCE_DIR=C:/git/googletest "-DMKL_DIR=C:/Program Files (x86)/Intel/oneAPI/mkl/2026.1/lib/cmake/mkl" "-DTBB_DIR=C:/Program Files (x86)/Intel/oneAPI/tbb/2023.1/lib/cmake/tbb" "-DHDF5_DIR=C:/Program Files/HDF_Group/HDF5/2.1.1/cmake"` | 0 | Approved MSVC x64 build tree regenerates with explicit dependencies | Configure/generate completed; oneMKL 2026.1 dynamic ILP64 backend resolved | +| VERIFY-build | `cmake --build .harness/build --config Debug` | 0 | Full Debug build passes without a new FESA warning | `fesa_solver.lib` and `fesa_unit_tests.exe` built under `/W4 /WX` | +| VERIFY-targeted | `ctest --test-dir .harness/build -C Debug -R DenseMath --output-on-failure` | 0 | Focused dense-math suite remains green | 2/2 tests passed | +| VERIFY-discovery | `ctest --test-dir .harness/build -C Debug --show-only=json-v1` | 0 | CTest discovers the accumulated tests and both exact DenseMath names | 7 tests discovered with feature/unit labels | +| VERIFY-full | `ctest --test-dir .harness/build -C Debug --output-on-failure` | 0 | Full accumulated C++ suite has zero failures | 7/7 tests passed | +| VERIFY-public-header | `$leaks = rg -n "mkl\.h\|MKL_INT\|CBLAS_" include/fesa; if ($LASTEXITCODE -eq 0) { throw "MKL API leaked into public headers:`n$leaks" }; if ($LASTEXITCODE -ne 1) { throw 'Public-header dependency scan failed' }` | 0 | Public headers remain independent of oneMKL API and types | `rg` returned the required no-match exit 1; wrapper scan completed with zero leaks | + +- contract_checks: `Vector` accepts zero size, owns contiguous doubles, deep-copies, + leaves moved-from objects empty and usable, checks every index, and rejects dot/axpy + mismatches; `Matrix` accepts zero dimensions, owns row-major contiguous doubles, + deep-copies, leaves moved-from objects empty and usable, checks every index, and + rejects GEMV/GEMM mismatches. Hand-derived nonsquare GEMV `[50,122]` and GEMM + `[[58,64],[139,154]]` pass. All copy, dot, norm, scale, axpy, GEMV, and GEMM calls + cross CBLAS only inside `.cpp` files. Windows tests stage the package-resolved + oneMKL/OpenMP runtime before post-build discovery because the exact shell does not + assume oneAPI in `PATH`. +- generated_evidence: `.harness/build/src/fesa/Debug/fesa_solver.lib`, + `.harness/build/tests/Debug/fesa_unit_tests.exe` +- reference_diff: unchanged; `git diff --exit-code -- reference/` exit 0 +- handoff: backend-free `fesa::Vector` and `fesa::Matrix` public APIs, with owning + dense storage and checked BLAS operations, are available to Step 10 and later + DOF/element/assembly tasks. diff --git a/include/fesa/math/matrix.hpp b/include/fesa/math/matrix.hpp new file mode 100644 index 0000000..1051295 --- /dev/null +++ b/include/fesa/math/matrix.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "fesa/math/vector.hpp" + +#include +#include + +namespace fesa { + +// Owns row-major contiguous dense storage independently of sparse matrices. +class Matrix { +public: + Matrix(std::size_t rows, std::size_t columns, double value = 0.0); + Matrix(const Matrix& other); + Matrix(Matrix&& other) noexcept; + Matrix& operator=(const Matrix& other); + Matrix& operator=(Matrix&& other) noexcept; + + std::size_t rows() const noexcept; + std::size_t columns() const noexcept; + double& operator()(std::size_t row, std::size_t column); + const double& operator()(std::size_t row, std::size_t column) const; + Vector multiply(const Vector& rhs) const; + Matrix multiply(const Matrix& rhs) const; + +private: + std::size_t rows_; + std::size_t columns_; + std::vector values_; +}; + +} // namespace fesa diff --git a/include/fesa/math/vector.hpp b/include/fesa/math/vector.hpp new file mode 100644 index 0000000..0842bcd --- /dev/null +++ b/include/fesa/math/vector.hpp @@ -0,0 +1,31 @@ +#pragma once + +#include +#include + +namespace fesa { + +// Owns a contiguous dense vector while keeping the MKL backend private. +class Vector { +public: + explicit Vector(std::size_t size, double value = 0.0); + Vector(const Vector& other); + Vector(Vector&& other) noexcept; + Vector& operator=(const Vector& other); + Vector& operator=(Vector&& other) noexcept; + + std::size_t size() const noexcept; + double* data() noexcept; + const double* data() const noexcept; + double& operator[](std::size_t index); + const double& operator[](std::size_t index) const; + double dot(const Vector& rhs) const; + double norm() const; + void scale(double alpha); + void axpy(double alpha, const Vector& x); + +private: + std::vector values_; +}; + +} // namespace fesa diff --git a/src/fesa/CMakeLists.txt b/src/fesa/CMakeLists.txt index a21a90a..59b9bec 100644 --- a/src/fesa/CMakeLists.txt +++ b/src/fesa/CMakeLists.txt @@ -4,6 +4,8 @@ add_library( build_info.cpp core/diagnostic.cpp core/status.cpp + math/matrix.cpp + math/vector.cpp ) target_include_directories( @@ -12,6 +14,12 @@ target_include_directories( "${PROJECT_SOURCE_DIR}/include" ) +target_link_libraries( + fesa_solver + PRIVATE + Fesa::MKL +) + # Product warnings are strict without imposing FESA policy on external targets. target_compile_options( fesa_solver diff --git a/src/fesa/math/matrix.cpp b/src/fesa/math/matrix.cpp new file mode 100644 index 0000000..c5fd90d --- /dev/null +++ b/src/fesa/math/matrix.cpp @@ -0,0 +1,150 @@ +#include "fesa/math/matrix.hpp" + +#include + +#include +#include +#include + +namespace fesa { +namespace { + +MKL_INT toMklSize(const std::size_t size) { + if (size > static_cast((std::numeric_limits::max)())) { + throw std::length_error{"Dense matrix dimension exceeds the MKL integer range."}; + } + return static_cast(size); +} + +void copyValues(const std::vector& source, std::vector& destination) { + if (source.empty()) { + return; + } + + cblas_dcopy(toMklSize(source.size()), source.data(), 1, destination.data(), 1); +} + +} // namespace + +Matrix::Matrix( + const std::size_t rows, + const std::size_t columns, + const double value) + : rows_(rows), columns_(columns), values_(rows * columns, value) {} + +Matrix::Matrix(const Matrix& other) + : rows_(other.rows_), columns_(other.columns_), values_(other.values_.size()) { + copyValues(other.values_, values_); +} + +Matrix::Matrix(Matrix&& other) noexcept + : rows_(other.rows_), + columns_(other.columns_), + values_(std::move(other.values_)) { + other.rows_ = 0; + other.columns_ = 0; + other.values_.clear(); +} + +Matrix& Matrix::operator=(const Matrix& other) { + if (this != &other) { + std::vector copied(other.values_.size()); + copyValues(other.values_, copied); + rows_ = other.rows_; + columns_ = other.columns_; + values_.swap(copied); + } + return *this; +} + +Matrix& Matrix::operator=(Matrix&& other) noexcept { + if (this != &other) { + rows_ = other.rows_; + columns_ = other.columns_; + values_ = std::move(other.values_); + other.rows_ = 0; + other.columns_ = 0; + other.values_.clear(); + } + return *this; +} + +std::size_t Matrix::rows() const noexcept { + return rows_; +} + +std::size_t Matrix::columns() const noexcept { + return columns_; +} + +double& Matrix::operator()(const std::size_t row, const std::size_t column) { + if (row >= rows_ || column >= columns_) { + throw std::out_of_range{"Matrix index is outside its dimensions."}; + } + return values_[row * columns_ + column]; +} + +const double& Matrix::operator()(const std::size_t row, const std::size_t column) const { + if (row >= rows_ || column >= columns_) { + throw std::out_of_range{"Matrix index is outside its dimensions."}; + } + return values_[row * columns_ + column]; +} + +Vector Matrix::multiply(const Vector& rhs) const { + if (columns_ != rhs.size()) { + throw std::invalid_argument{"Matrix-vector multiplication has incompatible dimensions."}; + } + + Vector result{rows_}; + if (rows_ == 0 || columns_ == 0) { + return result; + } + + // The owned layout is row-major, so the leading dimension is the column + // count for the adapter call and remains invisible to public consumers. + cblas_dgemv( + CblasRowMajor, + CblasNoTrans, + toMklSize(rows_), + toMklSize(columns_), + 1.0, + values_.data(), + toMklSize(columns_), + rhs.data(), + 1, + 0.0, + result.data(), + 1); + return result; +} + +Matrix Matrix::multiply(const Matrix& rhs) const { + if (columns_ != rhs.rows_) { + throw std::invalid_argument{"Matrix multiplication has incompatible dimensions."}; + } + + Matrix result{rows_, rhs.columns_}; + if (rows_ == 0 || columns_ == 0 || rhs.columns_ == 0) { + return result; + } + + cblas_dgemm( + CblasRowMajor, + CblasNoTrans, + CblasNoTrans, + toMklSize(rows_), + toMklSize(rhs.columns_), + toMklSize(columns_), + 1.0, + values_.data(), + toMklSize(columns_), + rhs.values_.data(), + toMklSize(rhs.columns_), + 0.0, + result.values_.data(), + toMklSize(rhs.columns_)); + return result; +} + +} // namespace fesa diff --git a/src/fesa/math/vector.cpp b/src/fesa/math/vector.cpp new file mode 100644 index 0000000..1c9c255 --- /dev/null +++ b/src/fesa/math/vector.cpp @@ -0,0 +1,119 @@ +#include "fesa/math/vector.hpp" + +#include + +#include +#include +#include + +namespace fesa { +namespace { + +MKL_INT toMklSize(const std::size_t size) { + if (size > static_cast((std::numeric_limits::max)())) { + throw std::length_error{"Dense vector size exceeds the MKL integer range."}; + } + return static_cast(size); +} + +void copyValues(const std::vector& source, std::vector& destination) { + if (source.empty()) { + return; + } + + // Keep the backend operation in this translation unit so public ownership + // remains independent of MKL headers and integer types. + cblas_dcopy(toMklSize(source.size()), source.data(), 1, destination.data(), 1); +} + +} // namespace + +Vector::Vector(const std::size_t size, const double value) + : values_(size, value) {} + +Vector::Vector(const Vector& other) + : values_(other.size()) { + copyValues(other.values_, values_); +} + +Vector::Vector(Vector&& other) noexcept + : values_(std::move(other.values_)) { + other.values_.clear(); +} + +Vector& Vector::operator=(const Vector& other) { + if (this != &other) { + std::vector copied(other.size()); + copyValues(other.values_, copied); + values_.swap(copied); + } + return *this; +} + +Vector& Vector::operator=(Vector&& other) noexcept { + if (this != &other) { + values_ = std::move(other.values_); + other.values_.clear(); + } + return *this; +} + +std::size_t Vector::size() const noexcept { + return values_.size(); +} + +double* Vector::data() noexcept { + return values_.data(); +} + +const double* Vector::data() const noexcept { + return values_.data(); +} + +double& Vector::operator[](const std::size_t index) { + return values_.at(index); +} + +const double& Vector::operator[](const std::size_t index) const { + return values_.at(index); +} + +double Vector::dot(const Vector& rhs) const { + if (size() != rhs.size()) { + throw std::invalid_argument{"Vector dot product requires equal dimensions."}; + } + if (values_.empty()) { + return 0.0; + } + + return cblas_ddot(toMklSize(size()), data(), 1, rhs.data(), 1); +} + +double Vector::norm() const { + if (values_.empty()) { + return 0.0; + } + + return cblas_dnrm2(toMklSize(size()), data(), 1); +} + +void Vector::scale(const double alpha) { + if (values_.empty()) { + return; + } + + cblas_dscal(toMklSize(size()), alpha, data(), 1); +} + +void Vector::axpy(const double alpha, const Vector& x) { + if (size() != x.size()) { + throw std::invalid_argument{"Vector axpy requires equal dimensions."}; + } + if (values_.empty()) { + return; + } + + cblas_daxpy(toMklSize(size()), alpha, x.data(), 1, data(), 1); +} + +} // namespace fesa diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9337a44..058f733 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -6,6 +6,8 @@ add_executable( unit/core/diagnostic_test.cpp unit/core/source_identity_test.cpp unit/core/status_test.cpp + unit/math/matrix_test.cpp + unit/math/vector_test.cpp ) target_link_libraries( @@ -15,6 +17,42 @@ target_link_libraries( GTest::gtest_main ) +if(WIN32) + if(NOT TARGET MKL::mkl_intel_thread OR NOT TARGET MKL::mkl_core OR + NOT OMP_DLL_DIR OR NOT OMP_DLLNAME) + message(FATAL_ERROR "The oneMKL runtime files required by tests were not resolved") + endif() + + file(GLOB _fesa_mkl_dispatch_dlls "${MKL_ROOT}/bin/mkl_def.*.dll") + if(NOT _fesa_mkl_dispatch_dlls) + message(FATAL_ERROR "The oneMKL default dispatch runtime was not resolved") + endif() + list(SORT _fesa_mkl_dispatch_dlls COMPARE NATURAL ORDER DESCENDING) + list(GET _fesa_mkl_dispatch_dlls 0 _fesa_mkl_dispatch_dll) + + # Stage dynamic backend dependencies before post-build GoogleTest + # discovery; the exact validation commands do not assume oneAPI in PATH. + add_custom_command( + TARGET fesa_unit_tests + POST_BUILD + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "$" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${OMP_DLL_DIR}/${OMP_DLLNAME}" + "$" + COMMAND "${CMAKE_COMMAND}" -E copy_if_different + "${_fesa_mkl_dispatch_dll}" + "$" + ) + + unset(_fesa_mkl_dispatch_dll) + unset(_fesa_mkl_dispatch_dlls) +endif() + gtest_discover_tests( fesa_unit_tests PROPERTIES diff --git a/tests/unit/math/matrix_test.cpp b/tests/unit/math/matrix_test.cpp new file mode 100644 index 0000000..11866dc --- /dev/null +++ b/tests/unit/math/matrix_test.cpp @@ -0,0 +1,104 @@ +#include "fesa/math/matrix.hpp" + +#include + +#include +#include + +namespace fesa { +namespace { + +TEST(DenseMath, RowMajorMatrixMatchesKnownGemvGemm) { + Matrix zeroRows{0, 3}; + Vector threeValues{3, 2.0}; + const Vector zeroRowProduct = zeroRows.multiply(threeValues); + EXPECT_EQ(zeroRowProduct.size(), 0U); + + Matrix zeroColumns{2, 0}; + const Vector zeroColumnProduct = zeroColumns.multiply(Vector{0}); + ASSERT_EQ(zeroColumnProduct.size(), 2U); + EXPECT_DOUBLE_EQ(zeroColumnProduct[0], 0.0); + EXPECT_DOUBLE_EQ(zeroColumnProduct[1], 0.0); + + Matrix zeroInnerRight{0, 3}; + const Matrix zeroInnerProduct = zeroColumns.multiply(zeroInnerRight); + EXPECT_EQ(zeroInnerProduct.rows(), 2U); + EXPECT_EQ(zeroInnerProduct.columns(), 3U); + for (std::size_t row = 0; row < zeroInnerProduct.rows(); ++row) { + for (std::size_t column = 0; column < zeroInnerProduct.columns(); ++column) { + EXPECT_DOUBLE_EQ(zeroInnerProduct(row, column), 0.0); + } + } + + Matrix left{2, 3}; + left(0, 0) = 1.0; + left(0, 1) = 2.0; + left(0, 2) = 3.0; + left(1, 0) = 4.0; + left(1, 1) = 5.0; + left(1, 2) = 6.0; + + EXPECT_EQ(&left(0, 0) + 1, &left(0, 1)); + EXPECT_EQ(&left(0, 0) + 2, &left(0, 2)); + EXPECT_EQ(&left(0, 0) + 3, &left(1, 0)); + const Matrix& constLeft = left; + EXPECT_DOUBLE_EQ(constLeft(1, 2), 6.0); + + Matrix copied{left}; + copied(0, 0) = 42.0; + EXPECT_DOUBLE_EQ(left(0, 0), 1.0); + + Matrix copyAssigned{0, 0}; + copyAssigned = left; + copyAssigned(1, 2) = -7.0; + EXPECT_DOUBLE_EQ(left(1, 2), 6.0); + + Matrix moved{std::move(copied)}; + EXPECT_EQ(copied.rows(), 0U); + EXPECT_EQ(copied.columns(), 0U); + EXPECT_EQ(moved.rows(), 2U); + EXPECT_EQ(moved.columns(), 3U); + EXPECT_DOUBLE_EQ(moved(0, 0), 42.0); + EXPECT_NO_THROW(static_cast(copied.multiply(Vector{0}))); + + Matrix moveAssigned{1, 1, -1.0}; + moveAssigned = std::move(copyAssigned); + EXPECT_EQ(copyAssigned.rows(), 0U); + EXPECT_EQ(copyAssigned.columns(), 0U); + EXPECT_EQ(moveAssigned.rows(), 2U); + EXPECT_EQ(moveAssigned.columns(), 3U); + EXPECT_DOUBLE_EQ(moveAssigned(1, 2), -7.0); + + Vector vector{3}; + vector[0] = 7.0; + vector[1] = 8.0; + vector[2] = 9.0; + const Vector matrixVectorProduct = left.multiply(vector); + ASSERT_EQ(matrixVectorProduct.size(), 2U); + EXPECT_DOUBLE_EQ(matrixVectorProduct[0], 50.0); + EXPECT_DOUBLE_EQ(matrixVectorProduct[1], 122.0); + + Matrix right{3, 2}; + right(0, 0) = 7.0; + right(0, 1) = 8.0; + right(1, 0) = 9.0; + right(1, 1) = 10.0; + right(2, 0) = 11.0; + right(2, 1) = 12.0; + const Matrix matrixProduct = left.multiply(right); + ASSERT_EQ(matrixProduct.rows(), 2U); + ASSERT_EQ(matrixProduct.columns(), 2U); + EXPECT_DOUBLE_EQ(matrixProduct(0, 0), 58.0); + EXPECT_DOUBLE_EQ(matrixProduct(0, 1), 64.0); + EXPECT_DOUBLE_EQ(matrixProduct(1, 0), 139.0); + EXPECT_DOUBLE_EQ(matrixProduct(1, 1), 154.0); + + EXPECT_THROW(static_cast(left(2, 0)), std::out_of_range); + EXPECT_THROW(static_cast(left(0, 3)), std::out_of_range); + EXPECT_THROW(static_cast(constLeft(2, 0)), std::out_of_range); + EXPECT_THROW(static_cast(left.multiply(Vector{2})), std::invalid_argument); + EXPECT_THROW(static_cast(left.multiply(Matrix{4, 1})), std::invalid_argument); +} + +} // namespace +} // namespace fesa diff --git a/tests/unit/math/vector_test.cpp b/tests/unit/math/vector_test.cpp new file mode 100644 index 0000000..00a0ba4 --- /dev/null +++ b/tests/unit/math/vector_test.cpp @@ -0,0 +1,86 @@ +#include "fesa/math/vector.hpp" + +#include + +#include +#include +#include + +namespace fesa { +namespace { + +TEST(DenseMath, VectorOwnsAndChecksContiguousStorage) { + Vector empty{0}; + EXPECT_EQ(empty.size(), 0U); + EXPECT_DOUBLE_EQ(empty.norm(), 0.0); + EXPECT_NO_THROW(empty.scale(3.0)); + EXPECT_NO_THROW(empty.axpy(-2.0, Vector{0})); + EXPECT_THROW(static_cast(empty[0]), std::out_of_range); + + Vector original{3}; + original[0] = 1.0; + original[1] = -2.0; + original[2] = 3.0; + + EXPECT_EQ(original.data() + 1, &original[1]); + EXPECT_EQ(original.data() + 2, &original[2]); + const Vector& constOriginal = original; + EXPECT_EQ(constOriginal.data() + 2, &constOriginal[2]); + + Vector copied{original}; + EXPECT_NE(copied.data(), original.data()); + copied[0] = 99.0; + EXPECT_DOUBLE_EQ(original[0], 1.0); + + Vector copyAssigned{0}; + copyAssigned = original; + EXPECT_NE(copyAssigned.data(), original.data()); + copyAssigned[1] = 17.0; + EXPECT_DOUBLE_EQ(original[1], -2.0); + + Vector moved{std::move(copied)}; + EXPECT_EQ(copied.size(), 0U); + EXPECT_EQ(moved.size(), 3U); + EXPECT_DOUBLE_EQ(moved[0], 99.0); + EXPECT_NO_THROW(copied.scale(4.0)); + + Vector moveAssigned{1, -1.0}; + moveAssigned = std::move(copyAssigned); + EXPECT_EQ(copyAssigned.size(), 0U); + EXPECT_EQ(moveAssigned.size(), 3U); + EXPECT_DOUBLE_EQ(moveAssigned[1], 17.0); + + Vector rhs{3}; + rhs[0] = 4.0; + rhs[1] = 5.0; + rhs[2] = -6.0; + EXPECT_DOUBLE_EQ(original.dot(rhs), -24.0); + EXPECT_NEAR(original.norm(), std::sqrt(14.0), 1.0e-15); + + Vector scaled{original}; + scaled.scale(-0.5); + EXPECT_DOUBLE_EQ(scaled[0], -0.5); + EXPECT_DOUBLE_EQ(scaled[1], 1.0); + EXPECT_DOUBLE_EQ(scaled[2], -1.5); + + Vector accumulated{3}; + accumulated[0] = 1.0; + accumulated[1] = 2.0; + accumulated[2] = 3.0; + Vector increment{3}; + increment[0] = 4.0; + increment[1] = -1.0; + increment[2] = 0.5; + accumulated.axpy(2.0, increment); + EXPECT_DOUBLE_EQ(accumulated[0], 9.0); + EXPECT_DOUBLE_EQ(accumulated[1], 0.0); + EXPECT_DOUBLE_EQ(accumulated[2], 4.0); + + EXPECT_THROW(static_cast(original[3]), std::out_of_range); + EXPECT_THROW(static_cast(constOriginal[3]), std::out_of_range); + EXPECT_THROW(static_cast(original.dot(Vector{2})), std::invalid_argument); + EXPECT_THROW(original.axpy(1.0, Vector{2}), std::invalid_argument); +} + +} // namespace +} // namespace fesa