Skip to content

Commit

Permalink
Merge pull request #2802 from lambday/develop
Browse files Browse the repository at this point in the history 
Added backend-type/matrix-type independent matrix-product in linalg
  • Loading branch information
@lambday
lambday committed Apr 7, 2015
2 parents e61044a + 4549994 commit d65b1e7
Show file tree
Hide file tree
Showing 3 changed files with 144 additions and 33 deletions.
50 changes: 50 additions & 0 deletions src/shogun/mathematics/linalg/internal/implementation/MatrixProduct.h
View file
Expand Up @@ -85,9 +85,34 @@ struct matrix_product<Backend::EIGEN3, Matrix>
/** Scalar type */
typedef typename Matrix::Scalar T;

/** Return type */
typedef SGMatrix<T> ReturnType;

/** Eigen3 matrix type */
typedef Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> MatrixXt;

/** Performs matrix multiplication
*
* @param A First matrix
* @param B Second matrix
* @param transpose_A Whether to the transpose of A should be used instead of A
* @param transpose_B Whether to the transpose of B should be used instead of B
* @return Result of the operation
*/
static ReturnType compute(SGMatrix<T> A, SGMatrix<T> B,
bool transpose_A, bool transpose_B)
{
REQUIRE(A.matrix, "Matrix A is not initialized!\n");
REQUIRE(B.matrix, "Matrix B is not initialized!\n");
REQUIRE(A.num_cols == B.num_rows, "Number of columns for A (%d) and "
"number of rows for B (%d) should be equal!\n", A.num_cols, B.num_rows);

ReturnType retMatrix(A.num_rows, B.num_cols);
compute(A, B, retMatrix, transpose_A, transpose_B, true);

return retMatrix;
}

/** Performs matrix multiplication
*
* @param A First matrix
Expand Down Expand Up @@ -146,6 +171,31 @@ struct matrix_product<Backend::VIENNACL, Matrix>
/** Scalar type */
typedef typename Matrix::Scalar T;

/** Return type */
typedef CGPUMatrix<T> ReturnType;

/** Performs matrix multiplication
*
* @param A First matrix
* @param B Second matrix
* @param transpose_A Whether to the transpose of A should be used instead of A
* @param transpose_B Whether to the transpose of B should be used instead of B
* @return Result of the operation
*/
static ReturnType compute(CGPUMatrix<T> A, CGPUMatrix<T> B,
bool transpose_A, bool transpose_B)
{
REQUIRE(A.matrix, "Matrix A is not initialized!\n");
REQUIRE(B.matrix, "Matrix B is not initialized!\n");
REQUIRE(A.num_cols == B.num_rows, "Number of columns for A (%d) and "
"number of rows for B (%d) should be equal!\n", A.num_cols, B.num_rows);

ReturnType retMatrix(A.num_rows, B.num_cols);
compute(A, B, retMatrix, transpose_A, transpose_B, true);

return retMatrix;
}

/** Performs matrix multiplication
*
* @param A First matrix
Expand Down
23 changes: 22 additions & 1 deletion src/shogun/mathematics/linalg/internal/modules/Core.h
View file
Expand Up @@ -76,7 +76,11 @@ void scale(Matrix A, typename Matrix::Scalar alpha)
}

#ifdef HAVE_LINALG_LIB
/** Performs matrix multiplication
/** Performs matrix multiplication. This version should be used for backend specific
* code requirements. For example, use this with CGPUMatrix and explicitly set
* ViennaCL backend, or SGMatrix and explicitly set Eigen3 backend. If matrix-type/
* backend-type independent code is desired, use the version that does not support
* preallocated result matrix but returns the result in a newly created matrix instead.
*
* @param A First matrix
* @param B Second matrix
Expand All @@ -93,6 +97,23 @@ void matrix_product(Matrix A, Matrix B, Matrix C,
implementation::matrix_product<backend, Matrix>::compute(A, B, C, transpose_A, transpose_B, overwrite);
}

/** Performs matrix multiplication. This version returns the result in a newly
* created matrix. If matrix-product is desired that will work irrespective of the
* backend and the matrix type used, then this method should be used.
*
* @param A First matrix
* @param B Second matrix
* @param transpose_A Whether to the transpose of A should be used instead of A
* @param transpose_B Whether to the transpose of B should be used instead of B
* @return Result of the operation
*/
template <Backend backend=linalg_traits<Core>::backend,class Matrix>
typename implementation::matrix_product<backend,Matrix>::ReturnType matrix_product(Matrix A, Matrix B,
bool transpose_A=false, bool transpose_B=false)
{
return implementation::matrix_product<backend, Matrix>::compute(A, B, transpose_A, transpose_B);
}

