Merge pull request #3387 from OXPHOS/linalg_add

LinalgRefactor - SGMatrix - add

src/shogun/mathematics/linalg/LinalgBackendBase.h

@@ -94,6 +94,7 @@ class LinalgBackendBase
 		return 0; \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_ADD
 
 	/**
@@ -107,6 +108,7 @@ class LinalgBackendBase
 		SG_SNOTIMPLEMENTED; \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_IN_PLACE_ADD
 
 	/**

src/shogun/mathematics/linalg/LinalgBackendEigen.h

@@ -82,6 +82,7 @@ class LinalgBackendEigen : public LinalgBackendBase
 		return add_impl(a, b, alpha, beta); \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_ADD
 
 	/** Implementation of @see LinalgBackendBase::add */
@@ -91,6 +92,7 @@ class LinalgBackendEigen : public LinalgBackendBase
 		add_impl(a, b, alpha, beta, result); \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_IN_PLACE_ADD
 
 	/** Implementation of @see LinalgBackendBase::dot */
@@ -193,6 +195,19 @@ class LinalgBackendEigen : public LinalgBackendBase
 		return c;
 	}
 
+	/** Eigen3 matrix C = alpha*A + beta*B method */
+	template <typename T>
+	SGMatrix<T> add_impl(const SGMatrix<T>& a, const SGMatrix<T>& b, T alpha, T beta) const
+	{
+		SGMatrix<T> c(a.num_rows, a.num_cols);
+		typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
+		typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
+		typename SGMatrix<T>::EigenMatrixXtMap c_eig = c;
+
+		c_eig = alpha * a_eig + beta * b_eig;
+		return c;
+	}
+
 	/** Eigen3 vector result = alpha*A + beta*B method */
 	template <typename T>
 	void add_impl(SGVector<T>& a, SGVector<T>& b, T alpha, T beta, SGVector<T>& result) const
@@ -204,6 +219,17 @@ class LinalgBackendEigen : public LinalgBackendBase
 		result_eig = alpha * a_eig + beta * b_eig;
 	}
 
+	/** Eigen3 matrix result = alpha*A + beta*B method */
+	template <typename T>
+	void add_impl(SGMatrix<T>& a, SGMatrix<T>& b, T alpha, T beta, SGMatrix<T>& result) const
+	{
+		typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
+		typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
+		typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
+
+		result_eig = alpha * a_eig + beta * b_eig;
+	}
+
 	/** Eigen3 vector dot-product method */
 	template <typename T>
 	T dot_impl(const SGVector<T>& a, const SGVector<T>& b) const

src/shogun/mathematics/linalg/LinalgBackendViennacl.h

@@ -72,6 +72,7 @@ class LinalgBackendViennaCL : public LinalgBackendGPUBase
 		return add_impl(a, b, alpha, beta); \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_ADD
 
 	#define BACKEND_GENERIC_IN_PLACE_ADD(Type, Container) \
@@ -80,6 +81,7 @@ class LinalgBackendViennaCL : public LinalgBackendGPUBase
 		add_impl(a, b, alpha, beta, result); \
 	}
 	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGVector)
+	DEFINE_FOR_ALL_PTYPE(BACKEND_GENERIC_IN_PLACE_ADD, SGMatrix)
 	#undef BACKEND_GENERIC_ADD
 
 	/** Implementation of @see LinalgBackendBase::dot */
@@ -171,6 +173,19 @@ class LinalgBackendViennaCL : public LinalgBackendGPUBase
 		return SGVector<T>(c_gpu, a.size());
 	}
 
+	/** ViennaCL matrix C = alpha*A + beta*B method */
+	template <typename T>
+	SGMatrix<T> add_impl(const SGMatrix<T>& a, const SGMatrix<T>& b, T alpha, T beta) const
+	{
+		GPUMemoryViennaCL<T>* a_gpu = cast_to_viennacl(a);
+		GPUMemoryViennaCL<T>* b_gpu = cast_to_viennacl(b);
+		GPUMemoryViennaCL<T>* c_gpu = new GPUMemoryViennaCL<T>(a.size());
+
+		c_gpu->data_matrix(a.num_rows, a.num_cols) = alpha * a_gpu->data_matrix(a.num_rows, a.num_cols)
+			+ beta * b_gpu->data_matrix(b.num_rows, b.num_cols);
+		return SGMatrix<T>(c_gpu, a.num_rows, a.num_cols);
+	}
+
 	/** ViennaCL vector result = alpha*A + beta*B method */
 	template <typename T>
 	void add_impl(SGVector<T>& a, SGVector<T>& b, T alpha, T beta, SGVector<T>& result) const
@@ -183,6 +198,19 @@ class LinalgBackendViennaCL : public LinalgBackendGPUBase
 			alpha * a_gpu->data_vector(a.size()) + beta * b_gpu->data_vector(b.size());
 	}
 
