Merge pull request #2809 from lambday/develop

Added matrix-type independent elementwise-product in linalg

src/shogun/lib/GPUMatrix.h

@@ -111,23 +111,24 @@ template <class T> class CGPUMatrix
 	/** Creates a gpu matrix using data from an SGMatrix */
 	CGPUMatrix(const SGMatrix<T>& cpu_mat);
 
+	/** Converts the matrix into an SGMatrix */
+	operator SGMatrix<T>() const;
+
 #ifndef SWIG // SWIG should skip this part
 #ifdef HAVE_EIGEN3
 	CGPUMatrix(const EigenMatrixXt& cpu_mat);
 
 	/** Converts the matrix into an Eigen3 matrix */
 	operator EigenMatrixXt() const;
-#endif
-#endif
-
-	/** Converts the matrix into an SGMatrix */
-	operator SGMatrix<T>() const;
+#endif // HAVE_EIGEN3
 
 	/** Returns a ViennaCL matrix wrapped around the data of this matrix. Can be
 	 * used to call native ViennaCL methods on this matrix
 	 */
 	VCLMatrixBase vcl_matrix();
 
+#endif // SWIG
+
 	/** Sets all the elements of the matrix to zero */
 	void zero();
 

src/shogun/mathematics/linalg/internal/implementation/ElementwiseProduct.h

@@ -74,13 +74,49 @@ struct elementwise_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;
 
-	/** Eigen3 vector type */
-	typedef Eigen::Matrix<T,Eigen::Dynamic,1> VectorXt;
+	/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication.
+	 *
+	 * This version returns the result in a newly created matrix. If elementwise-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
+	 * @return The result of the operation
+	 */
+	static ReturnType compute(SGMatrix<T> A, SGMatrix<T> B)
+	{
+		REQUIRE(A.matrix, "Matrix A is not initialized!\n");
+		REQUIRE(B.matrix, "Matrix A is not initialized!\n");
+
+		REQUIRE(A.num_rows == B.num_rows && A.num_cols == B.num_cols,
+				"Dimension mismatch! A(%d x %d) vs B(%d x %d)\n",
+				A.num_rows, A.num_cols, B.num_rows, B.num_cols);
 
-	/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication */
+		ReturnType retMatrix(A.num_rows, A.num_cols);
+		compute(A, B, retMatrix);
+
+		return retMatrix;
+	}
+
+	/** Performs the operation C = A .* B where ".*" denotes elementwise 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
+	 * @param C Result of the operation
+	 */
 	static void compute(SGMatrix<T> A, SGMatrix<T> B, SGMatrix<T> C)
 	{
 		Eigen::Map<MatrixXt> A_eig = A;
@@ -101,7 +137,46 @@ struct elementwise_product<Backend::VIENNACL, Matrix>
 	/** Scalar type */
 	typedef typename Matrix::Scalar T;
 
-	/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication */
+	/** Return type */
+	typedef CGPUMatrix<T> ReturnType;
+
+	/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication.
+	 *
+	 * This version returns the result in a newly created matrix. If elementwise-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
+	 * @return The result of the operation
+	 */
+	static ReturnType compute(CGPUMatrix<T> A, CGPUMatrix<T> B)
+	{
+		REQUIRE(A.matrix, "Matrix A is not initialized!\n");
+		REQUIRE(B.matrix, "Matrix A is not initialized!\n");
+
+		REQUIRE(A.num_rows == B.num_rows && A.num_cols == B.num_cols,
+				"Dimension mismatch! A(%d x %d) vs B(%d x %d)\n",
+				A.num_rows, A.num_cols, B.num_rows, B.num_cols);
+
+		ReturnType retMatrix(A.num_rows, A.num_cols);
+		compute(A, B, retMatrix);
+
+		return retMatrix;
+	}
+
+	/** Performs the operation C = A .* B where ".*" denotes elementwise 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
+	 * @param C Result of the operation
+	 */
 	static void compute(CGPUMatrix<T> A, CGPUMatrix<T> B, CGPUMatrix<T> C)
 	{
 		C.vcl_matrix() = viennacl::linalg::element_prod(A.vcl_matrix(), B.vcl_matrix());

src/shogun/mathematics/linalg/internal/modules/Core.h

@@ -122,13 +122,40 @@ void subtract(Matrix A, Matrix B, Matrix C,
 	implementation::add<backend, Matrix>::compute(A, B, C, alpha, -1*beta);
 }
 
-/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication */
+/** Performs the operation C = A .* B where ".*" denotes elementwise 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
+ * @param C Result of the operation
+ */
 template <Backend backend=linalg_traits<Core>::backend,class Matrix>
 void elementwise_product(Matrix A, Matrix B, Matrix C)
 {
 	implementation::elementwise_product<backend, Matrix>::compute(A, B, C);
 }
 
+/** Performs the operation C = A .* B where ".*" denotes elementwise multiplication.
+ *
+ * This version returns the result in a newly created matrix. If elementwise-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
+ * @return The result of the operation
+ */
+template <Backend backend=linalg_traits<Core>::backend,class Matrix>
+typename implementation::elementwise_product<backend,Matrix>::ReturnType elementwise_product(Matrix A, Matrix B)
+{
+	return implementation::elementwise_product<backend,Matrix>::compute(A, B);
+}
+
 /**
  * Wrapper method for internal implementation of square of co-efficients that works
  * with generic dense matrices.

tests/unit/mathematics/linalg/ElementwiseProduct_unittest.cc

@@ -46,43 +46,81 @@
 using namespace shogun;
 
 #ifdef HAVE_EIGEN3
-TEST(ElementwiseProduct, eigen3_backend)
+TEST(ElementwiseProduct, SGMatrix_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::elementwise_product<linalg::Backend::EIGEN3>(A, B, C);
-
+
+	for (int32_t i=0; i<9; i++)
+		EXPECT_NEAR(A[i]*B[i], C[i], 1e-15);
+}
+
+#ifdef HAVE_VIENNACL
+TEST(ElementwiseProduct, CGPUMatrix_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::elementwise_product<linalg::Backend::EIGEN3>(A, B);
+
 	for (int32_t i=0; i<9; i++)
 		EXPECT_NEAR(A[i]*B[i], C[i], 1e-15);
 }
+#endif // HAVE_VIENNACL
 #endif // HAVE_EIGEN3
 
 #ifdef HAVE_VIENNACL
-TEST(ElementwiseProduct, viennacl_backend)
+TEST(ElementwiseProduct, CGPUMatrix_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::elementwise_product<linalg::Backend::VIENNACL>(A, B, C);
-
+
+	for (int32_t i=0; i<9; i++)
+		EXPECT_NEAR(A[i]*B[i], C[i], 1e-15);
+}
+
+#ifdef HAVE_EIGEN3
+TEST(ElementwiseProduct, 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::elementwise_product<linalg::Backend::VIENNACL>(A, B);
+
 	for (int32_t i=0; i<9; i++)
 		EXPECT_NEAR(A[i]*B[i], C[i], 1e-15);
 }
+#endif // HAVE_EIGEN3
 #endif // HAVE_VIENNACL
 
 #endif // HAVE_LINALG_LIB