/** Performs the operation C = alpha*A - beta*B. Works for both matrices and vectors */
template <Backend backend=linalg_traits<Core>::backend,class Matrix>
void subtract(Matrix A, Matrix B, Matrix C,
Expand Down
104 changes: 72 additions & 32 deletions tests/unit/mathematics/linalg/MatrixProduct_unittest.cc
View file
Expand Up @@ -51,37 +51,58 @@ TEST(MatrixProduct, eigen3_backend)
SGMatrix<float64_t> A(3,3);
SGMatrix<float64_t> B(3,3);
SGMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::EIGEN3>(A, B, C);


float64_t ref[] = {7.5, 9, 10.5, 21, 27, 33, 34.5, 45, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}

#ifdef HAVE_VIENNACL
TEST(MatrixProduct, GPUMatrix_eigen3_backend)
{
CGPUMatrix<float64_t> A(3,3);
CGPUMatrix<float64_t> B(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

CGPUMatrix<float64_t> C = linalg::matrix_product<linalg::Backend::EIGEN3>(A, B);

float64_t ref[] = {7.5, 9, 10.5, 21, 27, 33, 34.5, 45, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
#endif // HAVE_VIENNACL

TEST(MatrixProduct, eigen3_backend_transpose_A)
{
SGMatrix<float64_t> A(3,3);
SGMatrix<float64_t> B(3,3);
SGMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::EIGEN3>(A, B, C, true);

float64_t ref[] = {2.5, 7, 11.5, 7, 25, 43, 11.5, 43, 74.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -91,17 +112,17 @@ TEST(MatrixProduct, eigen3_backend_transpose_B)
SGMatrix<float64_t> A(3,3);
SGMatrix<float64_t> B(3,3);
SGMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::EIGEN3>(A, B, C, false, true);

float64_t ref[] = {22.5, 27, 31.5, 27, 33, 39, 31.5, 39, 46.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -111,17 +132,17 @@ TEST(MatrixProduct, eigen3_backend_transpose_A_transpose_B)
SGMatrix<float64_t> A(3,3);
SGMatrix<float64_t> B(3,3);
SGMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::EIGEN3>(A, B, C, true, true);

float64_t ref[] = {7.5, 21, 34.5, 9, 27, 45, 10.5, 33, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -134,17 +155,36 @@ TEST(MatrixProduct, viennacl_backend)
CGPUMatrix<float64_t> A(3,3);
CGPUMatrix<float64_t> B(3,3);
CGPUMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::VIENNACL>(A, B, C);


float64_t ref[] = {7.5, 9, 10.5, 21, 27, 33, 34.5, 45, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}

TEST(MatrixProduct, SGMatrix_viennacl_backend)
{
SGMatrix<float64_t> A(3,3);
SGMatrix<float64_t> B(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

SGMatrix<float64_t> C = linalg::matrix_product<linalg::Backend::VIENNACL>(A, B);

float64_t ref[] = {7.5, 9, 10.5, 21, 27, 33, 34.5, 45, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -154,17 +194,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_A)
CGPUMatrix<float64_t> A(3,3);
CGPUMatrix<float64_t> B(3,3);
CGPUMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::VIENNACL>(A, B, C, true);

float64_t ref[] = {2.5, 7, 11.5, 7, 25, 43, 11.5, 43, 74.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -174,17 +214,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_B)
CGPUMatrix<float64_t> A(3,3);
CGPUMatrix<float64_t> B(3,3);
CGPUMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::VIENNACL>(A, B, C, false, true);

float64_t ref[] = {22.5, 27, 31.5, 27, 33, 39, 31.5, 39, 46.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand All @@ -194,17 +234,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_A_transpose_B)
CGPUMatrix<float64_t> A(3,3);
CGPUMatrix<float64_t> B(3,3);
CGPUMatrix<float64_t> C(3,3);

for (int32_t i=0; i<9; i++)
{
A[i] = i;
B[i] = 0.5*i;
}

linalg::matrix_product<linalg::Backend::VIENNACL>(A, B, C, true, true);

float64_t ref[] = {7.5, 21, 34.5, 9, 27, 45, 10.5, 33, 55.5};

for (int32_t i=0; i<9; i++)
EXPECT_NEAR(ref[i], C[i], 1e-15);
}
Expand Down

0 comments on commit d65b1e7

Please sign in to comment.