+	/** ViennaCL matrix result = alpha*A + beta*B method */
+	template <typename T>
+	void add_impl(SGMatrix<T>& a, SGMatrix<T>& b, T alpha, T beta, SGMatrix<T>& result) const
+	{
+		GPUMemoryViennaCL<T>* a_gpu = cast_to_viennacl(a);
+		GPUMemoryViennaCL<T>* b_gpu = cast_to_viennacl(b);
+		GPUMemoryViennaCL<T>* result_gpu = cast_to_viennacl(result);
+
+		result_gpu->data_matrix(a.num_rows, a.num_cols) =
+			alpha * a_gpu->data_matrix(a.num_rows, a.num_cols)
+			+ beta * b_gpu->data_matrix(b.num_rows, b.num_cols);
+	}
+
 	/** ViennaCL vector dot-product method. */
 	template <typename T>
 	T dot_impl(const SGVector<T>& a, const SGVector<T>& b) const

src/shogun/mathematics/linalg/LinalgNamespace.h

@@ -105,13 +105,31 @@ LinalgBackendBase* infer_backend(const Container<T>& a, const Container<T>& b)
  * @param beta constant to be multiplied by the second vector
  * @return The result vector
  */
-template <typename T, template <typename> class Container>
-Container<T> add(const Container<T>& a, const Container<T>& b, T alpha=1, T beta=1)
+template <typename T>
+SGVector<T> add(const SGVector<T>& a, const SGVector<T>& b, T alpha=1, T beta=1)
 {
 	REQUIRE(a.vlen == b.vlen, "Length of vector a (%d) doesn't match vector b (%d).\n", a.vlen, b.vlen);
 	return infer_backend(a, b)->add(a, b, alpha, beta);
 }
 
