From 392c074dbf172cf4f0e3989c031d6d830180d39d Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 3 Dec 2024 13:22:27 +0000 Subject: [PATCH 1/2] Refactor sparse_matrix_multiply using SparseMatrixView. Adapt interpolation to hold SparseMatrixStorage and create views Co-authored-by: Liam Adams --- src/atlas/CMakeLists.txt | 6 + src/atlas/interpolation/Cache.cc | 2 +- src/atlas/interpolation/Cache.h | 4 +- src/atlas/interpolation/method/Method.cc | 25 +- src/atlas/interpolation/method/Method.h | 18 +- .../interpolation/method/binning/Binning.cc | 26 +- .../method/cubedsphere/CubedSphereBilinear.cc | 3 +- .../method/knn/GridBoxMaximum.cc | 5 +- .../interpolation/method/knn/GridBoxMethod.cc | 3 +- .../method/knn/KNearestNeighbours.cc | 3 +- .../method/knn/NearestNeighbour.cc | 3 +- .../method/sphericalvector/SphericalVector.cc | 14 +- .../structured/StructuredInterpolation2D.tcc | 3 +- ...nservativeSphericalPolygonInterpolation.cc | 14 +- ...onservativeSphericalPolygonInterpolation.h | 4 +- .../method/unstructured/FiniteElement.cc | 7 +- .../UnstructuredBilinearLonLat.cc | 7 +- src/atlas/linalg/sparse.h | 4 + .../linalg/sparse/MakeEckitSparseMatrix.h | 85 +++++ .../sparse/MakeSparseMatrixStorageEckit.h | 135 +++++++ .../sparse/MakeSparseMatrixStorageEigen.h | 142 +++++++ .../linalg/sparse/SparseMatrixMultiply.h | 9 +- .../linalg/sparse/SparseMatrixMultiply.tcc | 52 ++- .../SparseMatrixMultiply_EckitLinalg.cc | 52 ++- .../sparse/SparseMatrixMultiply_EckitLinalg.h | 41 +- .../sparse/SparseMatrixMultiply_OpenMP.cc | 126 ++++--- .../sparse/SparseMatrixMultiply_OpenMP.h | 48 +-- .../linalg/sparse/SparseMatrixStorage.cc | 131 +++++++ src/atlas/linalg/sparse/SparseMatrixStorage.h | 319 ++++++++++++++++ src/atlas/linalg/sparse/SparseMatrixView.h | 64 ++++ ...est_interpolation_finite_element_cached.cc | 19 +- ...test_interpolation_k_nearest_neighbours.cc | 4 +- .../test_interpolation_non_linear.cc | 1 - src/tests/linalg/CMakeLists.txt | 17 + src/tests/linalg/test_linalg_sparse.cc | 4 +- src/tests/linalg/test_linalg_sparse_matrix.cc | 239 ++++++++++++ .../linalg/test_linalg_sparse_matrix_gpu.cc | 353 ++++++++++++++++++ 37 files changed, 1769 insertions(+), 223 deletions(-) create mode 100644 src/atlas/linalg/sparse/MakeEckitSparseMatrix.h create mode 100644 src/atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h create mode 100644 src/atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h create mode 100644 src/atlas/linalg/sparse/SparseMatrixStorage.cc create mode 100644 src/atlas/linalg/sparse/SparseMatrixStorage.h create mode 100644 src/atlas/linalg/sparse/SparseMatrixView.h create mode 100644 src/tests/linalg/test_linalg_sparse_matrix.cc create mode 100644 src/tests/linalg/test_linalg_sparse_matrix_gpu.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 60b93b9a1..1ea3b58dd 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -705,12 +705,18 @@ linalg/View.h linalg/sparse.h linalg/sparse/Backend.h linalg/sparse/Backend.cc +linalg/sparse/MakeEckitSparseMatrix.h +linalg/sparse/MakeSparseMatrixStorageEckit.h +linalg/sparse/MakeSparseMatrixStorageEigen.h linalg/sparse/SparseMatrixMultiply.h linalg/sparse/SparseMatrixMultiply.tcc linalg/sparse/SparseMatrixMultiply_EckitLinalg.h linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc linalg/sparse/SparseMatrixMultiply_OpenMP.h linalg/sparse/SparseMatrixMultiply_OpenMP.cc +linalg/sparse/SparseMatrixStorage.cc +linalg/sparse/SparseMatrixStorage.h +linalg/sparse/SparseMatrixView.h linalg/dense.h linalg/dense/Backend.h linalg/dense/Backend.cc diff --git a/src/atlas/interpolation/Cache.cc b/src/atlas/interpolation/Cache.cc index bae12b89c..d8051221c 100644 --- a/src/atlas/interpolation/Cache.cc +++ b/src/atlas/interpolation/Cache.cc @@ -54,7 +54,7 @@ MatrixCacheEntry::~MatrixCacheEntry() = default; class MatrixCacheEntryOwned : public MatrixCacheEntry { public: MatrixCacheEntryOwned(Matrix&& matrix): MatrixCacheEntry(&matrix_) { - const_cast(matrix_).swap(reinterpret_cast(matrix)); + const_cast(matrix_) = std::move(matrix); } private: diff --git a/src/atlas/interpolation/Cache.h b/src/atlas/interpolation/Cache.h index 2e70f827c..96e0447f8 100644 --- a/src/atlas/interpolation/Cache.h +++ b/src/atlas/interpolation/Cache.h @@ -17,8 +17,8 @@ #include "eckit/filesystem/PathName.h" #include "eckit/io/Buffer.h" -#include "eckit/linalg/SparseMatrix.h" +#include "atlas/linalg/sparse.h" #include "atlas/runtime/Exception.h" #include "atlas/util/KDTree.h" @@ -82,7 +82,7 @@ class Cache { class MatrixCacheEntry : public InterpolationCacheEntry { public: - using Matrix = eckit::linalg::SparseMatrix; + using Matrix = atlas::linalg::SparseMatrixStorage; ~MatrixCacheEntry() override; MatrixCacheEntry(const Matrix* matrix, const std::string& uid = ""): matrix_{matrix}, uid_(uid) { ATLAS_ASSERT(matrix_ != nullptr); diff --git a/src/atlas/interpolation/method/Method.cc b/src/atlas/interpolation/method/Method.cc index 1195a7d53..a11fc1c65 100644 --- a/src/atlas/interpolation/method/Method.cc +++ b/src/atlas/interpolation/method/Method.cc @@ -105,12 +105,12 @@ void Method::interpolate_field_rank1(const Field& src, Field& tgt, const Matrix& auto tgt_v = array::make_view(tgt); if (nonLinear_(src)) { - Matrix W_nl(W); // copy (a big penalty -- copy-on-write would definitely be better) + eckit::linalg::SparseMatrix W_nl = atlas::linalg::make_eckit_sparse_matrix(W); nonLinear_->execute(W_nl, src); sparse_matrix_multiply(W_nl, src_v, tgt_v, backend); } else { - sparse_matrix_multiply(W, src_v, tgt_v, backend); + sparse_matrix_multiply(make_host_view(W), src_v, tgt_v, backend); } } @@ -150,7 +150,7 @@ void Method::interpolate_field_rank2(const Field& src, Field& tgt, const Matrix& } } else { - sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); + sparse_matrix_multiply(make_host_view(W), src_v, tgt_v, sparse::backend::openmp()); } } @@ -163,7 +163,7 @@ void Method::interpolate_field_rank3(const Field& src, Field& tgt, const Matrix& if (not W.empty() && nonLinear_(src)) { ATLAS_ASSERT(false, "nonLinear interpolation not supported for rank-3 fields."); } - sparse_matrix_multiply(W, src_v, tgt_v, sparse::backend::openmp()); + sparse_matrix_multiply(make_host_view(W), src_v, tgt_v, sparse::backend::openmp()); } template @@ -173,7 +173,7 @@ void Method::adjoint_interpolate_field_rank1(Field& src, const Field& tgt, const auto src_v = array::make_view(src); auto tgt_v = array::make_view(tgt); - sparse_matrix_multiply_add(W, tgt_v, src_v, backend); + sparse_matrix_multiply_add(make_host_view(W), tgt_v, src_v, backend); } template @@ -182,7 +182,7 @@ void Method::adjoint_interpolate_field_rank2(Field& src, const Field& tgt, const auto src_v = array::make_view(src); auto tgt_v = array::make_view(tgt); - sparse_matrix_multiply_add(W, tgt_v, src_v, sparse::backend::openmp()); + sparse_matrix_multiply_add(make_host_view(W), tgt_v, src_v, sparse::backend::openmp()); } template @@ -191,7 +191,7 @@ void Method::adjoint_interpolate_field_rank3(Field& src, const Field& tgt, const auto src_v = array::make_view(src); auto tgt_v = array::make_view(tgt); - sparse_matrix_multiply_add(W, tgt_v, src_v, sparse::backend::openmp()); + sparse_matrix_multiply_add(make_host_view(W), tgt_v, src_v, sparse::backend::openmp()); } void Method::check_compatibility(const Field& src, const Field& tgt, const Matrix& W) const { @@ -266,12 +266,11 @@ void Method::setup(const FunctionSpace& source, const FunctionSpace& target) { ATLAS_TRACE("atlas::interpolation::method::Method::setup(FunctionSpace, FunctionSpace)"); this->do_setup(source, target); - if (adjoint_ && target.size() > 0) { - Matrix tmp(*matrix_); - - // if interpolation is matrix free then matrix->nonZeros() will be zero. - if (tmp.nonZeros() > 0) { - matrix_transpose_ = tmp.transpose(); + if (adjoint_ && target.size() > 0 && matrixAllocated()) { + if (not matrix_->empty()) { + eckit::linalg::SparseMatrix matrix_copy = make_eckit_sparse_matrix(*matrix_); // Makes a copy! + matrix_copy.transpose(); // transpose the copy in place + matrix_transpose_ = linalg::make_sparse_matrix_storage(std::move(matrix_copy)); // Move the copy into storage } } } diff --git a/src/atlas/interpolation/method/Method.h b/src/atlas/interpolation/method/Method.h index 7e52d5bf3..994ad9f7c 100644 --- a/src/atlas/interpolation/method/Method.h +++ b/src/atlas/interpolation/method/Method.h @@ -16,10 +16,12 @@ #include "atlas/interpolation/Cache.h" #include "atlas/interpolation/NonLinear.h" +#include "atlas/linalg/sparse/SparseMatrixStorage.h" #include "atlas/util/Metadata.h" #include "atlas/util/Object.h" #include "eckit/config/Configuration.h" #include "eckit/linalg/SparseMatrix.h" +#include "atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h" namespace atlas { class Field; @@ -87,7 +89,7 @@ class Method : public util::Object { using Triplet = eckit::linalg::Triplet; using Triplets = std::vector; - using Matrix = eckit::linalg::SparseMatrix; + using Matrix = atlas::linalg::SparseMatrixStorage; static void normalise(Triplets& triplets); @@ -102,13 +104,15 @@ class Method : public util::Object { friend class interpolation::MatrixCache; protected: - void setMatrix(Matrix& m, const std::string& uid = "") { + + void setMatrix(Matrix&& m, const std::string& uid = "") { if (not matrix_shared_) { matrix_shared_ = std::make_shared(); } - matrix_shared_->swap(m); + *matrix_shared_ = std::move(m); matrix_cache_ = interpolation::MatrixCache(matrix_shared_, uid); matrix_ = &matrix_cache_.matrix(); + } void setMatrix(interpolation::MatrixCache matrix_cache) { @@ -118,6 +122,14 @@ class Method : public util::Object { matrix_shared_.reset(); } + void setMatrix(eckit::linalg::SparseMatrix&& m, const std::string& uid = "") { + setMatrix( linalg::make_sparse_matrix_storage(std::move(m)), uid ); + } + + void setMatrix(std::size_t rows, std::size_t cols, const Triplets& triplets, const std::string& uid = "") { + setMatrix( eckit::linalg::SparseMatrix{rows, cols, triplets}, uid); + } + bool matrixAllocated() const { return matrix_shared_.use_count(); } const Matrix& matrix() const { return *matrix_; } diff --git a/src/atlas/interpolation/method/binning/Binning.cc b/src/atlas/interpolation/method/binning/Binning.cc index f5a6b0743..b3bcfef86 100644 --- a/src/atlas/interpolation/method/binning/Binning.cc +++ b/src/atlas/interpolation/method/binning/Binning.cc @@ -5,6 +5,10 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ +#include "eckit/config/LocalConfiguration.h" +#include "eckit/linalg/Triplet.h" +#include "eckit/linalg/SparseMatrix.h" +#include "eckit/mpi/Comm.h" #include "atlas/functionspace/NodeColumns.h" #include "atlas/functionspace/StructuredColumns.h" @@ -12,14 +16,12 @@ #include "atlas/interpolation/Interpolation.h" #include "atlas/interpolation/method/binning/Binning.h" #include "atlas/interpolation/method/MethodFactory.h" +#include "atlas/linalg/sparse/SparseMatrixStorage.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" #include "atlas/mesh.h" #include "atlas/mesh/actions/GetCubedSphereNodalArea.h" #include "atlas/runtime/Trace.h" -#include "eckit/config/LocalConfiguration.h" -#include "eckit/linalg/SparseMatrix.h" -#include "eckit/linalg/Triplet.h" -#include "eckit/mpi/Comm.h" namespace atlas { @@ -56,7 +58,8 @@ void Binning::do_setup(const FunctionSpace& source, const FunctionSpace& target) { ATLAS_TRACE("atlas::interpolation::method::Binning::do_setup()"); - using Index = eckit::linalg::Index; + using Scalar = eckit::linalg::Scalar; + using Index = eckit::linalg::Index; using Triplet = eckit::linalg::Triplet; using SMatrix = eckit::linalg::SparseMatrix; @@ -80,14 +83,16 @@ void Binning::do_setup(const FunctionSpace& source, auto smx_interp = smx_interp_cache.matrix(); - auto smx_interp_tr = smx_interp.transpose(); + auto eckit_smx_interp = make_non_owning_eckit_sparse_matrix(smx_interp); + SMatrix smx_interp_tr(eckit_smx_interp); // copy + smx_interp_tr.transpose(); // transpose the copy in-place const auto rows_tamx = smx_interp_tr.rows(); const auto cols_tamx = smx_interp_tr.cols(); - const double* ptr_tamx_data = smx_interp_tr.data(); - const Index* ptr_tamx_idxs_col = smx_interp_tr.inner(); - const Index* ptr_tamx_o = smx_interp_tr.outer(); + const Scalar* ptr_tamx_data = smx_interp_tr.data(); + const Index* ptr_tamx_idxs_col = smx_interp_tr.inner(); + const Index* ptr_tamx_o = smx_interp_tr.outer(); // diagonal of 'area weights matrix', W auto ds_aweights = getAreaWeights(source_); @@ -123,8 +128,7 @@ void Binning::do_setup(const FunctionSpace& source, } // 'binning matrix' (sparse matrix), B = N A^T W - SMatrix smx_binning{rows_tamx, cols_tamx, smx_binning_els}; - setMatrix(smx_binning); + setMatrix(rows_tamx, cols_tamx, smx_binning_els); } diff --git a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc index 366df8cfe..5966861f0 100644 --- a/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc +++ b/src/atlas/interpolation/method/cubedsphere/CubedSphereBilinear.cc @@ -99,8 +99,7 @@ void CubedSphereBilinear::do_setup(const FunctionSpace& source, const FunctionSp } // fill sparse matrix and return. - Matrix A(target_.size(), source_.size(), weights); - setMatrix(A); + setMatrix(target_.size(), source_.size(), weights); // Add tile index metadata to target. target_->metadata().set("tile index", tileIndex); diff --git a/src/atlas/interpolation/method/knn/GridBoxMaximum.cc b/src/atlas/interpolation/method/knn/GridBoxMaximum.cc index 8ca86672a..2db2e50d9 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMaximum.cc +++ b/src/atlas/interpolation/method/knn/GridBoxMaximum.cc @@ -22,6 +22,7 @@ #include "atlas/functionspace/PointCloud.h" #include "atlas/interpolation/method/MethodFactory.h" #include "atlas/runtime/Exception.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" namespace atlas { @@ -60,8 +61,8 @@ void GridBoxMaximum::do_execute(const Field& source, Field& target, Metadata&) c if (!matrixFree_) { - const Matrix& m = matrix(); - Matrix::const_iterator k(m); + const auto m = atlas::linalg::make_non_owning_eckit_sparse_matrix(matrix()); + auto k = m.begin(); for (decltype(m.rows()) i = 0, j = 0; i < m.rows(); ++i) { double max = std::numeric_limits::lowest(); diff --git a/src/atlas/interpolation/method/knn/GridBoxMethod.cc b/src/atlas/interpolation/method/knn/GridBoxMethod.cc index 729146e20..64e70e3c1 100644 --- a/src/atlas/interpolation/method/knn/GridBoxMethod.cc +++ b/src/atlas/interpolation/method/knn/GridBoxMethod.cc @@ -177,8 +177,7 @@ void GridBoxMethod::do_setup(const Grid& source, const Grid& target, const Cache { ATLAS_TRACE("GridBoxMethod::setup: build interpolant matrix"); - Matrix A(targetBoxes_.size(), sourceBoxes_.size(), allTriplets); - setMatrix(A); + setMatrix(targetBoxes_.size(), sourceBoxes_.size(), allTriplets); } } diff --git a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc index 4ac350962..fe36c8139 100644 --- a/src/atlas/interpolation/method/knn/KNearestNeighbours.cc +++ b/src/atlas/interpolation/method/knn/KNearestNeighbours.cc @@ -130,8 +130,7 @@ void KNearestNeighbours::do_setup(const FunctionSpace& source, const FunctionSpa } // fill sparse matrix and return - Matrix A(out_npts, inp_npts, weights_triplets); - setMatrix(A); + setMatrix(out_npts, inp_npts, weights_triplets); } } // namespace method diff --git a/src/atlas/interpolation/method/knn/NearestNeighbour.cc b/src/atlas/interpolation/method/knn/NearestNeighbour.cc index 65f36a46d..56b9432aa 100644 --- a/src/atlas/interpolation/method/knn/NearestNeighbour.cc +++ b/src/atlas/interpolation/method/knn/NearestNeighbour.cc @@ -102,8 +102,7 @@ void NearestNeighbour::do_setup(const FunctionSpace& source, const FunctionSpace } // fill sparse matrix and return - Matrix A(out_npts, inp_npts, weights_triplets); - setMatrix(A); + setMatrix( out_npts, inp_npts, weights_triplets ); } } // namespace method diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 8bcd0bef1..919547139 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -18,6 +18,7 @@ #include "atlas/interpolation/method/MethodFactory.h" #include "atlas/interpolation/method/sphericalvector/ComplexMatrixMultiply.h" #include "atlas/interpolation/method/sphericalvector/Types.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" #include "atlas/option/Options.h" #include "atlas/parallel/omp/omp.h" #include "atlas/runtime/Exception.h" @@ -65,12 +66,13 @@ void SphericalVector::do_setup(const FunctionSpace& source, setMatrix(Interpolation(interpolationScheme_, source_, target_)); // Get matrix data. - const auto nRows = static_cast(matrix().rows()); - const auto nCols = static_cast(matrix().cols()); - const auto nNonZeros = static_cast(matrix().nonZeros()); - const auto* outerIndices = matrix().outer(); - const auto* innerIndices = matrix().inner(); - const auto* baseWeights = matrix().data(); + const auto m = atlas::linalg::make_non_owning_eckit_sparse_matrix(matrix()); + const auto nRows = static_cast(m.rows()); + const auto nCols = static_cast(m.cols()); + const auto nNonZeros = static_cast(m.nonZeros()); + const auto* outerIndices = m.outer(); + const auto* innerIndices = m.inner(); + const auto* baseWeights = m.data(); // Note: need to store copy of weights as Eigen3 sorts compressed rows by j // whereas eckit does not. diff --git a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc index 2f2fa75a8..887f57a2c 100644 --- a/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc +++ b/src/atlas/interpolation/method/structured/StructuredInterpolation2D.tcc @@ -428,8 +428,7 @@ void StructuredInterpolation2D::setup( const FunctionSpace& source ) { // fill sparse matrix if( failed_points.empty() && out_npts_) { idx_t inp_npts = source.size(); - Matrix A( out_npts_, inp_npts, triplets ); - setMatrix(A); + setMatrix(out_npts_, inp_npts, triplets); } } } diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc index d1ea5addf..144c4fcd5 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.cc @@ -870,13 +870,11 @@ void ConservativeSphericalPolygonInterpolation::do_setup(const FunctionSpace& sr stopwatch.start(); switch (order_) { case 1: { - auto M = compute_1st_order_matrix(); - setMatrix(M, "1"); + setMatrix(n_tpoints_, n_spoints_, compute_1st_order_triplets(), "1"); break; } case 2: { - auto M = compute_2nd_order_matrix(); - setMatrix(M, "2"); + setMatrix(n_tpoints_, n_spoints_, compute_2nd_order_triplets(), "2"); break; } default: { @@ -1248,7 +1246,7 @@ void ConservativeSphericalPolygonInterpolation::intersect_polygons(const CSPolyg } } -eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1st_order_matrix() { +ConservativeSphericalPolygonInterpolation::Triplets ConservativeSphericalPolygonInterpolation::compute_1st_order_triplets() { ATLAS_TRACE("ConservativeMethod::setup: build cons-1 interpolant matrix"); ATLAS_ASSERT(not matrix_free_); Triplets triplets; @@ -1355,10 +1353,10 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_1 // } // } // } - return Matrix(n_tpoints_, n_spoints_, triplets); + return triplets; } -eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2nd_order_matrix() { +ConservativeSphericalPolygonInterpolation::Triplets ConservativeSphericalPolygonInterpolation::compute_2nd_order_triplets() { ATLAS_TRACE("ConservativeMethod::setup: build cons-2 interpolant matrix"); ATLAS_ASSERT(not matrix_free_); const auto& src_points_ = data_->src_points_; @@ -1562,7 +1560,7 @@ eckit::linalg::SparseMatrix ConservativeSphericalPolygonInterpolation::compute_2 } } sort_and_accumulate_triplets(triplets); // Very expensive!!! (90% of this routine). We need to avoid it - return Matrix(n_tpoints_, n_spoints_, triplets); + return triplets; } void ConservativeSphericalPolygonInterpolation::do_execute(const FieldSet& src_fields, FieldSet& tgt_fields, diff --git a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h index 6d4db4dba..9f5deff94 100644 --- a/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h +++ b/src/atlas/interpolation/method/unstructured/ConservativeSphericalPolygonInterpolation.h @@ -159,8 +159,8 @@ class ConservativeSphericalPolygonInterpolation : public Method { void do_setup_impl(const Grid& src_grid, const Grid& tgt_grid); void intersect_polygons(const CSPolygonArray& src_csp, const CSPolygonArray& tgt_scp); - Matrix compute_1st_order_matrix(); - Matrix compute_2nd_order_matrix(); + Triplets compute_1st_order_triplets(); + Triplets compute_2nd_order_triplets(); void dump_intersection(const std::string, const ConvexSphericalPolygon& plg_1, const CSPolygonArray& plg_2_array, const std::vector& plg_2_idx_array) const; template diff --git a/src/atlas/interpolation/method/unstructured/FiniteElement.cc b/src/atlas/interpolation/method/unstructured/FiniteElement.cc index 3e41003da..740ad7c12 100644 --- a/src/atlas/interpolation/method/unstructured/FiniteElement.cc +++ b/src/atlas/interpolation/method/unstructured/FiniteElement.cc @@ -25,6 +25,7 @@ #include "atlas/interpolation/element/Triag3D.h" #include "atlas/interpolation/method/MethodFactory.h" #include "atlas/interpolation/method/Ray.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" #include "atlas/mesh/ElementType.h" #include "atlas/mesh/Nodes.h" #include "atlas/mesh/actions/BuildCellCentres.h" @@ -152,7 +153,8 @@ void FiniteElement::print(std::ostream& out) const { auto stencil_size_loc = array::make_view(field_stencil_size_loc); stencil_size_loc.assign(0); - for (auto it = matrix().begin(); it != matrix().end(); ++it) { + const auto m = atlas::linalg::make_non_owning_eckit_sparse_matrix(matrix()); + for (auto it = m.begin(); it != m.end(); ++it) { idx_t p = idx_t(it.row()); idx_t& i = stencil_size_loc(p); stencil_points_loc(p, i) = gidx_src(it.col()); @@ -331,8 +333,7 @@ void FiniteElement::setup(const FunctionSpace& source) { } // fill sparse matrix and return - Matrix A(out_npts, inp_npts, weights_triplets); - setMatrix(A); + setMatrix(out_npts, inp_npts, weights_triplets); } struct ElementEdge { diff --git a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc index 840c851a8..6285cdde4 100644 --- a/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc +++ b/src/atlas/interpolation/method/unstructured/UnstructuredBilinearLonLat.cc @@ -22,6 +22,7 @@ #include "atlas/interpolation/element/Triag2D.h" #include "atlas/interpolation/method/MethodFactory.h" #include "atlas/interpolation/method/Ray.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" #include "atlas/mesh/ElementType.h" #include "atlas/mesh/Nodes.h" #include "atlas/mesh/actions/BuildCellCentres.h" @@ -148,7 +149,8 @@ void UnstructuredBilinearLonLat::print(std::ostream& out) const { auto stencil_size_loc = array::make_view(field_stencil_size_loc); stencil_size_loc.assign(0); - for (auto it = matrix().begin(); it != matrix().end(); ++it) { + const auto m = atlas::linalg::make_non_owning_eckit_sparse_matrix(matrix()); + for (auto it = m.begin(); it != m.end(); ++it) { idx_t p = idx_t(it.row()); idx_t& i = stencil_size_loc(p); stencil_points_loc(p, i) = gidx_src(it.col()); @@ -300,8 +302,7 @@ void UnstructuredBilinearLonLat::setup(const FunctionSpace& source) { } // fill sparse matrix and return - Matrix A(out_npts, inp_npts, weights_triplets); - setMatrix(A); + setMatrix(out_npts, inp_npts, weights_triplets); } Method::Triplets UnstructuredBilinearLonLat::projectPointToElements(size_t ip, const ElemIndex3::NodeList& elems, diff --git a/src/atlas/linalg/sparse.h b/src/atlas/linalg/sparse.h index b9398716a..4c6ebc02e 100644 --- a/src/atlas/linalg/sparse.h +++ b/src/atlas/linalg/sparse.h @@ -10,7 +10,11 @@ #pragma once #include "sparse/Backend.h" +#include "sparse/MakeEckitSparseMatrix.h" +#include "sparse/MakeSparseMatrixStorageEckit.h" #include "sparse/SparseMatrixMultiply.h" +#include "sparse/SparseMatrixStorage.h" +#include "sparse/SparseMatrixView.h" namespace atlas { namespace linalg {} // namespace linalg diff --git a/src/atlas/linalg/sparse/MakeEckitSparseMatrix.h b/src/atlas/linalg/sparse/MakeEckitSparseMatrix.h new file mode 100644 index 000000000..a639c649d --- /dev/null +++ b/src/atlas/linalg/sparse/MakeEckitSparseMatrix.h @@ -0,0 +1,85 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include + +#include "eckit/linalg/SparseMatrix.h" + +#include "atlas/linalg/sparse/SparseMatrixView.h" +#include "atlas/linalg/sparse/SparseMatrixStorage.h" + +namespace atlas { +namespace linalg { + +//---------------------------------------------------------------------------------------------------------------------- + +class EckitSparseMatrixNonOwningAllocator : public eckit::linalg::SparseMatrix::Allocator { +public: + EckitSparseMatrixNonOwningAllocator(const SparseMatrixView& view) : view_{view} {} + + eckit::linalg::SparseMatrix::Layout allocate(eckit::linalg::SparseMatrix::Shape& shape) override { + eckit::linalg::SparseMatrix::Layout p; + using Index = eckit::linalg::Index; + using OuterIndex = std::remove_pointer_t; + using InnerIndex = std::remove_pointer_t; + static_assert(sizeof(OuterIndex) == sizeof(Index)); + static_assert(sizeof(InnerIndex) == sizeof(Index)); + p.data_ = const_cast(view_.value()); + p.outer_ = reinterpret_cast(const_cast(view_.outer())); + p.inner_ = reinterpret_cast(const_cast(view_.inner())); + shape.size_ = view_.nnz(); + shape.rows_ = view_.rows(); + shape.cols_ = view_.cols(); + return p; + } + + void deallocate(eckit::linalg::SparseMatrix::Layout p, eckit::linalg::SparseMatrix::Shape) override { + std::cout << "do nothing" << std::endl; + } + + bool inSharedMemory() const override { return false; } + + void print(std::ostream& out) const override { out << "atlas::linalg::EckitSparseMatrixNonOwningAllocator[rows="< view_; +}; + +//---------------------------------------------------------------------------------------------------------------------- + +// Create new eckit sparsematrix from a SparseMatrixView. Warning, this matrix does not own the data! +inline eckit::linalg::SparseMatrix make_non_owning_eckit_sparse_matrix(const SparseMatrixView& view) { + eckit::linalg::SparseMatrix m(new EckitSparseMatrixNonOwningAllocator(view)); + return m; +} + +// Create new eckit sparsematrix from a SparseMatrixStorage. Warning, this matrix does not own the data! +inline eckit::linalg::SparseMatrix make_non_owning_eckit_sparse_matrix(const SparseMatrixStorage& m) { + return make_non_owning_eckit_sparse_matrix( + make_host_view(m) ); +} + +inline eckit::linalg::SparseMatrix make_eckit_sparse_matrix(const SparseMatrixView& view) { + eckit::linalg::SparseMatrix m(new EckitSparseMatrixNonOwningAllocator(view)); + return eckit::linalg::SparseMatrix{m}; // makes copy +} + +// Create new eckit sparsematrix from a SparseMatrixStorage. Warning, this creates a copy! +inline eckit::linalg::SparseMatrix make_eckit_sparse_matrix(const SparseMatrixStorage& m) { + return make_eckit_sparse_matrix( + make_host_view(m) ); +} + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace linalg +} // namespace atlas diff --git a/src/atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h b/src/atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h new file mode 100644 index 000000000..fffbf9eeb --- /dev/null +++ b/src/atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h @@ -0,0 +1,135 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include +#include +#include +#include + +#include "eckit/linalg/SparseMatrix.h" + +#include "atlas/array.h" +#include "atlas/linalg/sparse/SparseMatrixStorage.h" + +namespace atlas { +namespace linalg { + +//---------------------------------------------------------------------------------------------------------------------- + +template +void host_copy_eckit (const InputT* input_data, array::Array& output) { + auto size = output.size(); + OutputT* output_data = output.host_data(); + std::copy( input_data, input_data + size, output_data ); +} + +inline SparseMatrixStorage make_sparse_matrix_storage(eckit::linalg::SparseMatrix&& m) { + using Index = eckit::linalg::Index; + using Value = eckit::linalg::Scalar; + std::cout << "make_sparse_matrix_storage move" << std::endl; + + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = array::Array::wrap(const_cast(m.outer()), array::make_shape(rows+1)); + auto* inner = array::Array::wrap(const_cast(m.inner()), array::make_shape(nnz)); + auto* value = array::Array::wrap(const_cast(m.data()), array::make_shape(nnz)); + + // We now move the eckit::linalg::SparseMatrix into a generic storage so + // the wrapped array data does not go out of scope + auto storage = std::make_any(std::move(m)); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::move(storage) + ); +} + +template< typename value_type, typename index_type = eckit::linalg::Index> +SparseMatrixStorage make_sparse_matrix_storage(eckit::linalg::SparseMatrix&& m) { + std::cout << "make_sparse_matrix_storage move or convert" << std::endl; + + if (std::is_same_v && std::is_same_v) { + return make_sparse_matrix_storage(std::move(m)); + } + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = array::Array::create(rows+1); + auto* inner = array::Array::create(nnz); + auto* value = array::Array::create(nnz); + + host_copy_eckit(m.outer(), *outer); + host_copy_eckit(m.inner(), *inner); + host_copy_eckit(m.data(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + + +inline SparseMatrixStorage make_sparse_matrix_storage(const eckit::linalg::SparseMatrix& m) { + using Index = eckit::linalg::Index; + using Value = eckit::linalg::Scalar; + + std::cout << "make_sparse_matrix_storage copy" << std::endl; + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = array::Array::create(rows+1); + auto* inner = array::Array::create(nnz); + auto* value = array::Array::create(nnz); + + host_copy_eckit(m.outer(), *outer); + host_copy_eckit(m.inner(), *inner); + host_copy_eckit(m.data(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + +template< typename value_type, typename index_type = eckit::linalg::Index> +SparseMatrixStorage make_sparse_matrix_storage(const eckit::linalg::SparseMatrix& m) { + std::cout << "make_sparse_matrix_storage copy and convert" << std::endl; + + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = array::Array::create(rows+1); + auto* inner = array::Array::create(nnz); + auto* value = array::Array::create(nnz); + host_copy_eckit(m.outer(), *outer); + host_copy_eckit(m.inner(), *inner); + host_copy_eckit(m.data(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace linalg +} // namespace atlas \ No newline at end of file diff --git a/src/atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h b/src/atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h new file mode 100644 index 000000000..f306e85f0 --- /dev/null +++ b/src/atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h @@ -0,0 +1,142 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include +#include +#include +#include + +#include "atlas/library/defines.h" +#if ATLAS_HAVE_EIGEN +#include +#endif + +#include "atlas/array.h" +#include "atlas/linalg/sparse/SparseMatrixStorage.h" + +namespace atlas { +namespace linalg { + +//---------------------------------------------------------------------------------------------------------------------- + +#if ATLAS_HAVE_EIGEN + +template +void host_copy_eigen (const InputT* input_data, atlas::array::Array& output) { + auto size = output.size(); + OutputT* output_data = output.host_data(); + std::copy( input_data, input_data + size, output_data ); +} + +template +SparseMatrixStorage make_sparse_matrix_storage(Eigen::SparseMatrix&& m) { + std::cout << "make_sparse_matrix_storage move" << std::endl; + + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = atlas::array::Array::wrap(const_cast(m.outerIndexPtr()), atlas::array::make_shape(rows+1)); + auto* inner = atlas::array::Array::wrap(const_cast(m.innerIndexPtr()), atlas::array::make_shape(nnz)); + auto* value = atlas::array::Array::wrap(const_cast(m.valuePtr()), atlas::array::make_shape(nnz)); + + using EigenMatrix = Eigen::SparseMatrix; + auto m_ptr = std::make_shared(); + m_ptr->swap(m); + auto storage = std::make_any>(std::move(m_ptr)); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::move(storage) + ); +} + +template< typename value_type, typename index_type = eckit::linalg::Index, typename Value, typename Index, + typename = std::enable_if_t < !std::is_same_v>> +SparseMatrixStorage make_sparse_matrix_storage(Eigen::SparseMatrix&& m) { + std::cout << "make_sparse_matrix_storage move or convert" << std::endl; + + if (std::is_same_v && std::is_same_v) { + return make_sparse_matrix_storage(std::move(m)); + } + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = atlas::array::Array::create(rows+1); + auto* inner = atlas::array::Array::create(nnz); + auto* value = atlas::array::Array::create(nnz); + + host_copy_eigen(m.outerIndexPtr(), *outer); + host_copy_eigen(m.innerIndexPtr(), *inner); + host_copy_eigen(m.valuePtr(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + + +template +SparseMatrixStorage make_sparse_matrix_storage(const Eigen::SparseMatrix& m) { + std::cout << "make_sparse_matrix_storage copy" << std::endl; + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = atlas::array::Array::create(rows+1); + auto* inner = atlas::array::Array::create(nnz); + auto* value = atlas::array::Array::create(nnz); + + host_copy_eigen(m.outerIndexPtr(), *outer); + host_copy_eigen(m.innerIndexPtr(), *inner); + host_copy_eigen(m.valuePtr(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + +template< typename value_type, typename index_type = eckit::linalg::Index, typename Value, typename Index, + typename = std::enable_if_t < !std::is_same_v>> +SparseMatrixStorage make_sparse_matrix_storage(const Eigen::SparseMatrix& m) { + std::cout << "make_sparse_matrix_storage copy and convert" << std::endl; + + std::size_t rows = m.rows(); + std::size_t cols = m.cols(); + std::size_t nnz = m.nonZeros(); + auto* outer = atlas::array::Array::create(rows+1); + auto* inner = atlas::array::Array::create(nnz); + auto* value = atlas::array::Array::create(nnz); + host_copy_eigen(m.outerIndexPtr(), *outer); + host_copy_eigen(m.innerIndexPtr(), *inner); + host_copy_eigen(m.valuePtr(), *value); + + return SparseMatrixStorage::make( rows, cols, nnz, + std::unique_ptr(value), + std::unique_ptr(inner), + std::unique_ptr(outer), + std::any() + ); +} + +#endif + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace linalg +} // namespace atlas \ No newline at end of file diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply.h b/src/atlas/linalg/sparse/SparseMatrixMultiply.h index 9e337fe08..9f5823bdd 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply.h +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply.h @@ -11,18 +11,17 @@ #pragma once #include "eckit/config/Configuration.h" -#include "eckit/linalg/SparseMatrix.h" #include "atlas/linalg/Indexing.h" #include "atlas/linalg/View.h" #include "atlas/linalg/sparse/Backend.h" +#include "atlas/linalg/sparse/SparseMatrixView.h" #include "atlas/runtime/Exception.h" #include "atlas/util/Config.h" namespace atlas { namespace linalg { -using SparseMatrix = eckit::linalg::SparseMatrix; using Configuration = eckit::Configuration; template @@ -96,13 +95,13 @@ class SparseMatrixMultiply { namespace sparse { // Template class which needs (full or partial) specialization for concrete template parameters -template +template struct SparseMatrixMultiply { - static void multiply(const SparseMatrix&, const View&, View&, + static void multiply(const SparseMatrixView&, const View&, View&, const Configuration&) { throw_NotImplemented("SparseMatrixMultiply needs a template specialization with the implementation", Here()); } - static void multiply_add(const SparseMatrix&, const View&, View&, + static void multiply_add(const SparseMatrixView&, const View&, View&, const Configuration&) { throw_NotImplemented("SparseMatrixMultiply needs a template specialization with the implementation", Here()); } diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc b/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc index 4b44b984c..a41ad700d 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply.tcc @@ -12,10 +12,13 @@ #include "SparseMatrixMultiply.h" +#include "eckit/linalg/SparseMatrix.h" + #include "atlas/linalg/Indexing.h" #include "atlas/linalg/Introspection.h" #include "atlas/linalg/View.h" #include "atlas/linalg/sparse/Backend.h" +#include "atlas/linalg/sparse/SparseMatrixView.h" #include "atlas/runtime/Exception.h" #if ATLAS_ECKIT_HAVE_ECKIT_585 @@ -30,28 +33,67 @@ namespace linalg { namespace sparse { namespace { + +inline SparseMatrixView make_view_from_eckit_sparse_matrix(const eckit::linalg::SparseMatrix& m) { + return SparseMatrixView{ + m.rows(), + m.cols(), + m.nonZeros(), + m.data(), + m.inner(), + m.outer() + }; +} + template struct SparseMatrixMultiplyHelper { + template + static void multiply( const Matrix& W, const SourceView& src, TargetView& tgt, + const eckit::Configuration& config ) { + using MatrixValue = typename std::remove_const::type; + using SourceValue = const typename std::remove_const::type; + using TargetValue = typename std::remove_const::type; + constexpr int src_rank = introspection::rank(); + constexpr int tgt_rank = introspection::rank(); + static_assert( src_rank == tgt_rank, "src and tgt need same rank" ); + SparseMatrixMultiply::multiply( W, src, tgt, config ); + } + + template + static void multiply_add( const Matrix& W, const SourceView& src, TargetView& tgt, + const eckit::Configuration& config ) { + using MatrixValue = typename std::remove_const::type; + using SourceValue = const typename std::remove_const::type; + using TargetValue = typename std::remove_const::type; + constexpr int src_rank = introspection::rank(); + constexpr int tgt_rank = introspection::rank(); + static_assert( src_rank == tgt_rank, "src and tgt need same rank" ); + SparseMatrixMultiply::multiply_add( W, src, tgt, config ); + } template - static void multiply( const SparseMatrix& W, const SourceView& src, TargetView& tgt, + static void multiply( const eckit::linalg::SparseMatrix& W, const SourceView& src, TargetView& tgt, const eckit::Configuration& config ) { + using MatrixValue = double; + using MatrixIndex = int; using SourceValue = const typename std::remove_const::type; using TargetValue = typename std::remove_const::type; constexpr int src_rank = introspection::rank(); constexpr int tgt_rank = introspection::rank(); static_assert( src_rank == tgt_rank, "src and tgt need same rank" ); - SparseMatrixMultiply::multiply( W, src, tgt, config ); + SparseMatrixMultiply::multiply( make_view_from_eckit_sparse_matrix(W), src, tgt, config ); } template - static void multiply_add( const SparseMatrix& W, const SourceView& src, TargetView& tgt, - const eckit::Configuration& config ) { + static void multiply_add( const eckit::linalg::SparseMatrix& W, const SourceView& src, TargetView& tgt, + const eckit::Configuration& config ) { + using MatrixValue = double; + using MatrixIndex = int; using SourceValue = const typename std::remove_const::type; using TargetValue = typename std::remove_const::type; constexpr int src_rank = introspection::rank(); constexpr int tgt_rank = introspection::rank(); static_assert( src_rank == tgt_rank, "src and tgt need same rank" ); - SparseMatrixMultiply::multiply_add( W, src, tgt, config ); + SparseMatrixMultiply::multiply_add( make_view_from_eckit_sparse_matrix(W), src, tgt, config ); } }; diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc index 4c4ec1c58..a332f7cdd 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.cc @@ -1,15 +1,5 @@ /* - * (C) Copyright 2013 ECMWF. - * - * This software is licensed under the terms of the Apache Licence Version 2.0 - * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - * In applying this licence, ECMWF does not waive the privileges and immunities - * granted to it by virtue of its status as an intergovernmental organisation - * nor does it submit to any jurisdiction. - */ - -/* - * (C) Copyright 2013 ECMWF. + * (C) Copyright 2024- ECMWF. * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -33,6 +23,7 @@ #include "eckit/linalg/Vector.h" #include "atlas/runtime/Exception.h" +#include "atlas/linalg/sparse/MakeEckitSparseMatrix.h" namespace atlas { namespace linalg { @@ -71,54 +62,56 @@ auto linalg_make_view(atlas::array::ArrayT& array) { } // namespace -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { ATLAS_ASSERT(src.contiguous()); ATLAS_ASSERT(tgt.contiguous()); eckit::linalg::Vector v_src(src.data(), src.size()); eckit::linalg::Vector v_tgt(tgt.data(), tgt.size()); - eckit_linalg_backend(config).spmv(W, v_src, v_tgt); + eckit::linalg::SparseMatrix m_W = make_non_owning_eckit_sparse_matrix(W); + eckit_linalg_backend(config).spmv(m_W, v_src, v_tgt); } -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { ATLAS_ASSERT(src.contiguous()); ATLAS_ASSERT(tgt.contiguous()); ATLAS_ASSERT(src.shape(1) >= W.cols()); ATLAS_ASSERT(tgt.shape(1) >= W.rows()); eckit::linalg::Matrix m_src(src.data(), src.shape(1), src.shape(0)); eckit::linalg::Matrix m_tgt(tgt.data(), tgt.shape(1), tgt.shape(0)); - eckit_linalg_backend(config).spmm(W, m_src, m_tgt); + eckit::linalg::SparseMatrix m_W = make_non_owning_eckit_sparse_matrix(W); + eckit_linalg_backend(config).spmm(m_W, m_src, m_tgt); } -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { - SparseMatrixMultiply::multiply(W, src, tgt, +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { + SparseMatrixMultiply::multiply(W, src, tgt, config); } -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { array::ArrayT tmp(src.shape(0)); auto v_tmp = linalg_make_view(tmp); v_tmp.assign(0.); - SparseMatrixMultiply::multiply(W, src, v_tmp, config); + SparseMatrixMultiply::multiply(W, src, v_tmp, config); for (idx_t t = 0; t < tmp.shape(0); ++t) { tgt(t) += v_tmp(t); } } -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { array::ArrayT tmp(src.shape(0), src.shape(1)); auto v_tmp = linalg_make_view(tmp); v_tmp.assign(0.); - SparseMatrixMultiply::multiply(W, src, v_tmp, config); + SparseMatrixMultiply::multiply(W, src, v_tmp, config); for (idx_t t = 0; t < tmp.shape(0); ++t) { for (idx_t k = 0; k < tmp.shape(1); ++k) { @@ -127,10 +120,9 @@ void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { - SparseMatrixMultiply::multiply_add(W, src, tgt, - config); +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { + SparseMatrixMultiply::multiply_add(W, src, tgt, config); } } // namespace sparse diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.h b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.h index 6915a10f3..2f470e786 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.h +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply_EckitLinalg.h @@ -1,5 +1,5 @@ /* - * (C) Copyright 2013 ECMWF. + * (C) Copyright 2024- ECMWF. * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -8,15 +8,6 @@ * nor does it submit to any jurisdiction. */ -/* - * (C) Copyright 2013 ECMWF. - * - * This software is licensed under the terms of the Apache Licence Version 2.0 - * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. - * In applying this licence, ECMWF does not waive the privileges and immunities - * granted to it by virtue of its status as an intergovernmental organisation - * nor does it submit to any jurisdiction. - */ #pragma once @@ -27,28 +18,28 @@ namespace linalg { namespace sparse { template <> -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); - static void multiply_add(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); + static void multiply_add(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); }; template <> -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); - static void multiply_add(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); + static void multiply_add(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); }; template <> -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); - static void multiply_add(const SparseMatrix&, const View& src, View& tgt, - const Configuration&); +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); + static void multiply_add(const SparseMatrixView&, const View& src, View& tgt, + const Configuration&); }; } // namespace sparse diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.cc b/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.cc index 41d01737f..b4a96c638 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.cc +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.cc @@ -17,12 +17,12 @@ namespace atlas { namespace linalg { namespace sparse { -template -void spmv_layout_left(const SparseMatrix& W, const View& src, View& tgt) { +template +void spmv_layout_left(const SparseMatrixView& W, const View& src, View& tgt) { using Value = TargetValue; const auto outer = W.outer(); const auto index = W.inner(); - const auto weight = W.data(); + const auto weight = W.value(); const idx_t rows = static_cast(W.rows()); ATLAS_ASSERT(src.shape(0) >= W.cols()); @@ -40,24 +40,24 @@ void spmv_layout_left(const SparseMatrix& W, const View& src, Vi } } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmv_layout_left(W, src, tgt); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmv_layout_left(W, src, tgt); } -template -void spmm_layout_left(const SparseMatrix& W, const View& src, View& tgt) { +template +void spmm_layout_left(const SparseMatrixView& W, const View& src, View& tgt) { using Value = TargetValue; const auto outer = W.outer(); const auto index = W.inner(); - const auto weight = W.data(); + const auto weight = W.value(); const idx_t rows = static_cast(W.rows()); const idx_t Nk = src.shape(1); @@ -80,20 +80,20 @@ void spmm_layout_left(const SparseMatrix& W, const View& src, Vi } } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmm_layout_left(W, src, tgt); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmm_layout_left(W, src, tgt); } -template -void spmt_layout_left(const SparseMatrix& W, const View& src, View& tgt) { +template +void spmt_layout_left(const SparseMatrixView& W, const View& src, View& tgt) { if (src.contiguous() && tgt.contiguous()) { // We can take a more optimized route by reducing rank auto src_v = View(src.data(), array::make_shape(src.shape(0), src.stride(0))); @@ -104,7 +104,7 @@ void spmt_layout_left(const SparseMatrix& W, const View& src, Vi using Value = TargetValue; const auto outer = W.outer(); const auto index = W.inner(); - const auto weight = W.data(); + const auto weight = W.value(); const idx_t rows = static_cast(W.rows()); const idx_t Nk = src.shape(1); const idx_t Nl = src.shape(2); @@ -129,36 +129,36 @@ void spmt_layout_left(const SparseMatrix& W, const View& src, Vi } } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { spmt_layout_left(W, src, tgt); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { spmt_layout_left(W, src, tgt); } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { - return SparseMatrixMultiply::multiply(W, src, tgt, config); +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { + return SparseMatrixMultiply::multiply(W, src, tgt, config); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { - return SparseMatrixMultiply::multiply_add(W, src, tgt, config); +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { + return SparseMatrixMultiply::multiply_add(W, src, tgt, config); } -template -void spmm_layout_right(const SparseMatrix& W, const View& src, View& tgt) { +template +void spmm_layout_right(const SparseMatrixView& W, const View& src, View& tgt) { using Value = TargetValue; const auto outer = W.outer(); const auto index = W.inner(); - const auto weight = W.data(); + const auto weight = W.value(); const idx_t rows = static_cast(W.rows()); const idx_t Nk = src.shape(0); @@ -181,20 +181,20 @@ void spmm_layout_right(const SparseMatrix& W, const View& src, V } } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmm_layout_right(W, src, tgt); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration&) { +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration&) { spmm_layout_right(W, src, tgt); } -template -void spmt_layout_right(const SparseMatrix& W, const View& src, View& tgt) { +template +void spmt_layout_right(const SparseMatrixView& W, const View& src, View& tgt) { if (src.contiguous() && tgt.contiguous()) { // We can take a more optimized route by reducing rank auto src_v = View(src.data(), array::make_shape(src.shape(0), src.stride(0))); @@ -205,7 +205,7 @@ void spmt_layout_right(const SparseMatrix& W, const View& src, V using Value = TargetValue; const auto outer = W.outer(); const auto index = W.inner(); - const auto weight = W.data(); + const auto weight = W.value(); const idx_t rows = static_cast(W.rows()); const idx_t Nk = src.shape(1); const idx_t Nl = src.shape(0); @@ -230,25 +230,31 @@ void spmt_layout_right(const SparseMatrix& W, const View& src, V } } -template -void SparseMatrixMultiply::multiply( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +template +void SparseMatrixMultiply::multiply( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { spmt_layout_right(W, src, tgt); } -template -void SparseMatrixMultiply::multiply_add( - const SparseMatrix& W, const View& src, View& tgt, const Configuration& config) { +template +void SparseMatrixMultiply::multiply_add( + const SparseMatrixView& W, const View& src, View& tgt, const Configuration& config) { spmt_layout_right(W, src, tgt); } -#define EXPLICIT_TEMPLATE_INSTANTIATION(TYPE) \ - template struct SparseMatrixMultiply; \ - template struct SparseMatrixMultiply; \ - template struct SparseMatrixMultiply; \ - template struct SparseMatrixMultiply; \ - template struct SparseMatrixMultiply; \ - template struct SparseMatrixMultiply; +#define EXPLICIT_TEMPLATE_INSTANTIATION(TYPE) \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; \ + template struct SparseMatrixMultiply; EXPLICIT_TEMPLATE_INSTANTIATION(double); EXPLICIT_TEMPLATE_INSTANTIATION(float); diff --git a/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.h b/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.h index a3bdcea44..590dc17b3 100644 --- a/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.h +++ b/src/atlas/linalg/sparse/SparseMatrixMultiply_OpenMP.h @@ -17,51 +17,51 @@ namespace linalg { namespace sparse { -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; -template -struct SparseMatrixMultiply { - static void multiply(const SparseMatrix& W, const View& src, View& tgt, +template +struct SparseMatrixMultiply { + static void multiply(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); - static void multiply_add(const SparseMatrix& W, const View& src, View& tgt, + static void multiply_add(const SparseMatrixView& W, const View& src, View& tgt, const Configuration&); }; diff --git a/src/atlas/linalg/sparse/SparseMatrixStorage.cc b/src/atlas/linalg/sparse/SparseMatrixStorage.cc new file mode 100644 index 000000000..2e956558a --- /dev/null +++ b/src/atlas/linalg/sparse/SparseMatrixStorage.cc @@ -0,0 +1,131 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "SparseMatrixStorage.h" + +namespace atlas::linalg { + + +SparseMatrixStorage::SparseMatrixStorage(SparseMatrixStorage&& other) { + nnz_ = other.nnz_; + rows_ = other.rows_; + cols_ = other.cols_; + outer_ = std::move(other.outer_); + inner_ = std::move(other.inner_); + value_ = std::move(other.value_); + storage_ = std::move(other.storage_); + + // Invalidate other, just in case it is being used (IT SHOULDNT) + other.nnz_ = 0; + other.rows_ = 0; + other.cols_ = 0; +} + +SparseMatrixStorage::SparseMatrixStorage(const SparseMatrixStorage& other) { + nnz_ = other.nnz_; + rows_ = other.rows_; + cols_ = other.cols_; + outer_.reset(atlas::array::Array::create(other.outer_->datatype(), atlas::array::make_shape(other.outer_->size()))); + inner_.reset(atlas::array::Array::create(other.inner_->datatype(), atlas::array::make_shape(other.inner_->size()))); + value_.reset(atlas::array::Array::create(other.value_->datatype(), atlas::array::make_shape(other.value_->size()))); + outer_->copy(*other.outer_); + inner_->copy(*other.inner_); + value_->copy(*other.value_); +} + +SparseMatrixStorage& SparseMatrixStorage::operator=(SparseMatrixStorage&& other) { + nnz_ = other.nnz_; + rows_ = other.rows_; + cols_ = other.cols_; + outer_ = std::move(other.outer_); + inner_ = std::move(other.inner_); + value_ = std::move(other.value_); + storage_ = std::move(other.storage_); + + // Invalidate other, just in case it is being used (IT SHOULDNT) + other.nnz_ = 0; + other.rows_ = 0; + other.cols_ = 0; + return *this; +} + + +SparseMatrixStorage& SparseMatrixStorage::operator=(const SparseMatrixStorage& other) { + nnz_ = other.nnz_; + rows_ = other.rows_; + cols_ = other.cols_; + outer_.reset(atlas::array::Array::create(other.outer_->datatype(), atlas::array::make_shape(other.outer_->size()))); + inner_.reset(atlas::array::Array::create(other.inner_->datatype(), atlas::array::make_shape(other.inner_->size()))); + value_.reset(atlas::array::Array::create(other.value_->datatype(), atlas::array::make_shape(other.value_->size()))); + outer_->copy(*other.outer_); + inner_->copy(*other.inner_); + value_->copy(*other.value_); + return *this; +} + + +void SparseMatrixStorage::updateDevice() const { + outer_->updateDevice(); + inner_->updateDevice(); + value_->updateDevice(); +} + +void SparseMatrixStorage::updateHost() const { + outer_->updateHost(); + inner_->updateHost(); + value_->updateHost(); +} + +bool SparseMatrixStorage::hostNeedsUpdate() const { + return outer_->hostNeedsUpdate() || + inner_->hostNeedsUpdate() || + value_->hostNeedsUpdate(); +} + +bool SparseMatrixStorage::deviceNeedsUpdate() const { + return outer_->deviceNeedsUpdate() || + inner_->deviceNeedsUpdate() || + value_->deviceNeedsUpdate(); +} + +void SparseMatrixStorage::setHostNeedsUpdate(bool v) const { + outer_->setHostNeedsUpdate(v); + inner_->setHostNeedsUpdate(v); + value_->setHostNeedsUpdate(v); +} + +void SparseMatrixStorage::setDeviceNeedsUpdate(bool v) const { + outer_->setDeviceNeedsUpdate(v); + inner_->setDeviceNeedsUpdate(v); + value_->setDeviceNeedsUpdate(v); +} + +bool SparseMatrixStorage::deviceAllocated() const { + return outer_->deviceAllocated() && + inner_->deviceAllocated() && + value_->deviceAllocated(); +} + +void SparseMatrixStorage::allocateDevice() const { + outer_->allocateDevice(); + inner_->allocateDevice(); + value_->allocateDevice(); +} + +void SparseMatrixStorage::deallocateDevice() const { + outer_->deallocateDevice(); + inner_->deallocateDevice(); + value_->deallocateDevice(); +} + + + + +} \ No newline at end of file diff --git a/src/atlas/linalg/sparse/SparseMatrixStorage.h b/src/atlas/linalg/sparse/SparseMatrixStorage.h new file mode 100644 index 000000000..624c09664 --- /dev/null +++ b/src/atlas/linalg/sparse/SparseMatrixStorage.h @@ -0,0 +1,319 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include +#include +#include +#include + +#include "eckit/linalg/SparseMatrix.h" + +#include "atlas/array.h" +#include "atlas/runtime/Exception.h" +#include "atlas/linalg/sparse/SparseMatrixView.h" + +namespace atlas { +namespace linalg { + +//---------------------------------------------------------------------------------------------------------------------- + +/// @brief SparseMatrixStorage +/// Storage class based on atlas::array::Array with GPU offload capability +/// This class contains only the data and no operators, iterators, or special +/// sparse matrix construction from triplets etc. +/// +/// The construction is handled via SparseMatrixConvertor class. +/// Examples: +/// +/// To construct it, taking ownership of a constructed Eigen::SparseMatrix +/// +/// SparseMatrixStorage s = make_sparse_matrix_storage(std::move(eigen_matrix)); +/// +/// To construct it, taking a copy of a constructed Eigen::SparseMatrix +/// +/// SparseMatrixStorage s = make_sparse_matrix_storage(eigen_matrix); +/// +/// To construct it, taking ownership of a constructed eckit::linalg::SparseMatrix, avoiding copies if data types match +/// +/// SparseMatrixStorage s = make_sparse_matrix_storage(std::move(eckit_matrix)); +/// +/// +/// To construct it, taking a copy of a constructed eckit::linalg::SparseMatrix (no std::move) +/// +/// SparseMatrixStorage s = make_sparse_matrix_storage(eckit_matrix); +/// +/// +/// It is also possible to initialise empty and move into it at later stage: +/// +/// SparseMatrixStorage s; +/// s = make_sparse_matrix_storage(eckit_matrix); +/// +/// +/// To construct it, taking a single precision copy of a constructed eckit::linalg::SparseMatrix +/// +/// s = make_sparse_matrix_storage(eckit_matrix); +/// + +class SparseMatrixStorage { +public: + + // Virtual destructor + virtual ~SparseMatrixStorage() = default; + + /// Default constructor + SparseMatrixStorage() = default; + + /// Move constructor, takes ownership! + SparseMatrixStorage(SparseMatrixStorage&& other); + + /// Copy constructor, makes copy! + SparseMatrixStorage(const SparseMatrixStorage& other); + + /// Copy from a SparseMatrixView + template + SparseMatrixStorage(const SparseMatrixView& host_view); + + /// Move assign from other SparseMatrixStorage, takes ownership! + SparseMatrixStorage& operator=(SparseMatrixStorage&& other); + + /// Copy assign from other SparseMatrixStorage, takes ownership! + SparseMatrixStorage& operator=(const SparseMatrixStorage& other); + + /// Empty if rows and cols are zero. + bool empty() const { return rows_ == 0 && cols_ == 0; } + + /// Footprint in bytes in host memory space + std::size_t footprint() const { return value_->footprint() + outer_->footprint() + inner_->footprint(); } + + const atlas::array::Array& value() const { return *value_; } + const atlas::array::Array& outer() const { return *outer_; } + const atlas::array::Array& inner() const { return *inner_; } + std::size_t rows() const { return rows_;} + std::size_t cols() const { return cols_;} + std::size_t nnz() const { return nnz_;} + + void updateDevice() const; + + void updateHost() const; + + bool hostNeedsUpdate() const; + + bool deviceNeedsUpdate() const; + + void setHostNeedsUpdate(bool v) const; + + void setDeviceNeedsUpdate(bool v) const; + + bool deviceAllocated() const; + + void allocateDevice() const; + + void deallocateDevice() const; + + static SparseMatrixStorage make( + std::size_t rows, + std::size_t cols, + std::size_t nnz, + std::unique_ptr&& value, + std::unique_ptr&& inner, + std::unique_ptr&& outer, + std::any&& storage) { + SparseMatrixStorage S; + S.rows_ = rows; + S.cols_ = cols; + S.nnz_ = nnz; + S.outer_ = std::move(outer); + S.inner_ = std::move(inner); + S.value_ = std::move(value); + S.storage_ = std::move(storage); + return S; + } + + void swap(SparseMatrixStorage& other) { + std::swap(other.rows_, rows_); + std::swap(other.cols_, cols_); + std::swap(other.nnz_, nnz_); + std::swap(other.value_, value_); + std::swap(other.inner_, inner_); + std::swap(other.outer_, outer_); + std::swap(other.storage_, storage_); + } + + bool contains(DataType value_type, DataType index_type) { + if (empty()) { + return false; + } + return value_->datatype() == value_type && outer_->datatype() == index_type; + } + + template + bool contains() { + return contains(make_datatype(), make_datatype()); + } + +protected: + std::size_t nnz_{0}; + std::size_t rows_{0}; + std::size_t cols_{0}; + std::unique_ptr outer_; + std::unique_ptr inner_; + std::unique_ptr value_; + + std::any storage_; + // This storage is usually empty. + // It is used to allow to move alternative sparse matrix formats into it, + // and then wrap this data using the atlas Arrays + +private: + template + void host_copy(const ValueT* input_data, atlas::array::Array& output) { + auto size = output.size(); + ValueT* output_data = output.host_data(); + std::copy( input_data, input_data + size, output_data ); + } +}; + +//---------------------------------------------------------------------------------------------------------------------- + +template +SparseMatrixStorage::SparseMatrixStorage(const SparseMatrixView& host_view) { + nnz_ = host_view.nnz_; + rows_ = host_view.rows_; + cols_ = host_view.cols_; + outer_.reset(atlas::array::Array::create(host_view.outer_size())); + inner_.reset(atlas::array::Array::create(host_view.inner_size())); + value_.reset(atlas::array::Array::create(host_view.value_size())); + host_copy(host_view.outer(),*outer_); + host_copy(host_view.inner(),*inner_); + host_copy(host_view.value(),*value_); +} + +//---------------------------------------------------------------------------------------------------------------------- + +namespace { +template +void host_copy(const InputT* input_data, array::Array& output) { + auto size = output.size(); + OutputT* output_data = output.host_data(); + std::copy( input_data, input_data + size, output_data ); +} + +template +void host_copy(const array::Array& input, array::Array& output) { + host_copy( input.host_data(), output ); +} + +template +void host_copy(const array::Array& input, array::Array& output) { + switch(input.datatype().kind()) { + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + default: ATLAS_NOTIMPLEMENTED; + } +} +void host_copy(const array::Array& input, array::Array& output) { + switch(output.datatype().kind()) { + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + case DataType::kind(): return host_copy( input, output ); + default: ATLAS_NOTIMPLEMENTED; + } +} + +} + +template +SparseMatrixStorage make_sparse_matrix_storage(const SparseMatrixStorage& other) { + auto rows = other.rows(); + auto cols = other.cols(); + auto nnz = other.nnz(); + std::unique_ptr value(array::Array::create(nnz)); + std::unique_ptr inner(array::Array::create(nnz)); + std::unique_ptr outer(array::Array::create(rows+1)); + host_copy(other.value(), *value); + host_copy(other.inner(), *inner); + host_copy(other.outer(), *outer); + return SparseMatrixStorage::make(rows,cols,nnz,std::move(value), std::move(inner), std::move(outer), std::any()); +} + +template +SparseMatrixStorage make_sparse_matrix_storage(SparseMatrixStorage&& other) { + SparseMatrixStorage S; + + if (other.contains()) { + S = std::move(other); + } + else { + auto rows = other.rows(); + auto cols = other.cols(); + auto nnz = other.nnz(); + std::unique_ptr value(array::Array::create(nnz)); + std::unique_ptr inner(array::Array::create(nnz)); + std::unique_ptr outer(array::Array::create(rows+1)); + host_copy(other.value(), *value); + host_copy(other.inner(), *inner); + host_copy(other.outer(), *outer); + S = SparseMatrixStorage::make(rows,cols,nnz,std::move(value), std::move(inner), std::move(outer), std::any()); + } + return S; +} + +//---------------------------------------------------------------------------------------------------------------------- + +template +inline SparseMatrixView make_host_view(const SparseMatrixStorage& m) { + if( m.value().datatype().kind() != DataType::kind() || m.outer().datatype().kind() != DataType::kind() ) { + ATLAS_THROW_EXCEPTION("Cannot make_host_view<" + DataType::str() + "," << DataType::str() + + ">(const SparseMatrixStorage&) from SparseMatrixStorage containing values of type <" + m.value().datatype().str() + "> and indices of type <" + m.outer().datatype().str() +">" ); + } + return SparseMatrixView { + m.rows(), + m.cols(), + m.nnz(), + m.value().host_data(), + m.inner().host_data(), + m.outer().host_data() + }; +} + +//---------------------------------------------------------------------------------------------------------------------- + +template +inline SparseMatrixView make_device_view(const SparseMatrixStorage& m) { + if( m.value().datatype().kind() != DataType::kind() || m.outer().datatype().kind() != DataType::kind() ) { + ATLAS_THROW_EXCEPTION("Cannot make_device_view<" + DataType::str() + "," << DataType::str() + + ">(const SparseMatrixStorage&) from SparseMatrixStorage containing values of type <" + m.value().datatype().str() + "> and indices of type <" + m.outer().datatype().str() +">" ); + } + if( ! m.deviceAllocated() ) { + ATLAS_THROW_EXCEPTION("Cannot make_device_view(const SparseMatrixStorage&) as the device data is not allocated"); + } + return SparseMatrixView{ + m.rows(), + m.cols(), + m.nnz(), + m.value().device_data(), + m.inner().device_data(), + m.outer().device_data() + }; +} + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace linalg +} // namespace atlas \ No newline at end of file diff --git a/src/atlas/linalg/sparse/SparseMatrixView.h b/src/atlas/linalg/sparse/SparseMatrixView.h new file mode 100644 index 000000000..875a4164f --- /dev/null +++ b/src/atlas/linalg/sparse/SparseMatrixView.h @@ -0,0 +1,64 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#pragma once + +#include "eckit/linalg/types.h" + +namespace atlas { +namespace linalg { + +//---------------------------------------------------------------------------------------------------------------------- + +template +class SparseMatrixView { +public: + using value_type = Value; + using index_type = Index; + using size_type = std::size_t; + + ~SparseMatrixView() = default; + + SparseMatrixView() = default; + + SparseMatrixView( + size_type rows, + size_type cols, + size_type nnz, + value_type const* value, + index_type const* inner, + index_type const* outer) : + rows_(rows), cols_(cols), nnz_(nnz), outer_(outer), inner_(inner), value_(value) { + } + SparseMatrixView(const SparseMatrixView& other) = default; + + size_type rows() const { return rows_; } + size_type cols() const { return cols_; } + size_type nnz() const { return nnz_; } + size_type outer_size() const { return rows() + 1; } + size_type inner_size() const { return nnz(); } + size_type value_size() const { return nnz(); } + index_type const* outer() const { return outer_; } + index_type const* inner() const { return inner_; } + value_type const* value() const { return value_; } + bool empty() const { return nnz() == 0; } + +private: + size_type rows_{0}; + size_type cols_{0}; + size_type nnz_{0}; + index_type const* outer_{nullptr}; + index_type const* inner_{nullptr}; + value_type const* value_{nullptr}; +}; + +//---------------------------------------------------------------------------------------------------------------------- +} // namespace linalg +} // namespace atlas diff --git a/src/tests/interpolation/test_interpolation_finite_element_cached.cc b/src/tests/interpolation/test_interpolation_finite_element_cached.cc index 99d9803c6..605e44021 100644 --- a/src/tests/interpolation/test_interpolation_finite_element_cached.cc +++ b/src/tests/interpolation/test_interpolation_finite_element_cached.cc @@ -106,13 +106,18 @@ CASE("extract cache, copy it, and move it for use") { set_field(field_source, grid_source, func); - eckit::linalg::SparseMatrix matrix = get_or_create_cache(grid_source, grid_target).matrix(); + interpolation::MatrixCache created_cache = get_or_create_cache(grid_source, grid_target); + eckit::linalg::SparseMatrix eckit_matrix_view = atlas::linalg::make_non_owning_eckit_sparse_matrix(created_cache.matrix()); + eckit::linalg::SparseMatrix eckit_matrix_copy = atlas::linalg::make_eckit_sparse_matrix(created_cache.matrix()); - EXPECT(not matrix.empty()); + EXPECT(not created_cache.matrix().empty()); + EXPECT(not eckit_matrix_view.empty()); + EXPECT(not eckit_matrix_copy.empty()); - auto cache = interpolation::MatrixCache(std::move(matrix)); + auto cache = interpolation::MatrixCache(atlas::linalg::make_sparse_matrix_storage(std::move(eckit_matrix_copy))); - EXPECT(matrix.empty()); + EXPECT(not eckit_matrix_view.empty()); // We didn't touch the view + EXPECT(eckit_matrix_copy.empty()); // This has been moved into the new 'cache' ATLAS_TRACE_SCOPE("Interpolate with cache") { Interpolation interpolation_using_cache(option::type("finite-element"), grid_source, grid_target, cache); @@ -133,7 +138,9 @@ CASE("extract cache, copy it, and pass non-owning pointer") { set_field(field_source, grid_source, func); - eckit::linalg::SparseMatrix matrix = get_or_create_cache(grid_source, grid_target).matrix(); + const auto& matrix_storage = get_or_create_cache(grid_source, grid_target).matrix(); + atlas::linalg::SparseMatrixStorage matrix_storage_copy(matrix_storage); + auto matrix = atlas::linalg::make_host_view(matrix_storage_copy); EXPECT(not matrix.empty()); @@ -155,7 +162,7 @@ CASE("extract cache, copy it, and pass non-owning pointer") { check_field(field_target, grid_target, func, 1.e-4); set_field(field_target, 0.); - auto cache = interpolation::MatrixCache(&matrix); + auto cache = interpolation::MatrixCache( &matrix_storage_copy ); EXPECT(not matrix.empty()); diff --git a/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc b/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc index 641005fb6..83e1aacb3 100644 --- a/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc +++ b/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc @@ -40,10 +40,10 @@ class Access { std::string hash() { eckit::MD5 hash; - const auto& m = matrix(); + const auto m = atlas::linalg::make_host_view(matrix()); const auto outer = m.outer(); const auto index = m.inner(); - const auto weight = m.data(); + const auto weight = m.value(); idx_t rows = static_cast(m.rows()); hash.add(rows); diff --git a/src/tests/interpolation/test_interpolation_non_linear.cc b/src/tests/interpolation/test_interpolation_non_linear.cc index 61a879ebc..6ebb54a25 100644 --- a/src/tests/interpolation/test_interpolation_non_linear.cc +++ b/src/tests/interpolation/test_interpolation_non_linear.cc @@ -91,7 +91,6 @@ CASE("Interpolation with MissingValue") { } } - SECTION("missing-if-any-missing") { Interpolation interpolation(Config("type", "finite-element").set("non_linear", "missing-if-any-missing"), fsA, fsB); diff --git a/src/tests/linalg/CMakeLists.txt b/src/tests/linalg/CMakeLists.txt index d060dcc8a..ec70d891f 100644 --- a/src/tests/linalg/CMakeLists.txt +++ b/src/tests/linalg/CMakeLists.txt @@ -14,6 +14,23 @@ ecbuild_add_test( TARGET atlas_test_linalg_sparse ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) +ecbuild_add_test( TARGET atlas_test_linalg_sparse_matrix + SOURCES test_linalg_sparse_matrix.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + +ecbuild_add_test( TARGET atlas_test_linalg_sparse_matrix_gpu + SOURCES test_linalg_sparse_matrix_gpu.cc + LIBS atlas hicsparse + CONDITION HAVE_CUDA OR HAVE_HIP + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) +if( TARGET atlas_test_linalg_sparse_matrix_gpu ) + set_tests_properties( atlas_test_linalg_sparse_matrix_gpu PROPERTIES LABELS "gpu") +endif() + + ecbuild_add_test( TARGET atlas_test_linalg_dense SOURCES test_linalg_dense.cc LIBS atlas diff --git a/src/tests/linalg/test_linalg_sparse.cc b/src/tests/linalg/test_linalg_sparse.cc index 0e6f8c81f..0efe5d387 100644 --- a/src/tests/linalg/test_linalg_sparse.cc +++ b/src/tests/linalg/test_linalg_sparse.cc @@ -1,5 +1,5 @@ /* - * (C) Copyright 1996- ECMWF. + * (C) Copyright 2024- ECMWF. * * This software is licensed under the terms of the Apache Licence Version 2.0 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. @@ -25,6 +25,8 @@ using namespace atlas::linalg; namespace atlas { namespace test { +using SparseMatrix = eckit::linalg::SparseMatrix; + //---------------------------------------------------------------------------------------------------------------------- // strings to be used in the tests diff --git a/src/tests/linalg/test_linalg_sparse_matrix.cc b/src/tests/linalg/test_linalg_sparse_matrix.cc new file mode 100644 index 000000000..88a87b0c0 --- /dev/null +++ b/src/tests/linalg/test_linalg_sparse_matrix.cc @@ -0,0 +1,239 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "tests/AtlasTestEnvironment.h" + + +#include "atlas/linalg/sparse/SparseMatrixStorage.h" +#include "atlas/linalg/sparse/SparseMatrixView.h" +#include "atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h" +#include "atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h" + +using atlas::linalg::SparseMatrixStorage; +using atlas::linalg::make_sparse_matrix_storage; + +namespace atlas { +namespace test { + +template +void do_matrix_multiply(const atlas::linalg::SparseMatrixView& W, const XValue* x, YValue* y) { + const auto outer = W.outer(); + const auto inner = W.inner(); + const auto weight = W.value(); + int rows = W.rows(); + for (int r = 0; r < rows; ++r) { + y[r] = 0.; + for (IndexType c = outer[r]; c < outer[r + 1]; ++c) { + IndexType n = inner[c]; + y[r] += weight[c] * x[n]; + } + } +} + +#if ATLAS_HAVE_EIGEN +template +Eigen::SparseMatrix create_eigen_sparsematrix() { + Eigen::SparseMatrix matrix(3,3); + std::vector> triplets; + triplets.reserve(4); + triplets.emplace_back(0, 0, 2.); + triplets.emplace_back(0, 2, -3.); + triplets.emplace_back(1, 1, 2.); + triplets.emplace_back(2, 2, 2.); + matrix.setFromTriplets(triplets.begin(), triplets.end()); + return matrix; +} +#endif + +eckit::linalg::SparseMatrix create_eckit_sparsematrix() { + return eckit::linalg::SparseMatrix{3, 3, {{0, 0, 2.}, {0, 2, -3.}, {1, 1, 2.}, {2, 2, 2.}}}; +} + +template +SparseMatrixStorage create_sparsematrix_storage() { + return make_sparse_matrix_storage(create_eckit_sparsematrix()); +} + +template +void check_array(const array::Array& array, const std::vector& expected) { + EXPECT_EQ(array.size(), expected.size()); + auto view = array::make_view(array); + for( std::size_t i=0; i +void check_matrix(const SparseMatrixStorage& m) { + std::vector expected_value{2., -3., 2., 2.}; + std::vector expected_outer{0, 2, 3, 4}; + std::vector expected_inner{0, 2, 1, 2}; + + EXPECT_EQ(m.rows(), 3); + EXPECT_EQ(m.cols(), 3); + EXPECT_EQ(m.nnz(), 4); + + check_array(m.value(),expected_value); + check_array(m.outer(),expected_outer); + check_array(m.inner(),expected_inner); + + auto host_matrix_view = linalg::make_host_view(m); + + std::vector host_outer(host_matrix_view.outer_size()); + std::vector host_inner(host_matrix_view.inner_size()); + std::vector host_value(host_matrix_view.value_size()); + + array::ArrayT x(m.cols()); + array::ArrayT y(m.rows()); + array::make_host_view(x).assign({1.0, 2.0, 3.0}); + + do_matrix_multiply(host_matrix_view, x.host_data(), y.host_data()); + + // Check result + std::vector expected_y{-7.0, 4.0, 6.0}; + check_array(y, expected_y); +} + +void check_matrix(const SparseMatrixStorage& m) { + if (m.value().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + +CASE("Create SparseMatrix moving eckit::linalg::SparseMatrix") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eckit_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create SparseMatrix moving eckit::linalg::SparseMatrix using std::move") { + auto A = create_eckit_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(std::move(A)); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create SparseMatrix copying eckit::linalg::SparseMatrix") { + auto A = create_eckit_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(A); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create single precision SparseMatrix copying from double precision SparseMatrix ") { + auto Sdouble = create_sparsematrix_storage(); + SparseMatrixStorage S = make_sparse_matrix_storage(Sdouble); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create single precision SparseMatrix moving from double precision SparseMatrix which causes copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create double precision SparseMatrix copying from single precision SparseMatrix ") { + auto Sfloat = create_sparsematrix_storage(); + SparseMatrixStorage S = make_sparse_matrix_storage(Sfloat); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create double precision SparseMatrix moving from single precision SparseMatrix which causes copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with copy constructor from double precision") { + auto Sdouble = create_sparsematrix_storage(); + SparseMatrixStorage S(Sdouble); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with move constructor from double precision") { + SparseMatrixStorage S(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create base with copy constructor from single precision") { + auto Sfloat = create_sparsematrix_storage(); + SparseMatrixStorage S(Sfloat); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with move constructor from single precision") { + SparseMatrixStorage S(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} + +#if ATLAS_HAVE_EIGEN +CASE("Copy from Eigen double") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move from Eigen double, avoiding copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy and convert double from Eigen to single precision") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move and convert double from Eigen to single precision, should cause copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy from Eigen single") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move from Eigen single, avoiding copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy and convert single from Eigen to double precision") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move and convert single from Eigen to double precision, should cause copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create chain of moves from eigen double") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create chain of moves from eigen single") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} +#endif + +CASE("Test exceptions in make_host_view") { + SparseMatrixStorage S( create_sparsematrix_storage() ); + EXPECT_THROWS(linalg::make_host_view(S)); // value data type mismatch + EXPECT_THROWS((linalg::make_host_view(S))); // index data type mismatch +} + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} diff --git a/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc b/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc new file mode 100644 index 000000000..eaf003f31 --- /dev/null +++ b/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc @@ -0,0 +1,353 @@ +/* + * (C) Copyright 2024- ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +// #include "atlas/linalg/SparseMatrix.h" + +#include "hic/hic.h" +#include "hic/hicsparse.h" + +#include "tests/AtlasTestEnvironment.h" + + +#include "atlas/linalg/sparse/SparseMatrixStorage.h" +#include "atlas/linalg/sparse/SparseMatrixView.h" +#include "atlas/linalg/sparse/MakeSparseMatrixStorageEigen.h" +#include "atlas/linalg/sparse/MakeSparseMatrixStorageEckit.h" + +using atlas::linalg::SparseMatrixStorage; +using atlas::linalg::make_sparse_matrix_storage; +using atlas::linalg::make_host_view; +using atlas::linalg::make_device_view; + +namespace atlas { +namespace test { + +template +void do_hicsparse_matrix_multiply(const atlas::linalg::SparseMatrixView& device_A, const XValue* device_x, YValue* device_y) { + // Create sparse library handle + hicsparseHandle_t handle; + HICSPARSE_CALL( hicsparseCreate(&handle) ); + + auto get_hic_value_type = [](const auto& dummy) -> hicDataType { + if(std::is_same_v,double> ) { + return HIC_R_64F; + } + else if (std::is_same_v,float> ) { + return HIC_R_32F; + } + else if (std::is_same_v,int> ) { + return HIC_R_32F; + } + ATLAS_NOTIMPLEMENTED; + }; + auto get_hic_index_type = [](const auto& dummy) -> hicsparseIndexType_t { + if (std::is_same_v,int> ) { + return HICSPARSE_INDEX_32I; + } + else if (std::is_same_v,long> ) { + return HICSPARSE_INDEX_64I; + } + ATLAS_NOTIMPLEMENTED; + }; + auto compute_type = get_hic_value_type(MatrixValue()); + auto index_type = get_hic_index_type(IndexType()); + + hicsparseConstSpMatDescr_t matA; + HICSPARSE_CALL( hicsparseCreateConstCsr( + &matA, + device_A.rows(), device_A.cols(), device_A.nnz(), + device_A.outer(), // row_offsets + device_A.inner(), // column_indices + device_A.value(), // values + index_type, + index_type, + HICSPARSE_INDEX_BASE_ZERO, + get_hic_value_type(MatrixValue())) ); + + // Create dense matrix descriptors + hicsparseConstDnVecDescr_t vecX; + HICSPARSE_CALL( hicsparseCreateConstDnVec( + &vecX, + device_A.cols(), + device_x, + get_hic_value_type(XValue())) ); + + hicsparseDnVecDescr_t vecY; + HICSPARSE_CALL( hicsparseCreateDnVec( + &vecY, + device_A.rows(), + device_y, + get_hic_value_type(YValue())) ); + + // Set parameters in compute_type + const MatrixValue alpha = 1.0; + const MatrixValue beta = 0.0; + + // Determine buffer size + size_t bufferSize = 0; + HICSPARSE_CALL( hicsparseSpMV_bufferSize( + handle, + HICSPARSE_OPERATION_NON_TRANSPOSE, + &alpha, + matA, + vecX, + &beta, + vecY, + compute_type, + HICSPARSE_SPMV_ALG_DEFAULT, + &bufferSize) ); + + // Allocate buffer + char* buffer; + HIC_CALL( hicMalloc(&buffer, bufferSize) ); + + // Perform SpMV + HICSPARSE_CALL( hicsparseSpMV( + handle, + HICSPARSE_OPERATION_NON_TRANSPOSE, + &alpha, + matA, + vecX, + &beta, + vecY, + compute_type, + HICSPARSE_SPMV_ALG_DEFAULT, + buffer) ); + + HIC_CALL( hicFree(buffer) ); + HICSPARSE_CALL( hicsparseDestroyDnVec(vecY) ); + HICSPARSE_CALL( hicsparseDestroyDnVec(vecX) ); + HICSPARSE_CALL( hicsparseDestroySpMat(matA) ); + HICSPARSE_CALL( hicsparseDestroy(handle) ); +} + +#if ATLAS_HAVE_EIGEN +template +Eigen::SparseMatrix create_eigen_sparsematrix() { + Eigen::SparseMatrix matrix(3,3); + std::vector> triplets; + triplets.reserve(4); + triplets.emplace_back(0, 0, 2.); + triplets.emplace_back(0, 2, -3.); + triplets.emplace_back(1, 1, 2.); + triplets.emplace_back(2, 2, 2.); + matrix.setFromTriplets(triplets.begin(), triplets.end()); + return matrix; +} +#endif + +eckit::linalg::SparseMatrix create_eckit_sparsematrix() { + return eckit::linalg::SparseMatrix{3, 3, {{0, 0, 2.}, {0, 2, -3.}, {1, 1, 2.}, {2, 2, 2.}}}; +} + +template +SparseMatrixStorage create_sparsematrix_storage() { + return make_sparse_matrix_storage(create_eckit_sparsematrix()); +} + +template +void check_array(const array::Array& array, const std::vector& expected) { + EXPECT_EQ(array.size(), expected.size()); + auto view = array::make_view(array); + for( std::size_t i=0; i +void check_matrix(const SparseMatrixStorage& m) { + std::vector expected_value{2., -3., 2., 2.}; + std::vector expected_outer{0, 2, 3, 4}; + std::vector expected_inner{0, 2, 1, 2}; + + EXPECT_EQ(m.rows(), 3); + EXPECT_EQ(m.cols(), 3); + EXPECT_EQ(m.nnz(), 4); + + check_array(m.value(),expected_value); + check_array(m.outer(),expected_outer); + check_array(m.inner(),expected_inner); + + m.updateDevice(); + auto host_matrix_view = make_host_view(m); + auto device_matrix_view = make_device_view(m); + + std::vector host_outer(host_matrix_view.outer_size()); + std::vector host_inner(host_matrix_view.inner_size()); + std::vector host_value(host_matrix_view.value_size()); + + hicMemcpy(host_outer.data(), device_matrix_view.outer(), device_matrix_view.outer_size() * sizeof(eckit::linalg::Index), hicMemcpyDeviceToHost); + hicMemcpy(host_inner.data(), device_matrix_view.inner(), device_matrix_view.inner_size() * sizeof(eckit::linalg::Index), hicMemcpyDeviceToHost); + hicMemcpy(host_value.data(), device_matrix_view.value(), device_matrix_view.value_size() * sizeof(Value), hicMemcpyDeviceToHost); + + EXPECT(host_outer == expected_outer); + EXPECT(host_inner == expected_inner); + EXPECT(host_value == expected_value); + + array::ArrayT x(m.cols()); + array::ArrayT y(m.rows()); + array::make_host_view(x).assign({1.0, 2.0, 3.0}); + + m.updateDevice(); + x.updateDevice(); + y.allocateDevice(); + + do_hicsparse_matrix_multiply(device_matrix_view, x.device_data(), y.device_data()); + + // Check result + y.updateHost(); + std::vector expected_y{-7.0, 4.0, 6.0}; + check_array(y, expected_y); +} + +void check_matrix(const SparseMatrixStorage& m) { + if (m.value().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else { + ATLAS_NOTIMPLEMENTED; + } +} + +CASE("Create SparseMatrix moving eckit::linalg::SparseMatrix") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eckit_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create SparseMatrix moving eckit::linalg::SparseMatrix using std::move") { + auto A = create_eckit_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(std::move(A)); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create SparseMatrix copying eckit::linalg::SparseMatrix") { + auto A = create_eckit_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(A); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create single precision SparseMatrix copying from double precision SparseMatrix ") { + auto Sdouble = create_sparsematrix_storage(); + SparseMatrixStorage S = make_sparse_matrix_storage(Sdouble); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create single precision SparseMatrix moving from double precision SparseMatrix which causes copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create double precision SparseMatrix copying from single precision SparseMatrix ") { + auto Sfloat = create_sparsematrix_storage(); + SparseMatrixStorage S = make_sparse_matrix_storage(Sfloat); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create double precision SparseMatrix moving from single precision SparseMatrix which causes copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with copy constructor from double precision") { + auto Sdouble = create_sparsematrix_storage(); + SparseMatrixStorage S(Sdouble); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with move constructor from double precision") { + SparseMatrixStorage S(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create base with copy constructor from single precision") { + auto Sfloat = create_sparsematrix_storage(); + SparseMatrixStorage S(Sfloat); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create base with move constructor from single precision") { + SparseMatrixStorage S(create_sparsematrix_storage()); + EXPECT_NO_THROW(check_matrix(S)); +} + +#if ATLAS_HAVE_EIGEN +CASE("Copy from Eigen double") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move from Eigen double, avoiding copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy and convert double from Eigen to single precision") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move and convert double from Eigen to single precision, should cause copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy from Eigen single") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move from Eigen single, avoiding copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Copy and convert single from Eigen to double precision") { + auto eigen_matrix = create_eigen_sparsematrix(); + SparseMatrixStorage S = make_sparse_matrix_storage(eigen_matrix); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Move and convert single from Eigen to double precision, should cause copy") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +CASE("Create chain of moves from eigen double") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} +CASE("Create chain of moves from eigen single") { + SparseMatrixStorage S = make_sparse_matrix_storage(create_eigen_sparsematrix()); + EXPECT_NO_THROW(check_matrix(S)); +} + +#endif + +CASE("Test exceptions in make_host_view") { + SparseMatrixStorage S( create_sparsematrix_storage() ); + EXPECT_THROWS(make_host_view(S)); // value data type mismatch + EXPECT_THROWS((make_host_view(S))); // index data type mismatch +} +CASE("Test exceptions in make_device_view") { + SparseMatrixStorage S( create_sparsematrix_storage() ); + EXPECT_THROWS(make_device_view(S)); // device is not allocated + S.updateDevice(); + EXPECT_THROWS(make_device_view(S)); // value data type mismatch + EXPECT_THROWS((make_device_view(S))); // index data type mismatch +} + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} From bd9ab9aa7e5bd666d0122fe9648d9dbb3efd511a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Thu, 5 Dec 2024 16:08:26 +0000 Subject: [PATCH 2/2] Adapt tests for the hypothetical case when eckit::linalg::Index would become an unsigned integer --- ...est_interpolation_finite_element_cached.cc | 2 +- ...test_interpolation_k_nearest_neighbours.cc | 2 +- src/tests/linalg/test_linalg_sparse_matrix.cc | 44 ++++++++++----- .../linalg/test_linalg_sparse_matrix_gpu.cc | 56 ++++++++++++------- 4 files changed, 70 insertions(+), 34 deletions(-) diff --git a/src/tests/interpolation/test_interpolation_finite_element_cached.cc b/src/tests/interpolation/test_interpolation_finite_element_cached.cc index 605e44021..5d75a8dfd 100644 --- a/src/tests/interpolation/test_interpolation_finite_element_cached.cc +++ b/src/tests/interpolation/test_interpolation_finite_element_cached.cc @@ -140,7 +140,7 @@ CASE("extract cache, copy it, and pass non-owning pointer") { const auto& matrix_storage = get_or_create_cache(grid_source, grid_target).matrix(); atlas::linalg::SparseMatrixStorage matrix_storage_copy(matrix_storage); - auto matrix = atlas::linalg::make_host_view(matrix_storage_copy); + auto matrix = atlas::linalg::make_host_view(matrix_storage_copy); EXPECT(not matrix.empty()); diff --git a/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc b/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc index 83e1aacb3..86677c1ae 100644 --- a/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc +++ b/src/tests/interpolation/test_interpolation_k_nearest_neighbours.cc @@ -40,7 +40,7 @@ class Access { std::string hash() { eckit::MD5 hash; - const auto m = atlas::linalg::make_host_view(matrix()); + const auto m = atlas::linalg::make_host_view(matrix()); const auto outer = m.outer(); const auto index = m.inner(); const auto weight = m.value(); diff --git a/src/tests/linalg/test_linalg_sparse_matrix.cc b/src/tests/linalg/test_linalg_sparse_matrix.cc index 88a87b0c0..f6c43645b 100644 --- a/src/tests/linalg/test_linalg_sparse_matrix.cc +++ b/src/tests/linalg/test_linalg_sparse_matrix.cc @@ -38,7 +38,7 @@ void do_matrix_multiply(const atlas::linalg::SparseMatrixView +template > Eigen::SparseMatrix create_eigen_sparsematrix() { Eigen::SparseMatrix matrix(3,3); std::vector> triplets; @@ -71,11 +71,11 @@ void check_array(const array::Array& array, const std::vector& expected) } -template +template void check_matrix(const SparseMatrixStorage& m) { - std::vector expected_value{2., -3., 2., 2.}; - std::vector expected_outer{0, 2, 3, 4}; - std::vector expected_inner{0, 2, 1, 2}; + std::vector expected_value{2., -3., 2., 2.}; + std::vector expected_outer{0, 2, 3, 4}; + std::vector expected_inner{0, 2, 1, 2}; EXPECT_EQ(m.rows(), 3); EXPECT_EQ(m.cols(), 3); @@ -85,11 +85,11 @@ void check_matrix(const SparseMatrixStorage& m) { check_array(m.outer(),expected_outer); check_array(m.inner(),expected_inner); - auto host_matrix_view = linalg::make_host_view(m); + auto host_matrix_view = linalg::make_host_view(m); - std::vector host_outer(host_matrix_view.outer_size()); - std::vector host_inner(host_matrix_view.inner_size()); - std::vector host_value(host_matrix_view.value_size()); + std::vector host_outer(host_matrix_view.outer_size()); + std::vector host_inner(host_matrix_view.inner_size()); + std::vector host_value(host_matrix_view.value_size()); array::ArrayT x(m.cols()); array::ArrayT y(m.rows()); @@ -103,11 +103,29 @@ void check_matrix(const SparseMatrixStorage& m) { } void check_matrix(const SparseMatrixStorage& m) { - if (m.value().datatype() == atlas::make_datatype()) { - check_matrix(m); + if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); } - else if (m.value().datatype() == atlas::make_datatype()) { - check_matrix(m); + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); } else { ATLAS_NOTIMPLEMENTED; diff --git a/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc b/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc index eaf003f31..5bb2da54d 100644 --- a/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc +++ b/src/tests/linalg/test_linalg_sparse_matrix_gpu.cc @@ -48,10 +48,10 @@ void do_hicsparse_matrix_multiply(const atlas::linalg::SparseMatrixView hicsparseIndexType_t { - if (std::is_same_v,int> ) { + if (std::is_same_v>,int> ) { return HICSPARSE_INDEX_32I; } - else if (std::is_same_v,long> ) { + else if (std::is_same_v>,long> ) { return HICSPARSE_INDEX_64I; } ATLAS_NOTIMPLEMENTED; @@ -129,7 +129,7 @@ void do_hicsparse_matrix_multiply(const atlas::linalg::SparseMatrixView +template > Eigen::SparseMatrix create_eigen_sparsematrix() { Eigen::SparseMatrix matrix(3,3); std::vector> triplets; @@ -162,11 +162,11 @@ void check_array(const array::Array& array, const std::vector& expected) } -template +template void check_matrix(const SparseMatrixStorage& m) { - std::vector expected_value{2., -3., 2., 2.}; - std::vector expected_outer{0, 2, 3, 4}; - std::vector expected_inner{0, 2, 1, 2}; + std::vector expected_value{2., -3., 2., 2.}; + std::vector expected_outer{0, 2, 3, 4}; + std::vector expected_inner{0, 2, 1, 2}; EXPECT_EQ(m.rows(), 3); EXPECT_EQ(m.cols(), 3); @@ -177,16 +177,16 @@ void check_matrix(const SparseMatrixStorage& m) { check_array(m.inner(),expected_inner); m.updateDevice(); - auto host_matrix_view = make_host_view(m); - auto device_matrix_view = make_device_view(m); + auto host_matrix_view = make_host_view(m); + auto device_matrix_view = make_device_view(m); - std::vector host_outer(host_matrix_view.outer_size()); - std::vector host_inner(host_matrix_view.inner_size()); - std::vector host_value(host_matrix_view.value_size()); + std::vector host_outer(host_matrix_view.outer_size()); + std::vector host_inner(host_matrix_view.inner_size()); + std::vector host_value(host_matrix_view.value_size()); - hicMemcpy(host_outer.data(), device_matrix_view.outer(), device_matrix_view.outer_size() * sizeof(eckit::linalg::Index), hicMemcpyDeviceToHost); - hicMemcpy(host_inner.data(), device_matrix_view.inner(), device_matrix_view.inner_size() * sizeof(eckit::linalg::Index), hicMemcpyDeviceToHost); - hicMemcpy(host_value.data(), device_matrix_view.value(), device_matrix_view.value_size() * sizeof(Value), hicMemcpyDeviceToHost); + hicMemcpy(host_outer.data(), device_matrix_view.outer(), device_matrix_view.outer_size() * sizeof(Index), hicMemcpyDeviceToHost); + hicMemcpy(host_inner.data(), device_matrix_view.inner(), device_matrix_view.inner_size() * sizeof(Index), hicMemcpyDeviceToHost); + hicMemcpy(host_value.data(), device_matrix_view.value(), device_matrix_view.value_size() * sizeof(Value), hicMemcpyDeviceToHost); EXPECT(host_outer == expected_outer); EXPECT(host_inner == expected_inner); @@ -209,11 +209,29 @@ void check_matrix(const SparseMatrixStorage& m) { } void check_matrix(const SparseMatrixStorage& m) { - if (m.value().datatype() == atlas::make_datatype()) { - check_matrix(m); + if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); + } + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); } - else if (m.value().datatype() == atlas::make_datatype()) { - check_matrix(m); + else if (m.value().datatype() == atlas::make_datatype() && m.outer().datatype() == atlas::make_datatype()) { + check_matrix(m); } else { ATLAS_NOTIMPLEMENTED;