Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 2 additions & 4 deletions Sofa/framework/LinearAlgebra/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -37,8 +37,7 @@ set(HEADER_FILES
${SOFALINEARALGEBRASRC_ROOT}/RotationMatrix.h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrix.h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[CompressedRowSparseMatrix].h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[EigenSparseMatrix].h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.inl
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder.h
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder[EigenSparseMatrix].h
${SOFALINEARALGEBRASRC_ROOT}/TriangularSystemSolver.h
Expand All @@ -59,8 +58,7 @@ set(SOURCE_FILES
${SOFALINEARALGEBRASRC_ROOT}/FullMatrix.cpp
${SOFALINEARALGEBRASRC_ROOT}/FullVector.cpp
${SOFALINEARALGEBRASRC_ROOT}/RotationMatrix.cpp
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[CompressedRowSparseMatrix].cpp
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct[EigenSparseMatrix].cpp
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixProduct.cpp
${SOFALINEARALGEBRASRC_ROOT}/SparseMatrixStorageOrder[EigenSparseMatrix].cpp
${SOFALINEARALGEBRASRC_ROOT}/init.cpp
)
Expand Down
1 change: 1 addition & 0 deletions Sofa/framework/LinearAlgebra/Testing/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ project(Sofa.LinearAlgebra.Testing LANGUAGES CXX)
set(HEADER_FILES
src/Sofa.LinearAlgebra.Testing/BaseMatrix_test.h
src/Sofa.LinearAlgebra.Testing/SparseMatrixTest.h
src/Sofa.LinearAlgebra.Testing/SparseMatrixProduct_test.h
)

add_library(${PROJECT_NAME} INTERFACE)
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
/******************************************************************************
* SOFA, Simulation Open-Framework Architecture *
* (c) 2006 INRIA, USTL, UJF, CNRS, MGH *
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by *
* the Free Software Foundation; either version 2.1 of the License, or (at *
* your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*******************************************************************************
* Authors: The SOFA Team and external contributors (see Authors.txt) *
* *
* Contact information: contact@sofa-framework.org *
******************************************************************************/
#pragma once
#include <gtest/gtest.h>
#include <sofa/testing/NumericTest.h>
#include <Sofa.LinearAlgebra.Testing/SparseMatrixTest.h>

#include <sofa/linearalgebra/BaseMatrix.h>
#include <Eigen/Sparse>
#include <sofa/helper/random.h>

namespace sofa::linearalgebra::testing
{

template<class T>
struct SparseMatrixProductInit
{
static void init(T& product)
{
SOFA_UNUSED(product);
};

static void cleanup(T& product)
{
SOFA_UNUSED(product);
}
};

/**
* Test the class SparseMatrixProduct
* The class is designed to use the templates defined in TestSparseMatrixProductTraits
* The type of matrix can be any of the types supported by SparseMatrixProduct.
*/
template <class T>
struct TestSparseMatrixProduct : public sofa::testing::SparseMatrixTest<typename T::LhsScalar>
{
using LHSMatrix = typename T::LhsCleaned;
using RHSMatrix = typename T::RhsCleaned;
using Real = typename T::LhsScalar;
using Base = sofa::testing::SparseMatrixTest<Real>;
using Base::generateRandomSparseMatrix;
using Base::copyFromEigen;
using Base::compareSparseMatrix;

bool checkMatrix(typename LHSMatrix::Index nbRowsA, typename LHSMatrix::Index nbColsA, typename RHSMatrix::Index nbColsB, Real sparsity)
{
Eigen::SparseMatrix<Real, Eigen::RowMajor> eigen_a;
Eigen::SparseMatrix<Real, Eigen::ColMajor> eigen_b;

generateRandomSparseMatrix(eigen_a, nbRowsA, nbColsA, sparsity);
generateRandomSparseMatrix(eigen_b, nbColsA, nbColsB, sparsity);

LHSMatrix A;
RHSMatrix B;
copyFromEigen(A, eigen_a);
copyFromEigen(B, eigen_b);

EXPECT_GT(eigen_a.outerSize(), 0);
EXPECT_GT(eigen_b.outerSize(), 0);

Eigen::SparseMatrix<Real, Eigen::RowMajor> eigen_c = eigen_a * eigen_b;

EXPECT_EQ(eigen_c.rows(), nbRowsA);
EXPECT_EQ(eigen_c.cols(), nbColsB);
EXPECT_GT(eigen_c.outerSize(), 0); //to make sure that there are non-zero values in the result matrix

T product(&A, &B);
SparseMatrixProductInit<T>::init(product);
product.computeProduct();

EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult()));

//the second time computeProduct is called uses the faster algorithm
product.computeProduct();
EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult()));

// force re-computing the intersection
product.computeProduct(true);
EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult()));

//modify the values of A, but not its pattern
for (int i = 0; i < eigen_a.nonZeros(); ++i)
{
eigen_a.valuePtr()[i] = static_cast<SReal>(sofa::helper::drand(1));
}

