diff --git a/src/shogun/mathematics/linalg/internal/implementation/MatrixProduct.h b/src/shogun/mathematics/linalg/internal/implementation/MatrixProduct.h index 05ceae218a3..8f0340b6a18 100644 --- a/src/shogun/mathematics/linalg/internal/implementation/MatrixProduct.h +++ b/src/shogun/mathematics/linalg/internal/implementation/MatrixProduct.h @@ -85,9 +85,34 @@ struct matrix_product /** Scalar type */ typedef typename Matrix::Scalar T; + /** Return type */ + typedef SGMatrix ReturnType; + /** Eigen3 matrix type */ typedef Eigen::Matrix 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 A, SGMatrix 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 @@ -146,6 +171,31 @@ struct matrix_product /** Scalar type */ typedef typename Matrix::Scalar T; + /** Return type */ + typedef CGPUMatrix 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 A, CGPUMatrix 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 diff --git a/src/shogun/mathematics/linalg/internal/modules/Core.h b/src/shogun/mathematics/linalg/internal/modules/Core.h index 1e62f57f706..e1973b85fc5 100644 --- a/src/shogun/mathematics/linalg/internal/modules/Core.h +++ b/src/shogun/mathematics/linalg/internal/modules/Core.h @@ -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 @@ -93,6 +97,23 @@ void matrix_product(Matrix A, Matrix B, Matrix C, implementation::matrix_product::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,class Matrix> +typename implementation::matrix_product::ReturnType matrix_product(Matrix A, Matrix B, + bool transpose_A=false, bool transpose_B=false) +{ + return implementation::matrix_product::compute(A, B, transpose_A, transpose_B); +} + /** Performs the operation C = alpha*A - beta*B. Works for both matrices and vectors */ template ::backend,class Matrix> void subtract(Matrix A, Matrix B, Matrix C, diff --git a/tests/unit/mathematics/linalg/MatrixProduct_unittest.cc b/tests/unit/mathematics/linalg/MatrixProduct_unittest.cc index c7194e062ad..c9e2c9c4a74 100644 --- a/tests/unit/mathematics/linalg/MatrixProduct_unittest.cc +++ b/tests/unit/mathematics/linalg/MatrixProduct_unittest.cc @@ -51,37 +51,58 @@ TEST(MatrixProduct, eigen3_backend) SGMatrix A(3,3); SGMatrix B(3,3); SGMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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 A(3,3); + CGPUMatrix B(3,3); + + for (int32_t i=0; i<9; i++) + { + A[i] = i; + B[i] = 0.5*i; + } + + CGPUMatrix C = linalg::matrix_product(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 A(3,3); SGMatrix B(3,3); SGMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); } @@ -91,17 +112,17 @@ TEST(MatrixProduct, eigen3_backend_transpose_B) SGMatrix A(3,3); SGMatrix B(3,3); SGMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); } @@ -111,17 +132,17 @@ TEST(MatrixProduct, eigen3_backend_transpose_A_transpose_B) SGMatrix A(3,3); SGMatrix B(3,3); SGMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); } @@ -134,17 +155,36 @@ TEST(MatrixProduct, viennacl_backend) CGPUMatrix A(3,3); CGPUMatrix B(3,3); CGPUMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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 A(3,3); + SGMatrix B(3,3); + + for (int32_t i=0; i<9; i++) + { + A[i] = i; + B[i] = 0.5*i; + } + + SGMatrix C = linalg::matrix_product(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); } @@ -154,17 +194,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_A) CGPUMatrix A(3,3); CGPUMatrix B(3,3); CGPUMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); } @@ -174,17 +214,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_B) CGPUMatrix A(3,3); CGPUMatrix B(3,3); CGPUMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); } @@ -194,17 +234,17 @@ TEST(MatrixProduct, viennacl_backend_transpose_A_transpose_B) CGPUMatrix A(3,3); CGPUMatrix B(3,3); CGPUMatrix C(3,3); - + for (int32_t i=0; i<9; i++) { A[i] = i; B[i] = 0.5*i; } - + linalg::matrix_product(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); }