+/**
+ * Performs the operation C = alpha*A + beta*B.
+ * @param a first matrix
+ * @param b second matrix
+ * @param alpha constant to be multiplied by the first matrix
+ * @param beta constant to be multiplied by the second matrix
+ * @return the result matrix
+ */
+template <typename T>
+SGMatrix<T> add(const SGMatrix<T>& a, const SGMatrix<T>& b, T alpha=1, T beta=1)
+{
+	REQUIRE((a.num_rows == b.num_rows), "Number of rows of matrix a (%d) must match matrix b (%d).\n",
+			a.num_rows, b.num_rows);
+	REQUIRE((a.num_cols == b.num_cols), "Number of columns of matrix a (%d) must match matrix b (%d).\n",
+			a.num_cols, b.num_cols);
+	return infer_backend(a, b)->add(a, b, alpha, beta);
+}
+
 /**
  * Performs the operation result = alpha*a + beta*b.
  *
@@ -135,6 +153,31 @@ void add(SGVector<T>& a, SGVector<T>& b, SGVector<T>& result, T alpha=1, T beta=
 	infer_backend(a, b)->add(a, b, alpha, beta, result);
 }
 
+/**
+ * Performs the operation result = alpha*a + beta*b.
+ *
+ * @param a first matrix
+ * @param b second matrix
+ * @param result the matrix that saves the result
+ * @param alpha constant to be multiplied by the first matrix
+ * @param beta constant to be multiplied by the second matrix
+ */
+template <typename T>
+void add(SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result, T alpha=1, T beta=1)
+{
+	REQUIRE((a.num_rows == b.num_rows), "Number of rows of matrix a (%d) must match matrix b (%d).\n",
+		a.num_rows, b.num_rows);
+	REQUIRE((a.num_cols == b.num_cols), "Number of columns of matrix a (%d) must match matrix b (%d).\n",
+		a.num_cols, b.num_cols);
+
+	REQUIRE(!(result.on_gpu()^a.on_gpu()),
+		"Cannot operate with matrix result on_gpu (%d) and matrix a on_gpu (%d).\n", result.on_gpu(), a.on_gpu());
+	REQUIRE(!(result.on_gpu()^b.on_gpu()),
+		"Cannot operate with matrix result on_gpu (%d) and matrix b on_gpu (%d).\n", result.on_gpu(), b.on_gpu());
+
+	infer_backend(a, b)->add(a, b, alpha, beta, result);
+}
+
 /**
  * Vector dot-product that works with generic vectors.
  *

tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc

@@ -26,7 +26,28 @@ TEST(LinalgBackendEigen, SGVector_add)
 		EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15);
 }
 
-TEST(LinalgBackendEigen, add_in_place)
+TEST(LinalgBackendEigen, SGMatrix_add)
+{
+	const float64_t alpha = 0.3;
+	const float64_t beta = -1.5;
+	const index_t nrows = 2, ncols = 3;
+
+	SGMatrix<float64_t> A(nrows, ncols);
+	SGMatrix<float64_t> B(nrows, ncols);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+	{
+		A[i] = i;
+		B[i] = 0.5*i;
+	}
+
+	auto result = add(A, B, alpha, beta);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+		EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15);
+}
+
+TEST(LinalgBackendEigen, SGVector_add_in_place)
 {
 	const float64_t alpha = 0.3;
 	const float64_t beta = -1.5;
@@ -46,6 +67,29 @@ TEST(LinalgBackendEigen, add_in_place)
 		EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15);
 }
 
+TEST(LinalgBackendEigen, SGMatrix_add_in_place)
+{
+	const float64_t alpha = 0.3;
+	const float64_t beta = -1.5;
+	const index_t nrows = 2, ncols = 3;
+
+	SGMatrix<float64_t> A(nrows, ncols);
+	SGMatrix<float64_t> B(nrows, ncols);
+	SGMatrix<float64_t> C(nrows, ncols);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+	{
+		A[i] = i;
+		B[i] = 0.5*i;
+		C[i] = i;
+	}
+
+	add(A, B, A, alpha, beta);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+		EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15);
+}
+
 TEST(LinalgBackendEigen, SGVector_dot)
 {
 	const index_t size = 3;

tests/unit/mathematics/linalg/operations/Viennacl_operations_unittest.cc

@@ -13,11 +13,11 @@ TEST(LinalgBackendViennaCL, SGVector_add)
 {
 	sg_linalg->set_gpu_backend(new LinalgBackendViennaCL());
 
-	const float64_t alpha = 0.3;
-	const float64_t beta = -1.5;
+	const float32_t alpha = 0.3;
+	const float32_t beta = -1.5;
 
-	SGVector<float64_t> A(9), A_gpu;
-	SGVector<float64_t> B(9), B_gpu;
+	SGVector<float32_t> A(9), A_gpu;
+	SGVector<float32_t> B(9), B_gpu;
 
 	for (index_t i = 0; i < 9; ++i)
 	{
@@ -28,13 +28,40 @@ TEST(LinalgBackendViennaCL, SGVector_add)
 	A_gpu = to_gpu(A);
 	B_gpu = to_gpu(B);
 
-	auto result = add(A, B, alpha, beta);
+	auto result_gpu = add(A_gpu, B_gpu, alpha, beta);
+	auto result = from_gpu(result_gpu);
 
 	for (index_t i = 0; i < 9; ++i)
 		EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15);
 }
 
-TEST(LinalgBackendViennaCL, add_in_place)
+TEST(LinalgBackendViennaCL, SGMatrix_add)
+{
+	sg_linalg->set_gpu_backend(new LinalgBackendViennaCL());
+
+	const float32_t alpha = 0.3;
+	const float32_t beta = -1.5;
+	const index_t nrows = 2, ncols = 3;
+
+	SGMatrix<float32_t> A(nrows, ncols), A_gpu;
+	SGMatrix<float32_t> B(nrows, ncols), B_gpu;
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+	{
+		A[i] = i;
+		B[i] = 0.5*i;
+	}
+	A_gpu = to_gpu(A);
+	B_gpu = to_gpu(B);
+
+	auto result_gpu = add(A_gpu, B_gpu, alpha, beta);
+	auto result = from_gpu(result_gpu);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+		EXPECT_NEAR(alpha*A[i]+beta*B[i], result[i], 1e-15);
+}
+
+TEST(LinalgBackendViennaCL, SGVector_add_in_place)
 {
 	sg_linalg->set_gpu_backend(new LinalgBackendViennaCL());
 
@@ -60,6 +87,34 @@ TEST(LinalgBackendViennaCL, add_in_place)
 		EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15);
 }
 
+TEST(LinalgBackendViennaCL, SGMatrix_add_in_place)
+{
+	sg_linalg->set_gpu_backend(new LinalgBackendViennaCL());
+
+	const float32_t alpha = 0.3;
+	const float32_t beta = -1.5;
+	const index_t nrows = 2, ncols = 3;
+
+	SGMatrix<float32_t> A(nrows, ncols), A_gpu;
+	SGMatrix<float32_t> B(nrows, ncols), B_gpu;
+	SGMatrix<float32_t> C(nrows, ncols);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+	{
+		A[i] = i;
+		B[i] = 0.5*i;
+		C[i] = i;
+	}
+	A_gpu = to_gpu(A);
+	B_gpu = to_gpu(B);
+
+	add(A_gpu, B_gpu, A_gpu, alpha, beta);
+	A = from_gpu(A_gpu);
+
+	for (index_t i = 0; i < nrows*ncols; ++i)
+		EXPECT_NEAR(alpha*C[i]+beta*B[i], A[i], 1e-15);
+}
+
 TEST(LinalgBackendViennaCL, SGVector_dot)
 {
 	sg_linalg->set_gpu_backend(new LinalgBackendViennaCL());