eigen_c = eigen_a * eigen_b; //result is updated using the regular matrix product
copyFromEigen(A, eigen_a);

product.m_lhs = &A;
product.computeProduct(); //intersection is already computed: uses the faster algorithm
EXPECT_TRUE(compareSparseMatrix(eigen_c, product.getProductResult()));

SparseMatrixProductInit<T>::cleanup(product);

return true;
}
};

TYPED_TEST_SUITE_P(TestSparseMatrixProduct);

TYPED_TEST_P(TestSparseMatrixProduct, squareMatrix )
{
EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 1. / 5. ) );
EXPECT_TRUE( this->checkMatrix( 5, 5, 5, 3. / 5. ) );

EXPECT_TRUE( this->checkMatrix( 100, 1000, 1000, 1. / 1000. ) );
EXPECT_TRUE( this->checkMatrix( 1000, 1000, 1000, 20. / 1000. ) );
// EXPECT_TRUE( this->checkMatrix( 1000, 1000, 1000, 150. / 1000. ) );

EXPECT_TRUE( this->checkMatrix( 20, 20, 20, 1. ) );
}

TYPED_TEST_P(TestSparseMatrixProduct, rectangularMatrix )
{
EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 1. / 5. ) );
EXPECT_TRUE( this->checkMatrix( 5, 10, 7, 3. / 5. ) );

EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 1. / 5. ) );
EXPECT_TRUE( this->checkMatrix( 10, 5, 7, 3. / 5. ) );

EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 1. / 1000. ) );
EXPECT_TRUE( this->checkMatrix( 1000, 3000, 2000, 20. / 1000. ) );

EXPECT_TRUE( this->checkMatrix( 20, 30, 10, 1. ) );
}

REGISTER_TYPED_TEST_SUITE_P(TestSparseMatrixProduct,
squareMatrix,rectangularMatrix);

}
Original file line number Diff line number Diff line change
Expand Up @@ -68,12 +68,11 @@ struct SparseMatrixTest : public virtual NumericTest<TReal>
eigenMatrix.setFromTriplets(first, last);
}

static void copyFromEigen(Eigen::SparseMatrix<TReal>& dst, const Eigen::SparseMatrix<TReal>& src)
{
dst = src;
}

static void copyFromEigen(Eigen::SparseMatrix<TReal, Eigen::RowMajor>& dst, const Eigen::SparseMatrix<TReal, Eigen::RowMajor>& src)
template<
typename _DstScalar, int _DstOptions, typename _DstStorageIndex,
typename _SrcScalar, int _SrcOptions, typename _SrcStorageIndex
>
static void copyFromEigen(Eigen::SparseMatrix<_DstScalar, _DstOptions, _DstStorageIndex>& dst, const Eigen::SparseMatrix<_SrcScalar, _SrcOptions, _SrcStorageIndex>& src)
{
dst = src;
}
Expand All @@ -93,48 +92,26 @@ struct SparseMatrixTest : public virtual NumericTest<TReal>
}
}

static bool compareSparseMatrix(const Eigen::SparseMatrix<TReal>& A, const Eigen::SparseMatrix<TReal>& B)
template<
typename _AScalar, int _AOptions, typename _AStorageIndex,
typename _BScalar, int _BOptions, typename _BStorageIndex
>
static bool compareSparseMatrix(const Eigen::SparseMatrix<_AScalar, _AOptions, _AStorageIndex>& A, const Eigen::SparseMatrix<_BScalar, _BOptions, _BStorageIndex>& B)
{
return compareEigenSparseMatrix(A, B);
}

static bool compareEigenSparseMatrix(const Eigen::SparseMatrix<TReal>& A, const Eigen::SparseMatrix<TReal>& B)
template<
typename _AScalar, int _AOptions, typename _AStorageIndex,
typename _BScalar, int _BOptions, typename _BStorageIndex
>
static bool compareEigenSparseMatrix(const Eigen::SparseMatrix<_AScalar, _AOptions, _AStorageIndex>& A, const Eigen::SparseMatrix<_BScalar, _BOptions, _BStorageIndex>& B)
{
if (A.outerSize() != B.outerSize())
return false;

for (int k = 0; k < A.outerSize(); ++k)
{
sofa::type::vector<Eigen::Triplet<TReal> > triplets_a, triplets_b;
for (typename Eigen::SparseMatrix<TReal>::InnerIterator it(A, k); it; ++it)
{
triplets_a.emplace_back(it.row(), it.col(), it.value());
}
for (typename Eigen::SparseMatrix<TReal>::InnerIterator it(B, k); it; ++it)
{
triplets_b.emplace_back(it.row(), it.col(), it.value());
}

if (triplets_a.size() != triplets_b.size())
return false;

for (size_t i = 0 ; i < triplets_a.size(); ++i)
{
const auto& a = triplets_a[i];
const auto& b = triplets_b[i];

if (a.row() != b.row() || a.col() != b.col() || !NumericTest<TReal>::isSmall(a.value() - b.value(), 1))
{
// ret = false;
return false;
}
}
}
return true;
return A.isApprox(B);
}


};


} // namespace sofa::testing
} // namespace sofa::testing
Original file line number Diff line number Diff line change
Expand Up @@ -19,17 +19,8 @@
* *
* Contact information: contact@sofa-framework.org *
******************************************************************************/
#pragma once

