diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc index 1a38a9ca6..4a88f5be2 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.cc @@ -14,6 +14,7 @@ #include #include "eckit/config/LocalConfiguration.h" +#include "eckit/linalg/types.h" #include "atlas/array/ArrayView.h" #include "atlas/array/helpers/ArrayForEach.h" @@ -53,7 +54,6 @@ using Complex = SphericalVector::Complex; template using SparseMatrix = SphericalVector::SparseMatrix; -using RealMatrixMap = Eigen::Map>; using ComplexTriplets = std::vector>; using RealTriplets = std::vector>; @@ -61,12 +61,6 @@ using EckitMatrix = eckit::linalg::SparseMatrix; namespace { -RealMatrixMap makeMatrixMap(const EckitMatrix& baseMatrix) { - return RealMatrixMap(baseMatrix.rows(), baseMatrix.cols(), - baseMatrix.nonZeros(), baseMatrix.outer(), - baseMatrix.inner(), baseMatrix.data()); -} - template auto getInnerIt(const Matrix& matrix, typename Matrix::Index k) { return typename Matrix::InnerIterator(matrix, k); @@ -148,7 +142,11 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto nRows = matrix().rows(); const auto nCols = matrix().cols(); const auto nNonZeros = matrix().nonZeros(); - const auto baseWeights = makeMatrixMap(matrix()); + const auto outerIndices = matrix().outer(); + const auto innerIndices = matrix().inner(); + const auto baseWeights = matrix().data(); + + using Index = eckit::linalg::Index; // Note: need to store copy of weights as Eigen3 sorts compressed rows by j // whereas eckit does not. @@ -162,24 +160,27 @@ void SphericalVector::do_setup(const FunctionSpace& source, const auto unitSphere = geometry::UnitSphere{}; - const auto setWeights = [&](auto i, auto j, const auto& baseWeight) { - const auto sourceLonLat = - PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); - const auto targetLonLat = - PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); + for (auto i = Index{}; i < nRows; ++i) { + for (auto k = outerIndices[i]; k < outerIndices[i + 1]; ++k) { + const auto j = innerIndices[k]; + const auto baseWeight = baseWeights[k]; - const auto alpha = unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); + const auto sourceLonLat = + PointLonLat(sourceLonLats(j, 0), sourceLonLats(j, 1)); + const auto targetLonLat = + PointLonLat(targetLonLats(i, 0), targetLonLats(i, 1)); - const auto deltaAlpha = - (alpha.first - alpha.second) * util::Constants::degreesToRadians(); + const auto alpha = + unitSphere.greatCircleCourse(sourceLonLat, targetLonLat); - const auto idx = std::distance(baseWeights.valuePtr(), &baseWeight); + const auto deltaAlpha = + (alpha.first - alpha.second) * util::Constants::degreesToRadians(); - complexTriplets[idx] = {int(i), int(j), std::polar(baseWeight, deltaAlpha)}; - realTriplets[idx] = {int(i), int(j), baseWeight}; - }; + complexTriplets[k] = {i, j, std::polar(baseWeight, deltaAlpha)}; + realTriplets[k] = {i, j, baseWeight}; + } + } - sparseMatrixForEach(setWeights, baseWeights); complexWeights_->setFromTriplets(complexTriplets.begin(), complexTriplets.end()); realWeights_->setFromTriplets(realTriplets.begin(), realTriplets.end()); diff --git a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h index 91ce2a1a2..b4e809999 100644 --- a/src/atlas/interpolation/method/sphericalvector/SphericalVector.h +++ b/src/atlas/interpolation/method/sphericalvector/SphericalVector.h @@ -18,7 +18,6 @@ #include "atlas/functionspace/FunctionSpace.h" #include "atlas/interpolation/method/Method.h" -#include "atlas/linalg/sparse.h" namespace atlas {