#include <sofa/linearalgebra/SparseMatrixProduct.h>
#include <sofa/linearalgebra/CompressedRowSparseMatrix.h>
#define SOFA_LINEARALGEBRA_SPARSEMATRIXPRODUCT_DEFINITION

namespace sofa::linearalgebra
{

#if !defined(SOFA_LINEARAGEBRA_SPARSEMATRIXPRODUCT_COMPRESSEDROWSPARSEMATRIX_CPP)
extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct<CompressedRowSparseMatrix<float> >;
extern template class SOFA_LINEARALGEBRA_API SparseMatrixProduct<CompressedRowSparseMatrix<double> >;
#endif

} //namespace sofa::linearalgebra
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@alxbilger is there a sense to keep this empty file?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

File removed in #4599

Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@

#include <sofa/linearalgebra/config.h>
#include <utility>

#include <sofa/type/vector_T.h>
#include <sofa/type/vector.h>
#include <Eigen/Core>

namespace sofa::linearalgebra
{
Expand All @@ -44,81 +44,63 @@ namespace sofa::linearalgebra
* and
* Saupin, G., 2008. Vers la simulation interactive réaliste de corps déformables virtuels (Doctoral dissertation, Lille 1).
*/
template<class TMatrix>
template<class Lhs, class Rhs, class ResultType>
class SparseMatrixProduct
{
public:

using LhsCleaned = std::decay_t<Lhs>;
using RhsCleaned = std::decay_t<Rhs>;
using ResultCleaned = std::decay_t<ResultType>;

using LhsScalar = typename LhsCleaned::Scalar;
using RhsScalar = typename RhsCleaned::Scalar;
using ResultScalar = typename ResultCleaned::Scalar;

/// Left side of the product A*B
const TMatrix* matrixA { nullptr };
const LhsCleaned* m_lhs { nullptr };
/// Right side of the product A*B
const TMatrix* matrixB { nullptr };
const RhsCleaned* m_rhs { nullptr };

using Index = Eigen::Index;

using ProductResult = ResultType;


void computeProduct(bool forceComputeIntersection = false);
void computeRegularProduct();

const TMatrix& getProductResult() const { return matrixC; }
[[nodiscard]] const ResultType& getProductResult() const { return m_productResult; }

void invalidateIntersection();

SparseMatrixProduct(TMatrix* a, TMatrix* b) : matrixA(a), matrixB(b) {}
SparseMatrixProduct(Lhs* lhs, Rhs* rhs) : m_lhs(lhs), m_rhs(rhs) {}
SparseMatrixProduct() = default;
virtual ~SparseMatrixProduct() = default;

struct Intersection
{
// Two indices: the first for the values vector of the matrix A, the second for the values vector of the matrix B
using PairIndex = std::pair<typename TMatrix::Index, typename TMatrix::Index>;
using PairIndex = std::pair<Index, Index>;
// A list of pairs of indices
using ListPairIndex = type::vector<PairIndex>;
using ListPairIndex = sofa::type::vector<PairIndex>;
/// Each element of this vector gives the list of values from matrix A and B to multiply together and accumulate
/// them into the matrix C at the same location in the values vector
using ValuesIntersection = type::vector<ListPairIndex>;
using ValuesIntersection = sofa::type::vector<ListPairIndex>;

ValuesIntersection intersection;
};

protected:
TMatrix matrixC; /// Result of A*B
ProductResult m_productResult; /// Result of LHS * RHS

bool m_hasComputedIntersection { false };
void computeIntersection();
void computeProductFromIntersection();
virtual void computeProductFromIntersection();

Intersection m_intersectionAB;

};

template <class TMatrix>
void SparseMatrixProduct<TMatrix>::computeProduct(bool forceComputeIntersection)
{
if (forceComputeIntersection)
{
m_hasComputedIntersection = false;
}

if (m_hasComputedIntersection == false)
{
computeIntersection();
m_hasComputedIntersection = true;
}
computeProductFromIntersection();
}

template <class TMatrix>
void SparseMatrixProduct<TMatrix>::computeIntersection()
{

}

template <class TMatrix>
void SparseMatrixProduct<TMatrix>::computeProductFromIntersection()
{
matrixC = (*matrixA) * (*matrixB);
}

template <class TMatrix>
void SparseMatrixProduct<TMatrix>::invalidateIntersection()
{
m_hasComputedIntersection = false;
}

}// sofa::linearalgebra
}// sofa::linearalgebra
Loading