add matmul and inverse

traph/include/traph/core/operation.h

@@ -80,6 +80,33 @@ namespace traph
 		}
 	};
 
+	class MatmulOp : public OpBase
+	{
+	public:
+		virtual TensorInterfacePtr forward(std::vector<TensorInterfacePtr> inputs) override
+		{
+			assert(inputs.size() == 2);
+
+			TensorInterfacePtr left_input = inputs[0];
+			TensorInterfacePtr right_input = inputs[1];
+			TensorInterfacePtr result = left_input->matmul(right_input);
+
+			context.save(left_input);
+			context.save(right_input);
+
+			return result;
+		}
+
+		virtual std::vector<TensorBasePtr<f32>> backward(TensorBasePtr<f32> output_grad) override
+		{
+			auto saved_tensors = context.get_saved_tensors();
+			assert(saved_tensors.size() == 1);
+			std::shared_ptr<TensorBase<f32>> left_out = std::dynamic_pointer_cast<TensorBase<f32>>(output_grad->matmul(saved_tensors[0]->inverse()));
+			std::shared_ptr<TensorBase<f32>> right_out = std::dynamic_pointer_cast<TensorBase<f32>>(saved_tensors[0]->inverse()->matmul(output_grad));
+			return { left_out, right_out };
+		}
+	};
+
 	class SelectOp : public OpBase
 	{
 	public:

traph/include/traph/core/tensor.h

@@ -33,7 +33,8 @@ namespace traph
         virtual void cos_() = 0;
         virtual std::shared_ptr<TensorBase<f32>> create_grad() = 0;
         virtual device_id device() = 0;
-        virtual std::shared_ptr<TensorInterface> matmul() const = 0;
+        virtual std::shared_ptr<TensorInterface> inverse() const = 0;
+        virtual std::shared_ptr<TensorInterface> matmul(std::shared_ptr<TensorInterface> mat) const = 0;
         virtual idx_type offset() const = 0;
 		virtual layout_type order() const = 0;
         virtual platform_type platform() = 0;
@@ -77,8 +78,9 @@ namespace traph
         virtual const T* data_ptr() const = 0;
         virtual device_id device() = 0;
         virtual void fill_(T value) = 0;
+        virtual std::shared_ptr<TensorInterface> inverse() const = 0;
         virtual T item() const = 0;
-        virtual std::shared_ptr<TensorInterface> matmul() const = 0;
+        virtual std::shared_ptr<TensorInterface> matmul(std::shared_ptr<TensorInterface> mat) const = 0;
         virtual idx_type offset() const = 0;
 		virtual layout_type order() const = 0;
         virtual platform_type platform() = 0;

traph/include/traph/tensor/arithmetic.h

@@ -3,6 +3,7 @@
 
 #include <utility>
 #include <cmath>
+#include <memory>
 
 #include <traph/core/type.h>
 #include <traph/core/index.h>
@@ -11,78 +12,11 @@
 
 namespace traph
 {
-	template<class T>
-	Tensor<T> abs(const Tensor<T> &t)
-	{
-		Tensor<T> result(t.size());
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = t.offset(); i < flat_size_end; ++i)
-		{
-			result.data()[i] = std::abs(t.data()[i]);
-		}
-
-		return result;
-	}
-
-	template<class T>
-	Tensor<T> acos(const Tensor<T> &t)
-	{
-		Tensor<T> result(t.size());
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = t.offset(); i < flat_size_end; ++i)
-		{
-			result.data()[i] = std::acos(t.data()[i]);
-		}
+	template<typename T>
+	class Tensor;
 
-		return result;
-	}
-
-	// add fallback
 	template<class T>
-	Tensor<T> add(const Tensor<T> &t, T v)
-	{
-		Tensor<T> result(t.size());
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = 0; i < flat_size_end; ++i)
-		{
-			result.data()[i] = t.data()[i] + v;
-		}
-
-		return result;
-	}
-
-	Tensor<f32> add(const Tensor<f32> &t, f32 v);
-
-	Tensor<f64> add(const Tensor<f64> &t, f64 v);
-
-	template<class T>
-	Tensor<T> asin(const Tensor<T> &t)
-	{
-		Tensor<T> result(t.size());
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = t.offset(); i < flat_size_end; ++i)
-		{
-			result.data()[i] = std::asin(t.data()[i]);
-		}
-
-		return result;
-	}
-
-	template<class T>
-	Tensor<T> atan(const Tensor<T> &t)
-	{
-		Tensor<T> result(t.size());
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = t.offset(); i < flat_size_end; ++i)
-		{
-			result.data()[i] = std::atan(t.data()[i]);
-		}
-
-		return result;
-	}
-
-	template<class T>
-	void matmul_check(const Tensor<T> &a, const Tensor<T> &b)
+	void matmul_check(const Tensor<T>& a, const Tensor<T>& b)
 	{
 		// check dimension
 		if (a.size().size() > 2 || b.size().size() > 2)
@@ -95,115 +29,29 @@ namespace traph
 			throw std::runtime_error("matmul: Dimension 0 of the first matrix shall be equal to dimension 1 of the second matrix.");
 		}
 		// check order
-		if (a.layout() != b.layout())
+		if (a.order() != b.order())
 		{
 			throw std::runtime_error("matmul: Two matrix shall have the same layout.");
 		}
 	}
 
-	// matmull fallback
-	template<class T>
-	Tensor<T> matmul(const Tensor<T> &a, const Tensor<T> &b)
-	{
-		// check
-		matmul_check(a, b);
-		// result
-		Tensor<T> result = zeros<T>({ a.size()[0], b.size()[1] });
-		// compute
-		idx_type m = a.size()[0];
-		idx_type n = b.size()[1];
-		idx_type k = a.size()[1];
-
-		for (idx_type j = 0; j < n; ++j)
-		{
-			idx_type jm = j * m;
-			idx_type jk = j * k;
-			for (idx_type i = 0; i < m; ++i)
-			{
-				idx_type jm_i = jm + i;
-				for (idx_type p = 0; p < k; ++p)
-				{
-					result.data()[jm_i] = a.data()[p*m + i] + b.data()[jk + p];
-				}
-			}
-		}
-
-		return result;
-	}
-
-	Tensor<u8> matmul(const Tensor<u8> &a, const Tensor<u8> &b);
-
-	Tensor<i8> matmul(const Tensor<i8> &a, const Tensor<i8> &b);
-
-	Tensor<i16> matmul(const Tensor<i16> &a, const Tensor<i16> &b);
-
-	Tensor<i32> matmul(const Tensor<i32> &a, const Tensor<i32> &b);
-
-	Tensor<i64> matmul(const Tensor<i64> &a, const Tensor<i64> &b);
-
-	Tensor<f32> matmul(const Tensor<f32> &a, const Tensor<f32> &b);
-
-	Tensor<f64> matmul(const Tensor<f64> &a, const Tensor<f64> &b);
-
-	template<class T>
-	Tensor<T> mean(const Tensor<T> &input)
-	{
-		T result{};
-		idx_type flat_size = input.size().flat_size();
-		idx_type offset = input.offset();
-		const T *input_data = input.data();
-		for (idx_type i = offset; i < offset + flat_size; ++i)
-		{
-			result += input_data[i];
-		}
+	std::shared_ptr<Tensor<u8>> matmul_impl(const Tensor<u8>& a, const Tensor<u8>& b);
 
-		return result / flat_size;
-	}
+	std::shared_ptr<Tensor<i8>> matmul_impl(const Tensor<i8>& a, const Tensor<i8>& b);
 
-	template<class T>
-	Tensor<T> mul(const Tensor<T> &input, T value)
-	{
-		Tensor<T> result(input.size());
-		idx_type flat_size = input.size().flat_size();
-		idx_type offset = input.offset();
-		const T *input_data = input.data();
-		T *result_data = result.data();
-		for (idx_type i = 0; i < flat_size; ++i)
-		{
-			result_data[i] += input_data[i + offset] + value;
-		}
+	std::shared_ptr<Tensor<i16>> matmul_impl(const Tensor<i16>& a, const Tensor<i16>& b);
 
-		return result;
-	}
+	std::shared_ptr<Tensor<i32>> matmul_impl(const Tensor<i32>& a, const Tensor<i32>& b);
 
-	template<class T>
-	void mul_check(const Tensor<T> &input, const Tensor<T> & other)
-	{
-		if(!strict_same_shape(input, other))
-			throw std::runtime_error("mul: Two tensor must have the same shape or be broadcastable.");
-	}
+	std::shared_ptr<Tensor<i64>> matmul_impl(const Tensor<i64>& a, const Tensor<i64>& b);
 
-	template<class T>
-	Tensor<T> mul(const Tensor<T> &input, const Tensor<T> & other)
-	{
-		// check
-		mul_check(input, other);
+	std::shared_ptr<Tensor<f32>> matmul_impl(const Tensor<f32>& a, const Tensor<f32>& b);
 
-		Tensor<T> result(input.size());
-		idx_type flat_size = input.size().flat_size();
-		idx_type input_offset = input.offset();
-		idx_type other_offset = other.offset();
-		const T *input_data = input.data();
-		const T *other_data = other.data();
-		T *result_data = result.data();
+	std::shared_ptr<Tensor<f64>> matmul_impl(const Tensor<f64>& a, const Tensor<f64>& b);
 
-		for (idx_type i = 0; i < flat_size; ++i)
-		{
-			result_data[i] += input_data[i + input_offset] + other_data[i + other_offset];
-		}
+	std::shared_ptr<Tensor<f32>> inverse_impl(const Tensor<f32>& a);
 
-		return result;
-	}
+	std::shared_ptr<Tensor<f64>> inverse_impl(const Tensor<f64>& a);
 }
 
 #endif

traph/include/traph/tensor/tensor.h

@@ -3,6 +3,7 @@
 
 #include <initializer_list>
 #include <cmath>
+#include <cfenv>
 #include <memory>
 #include <functional>
 #include <stdexcept>
@@ -16,6 +17,7 @@
 #include<traph/core/tensor.h>
 
 #include<traph/tensor/tensor_storage.h>
+#include<traph/tensor/arithmetic.h>
 
 namespace traph
 {
@@ -75,8 +77,9 @@ namespace traph
         virtual const T* data_ptr() const override;
         virtual device_id device() override;
         virtual void fill_(T value) override;
+        virtual std::shared_ptr<TensorInterface> inverse() const override;
         virtual T item() const override;
-        virtual std::shared_ptr<TensorInterface> matmul() const override;
+        virtual std::shared_ptr<TensorInterface> matmul(std::shared_ptr<TensorInterface> mat) const override;
 		virtual idx_type offset() const override;
 		virtual layout_type order() const override;
         virtual platform_type platform() override;
@@ -354,6 +357,13 @@ namespace traph
     {
         apply_([&value](T a)->T {return value; });
     }
+
+    template<typename T>
+    std::shared_ptr<TensorInterface> Tensor<T>::inverse() const
+    {
+        return std::dynamic_pointer_cast<TensorInterface>(inverse_impl(*this);
+    }
+
     template<typename T>
     T Tensor<T>::item() const
     {
@@ -367,9 +377,10 @@ namespace traph
         }
     }
     template<typename T>
-    std::shared_ptr<TensorInterface> Tensor<T>::matmul() const
+    std::shared_ptr<TensorInterface> Tensor<T>::matmul(std::shared_ptr<TensorInterface> mat) const
     {
-
+		auto right_matrix = std::dynamic_pointer_cast<Tensor<T>>(mat);
+		return matmul_impl(*this, *right_matrix);
     }
     template<typename T>
     idx_type Tensor<T>::offset() const { return _offset; }
@@ -414,11 +425,12 @@ namespace traph
 
         // dimension
         DimVector dim;
+		std::fesetround(FE_TONEAREST);
         for(idx_type i = 0; i<slice.size(); ++i)
         {
 			auto& each = slice[i];
             dim.push_back(
-				std::ceil((each.end.value_or(_dimensions[i]) - each.start.value_or(0))/(float)each.step.value_or(1))
+				std::lrint(std::ceil((each.end.value_or(_dimensions[i]) - each.start.value_or(0))/(float)each.step.value_or(1)))
 			);
         }
         result->_dimensions = dim;

traph/source/tensor/CMakeLists.txt

@@ -7,6 +7,8 @@ SET(SOURCE_PATH ${TRAPH_PATH_SOURCE}/${LIB_NAME})
 SET(TENSOR_LIST
 	${HEADER_PATH}/tensor.h
 	${SOURCE_PATH}/tensor.cpp
+	${HEADER_PATH}/arithmetic.h
+	${SOURCE_PATH}/arithmetic.cpp
 )
 
 ADD_LIBRARY(${LIB_OUTNAME} ${TENSOR_LIST})

traph/source/tensor/arithmetic.cpp

@@ -2,6 +2,7 @@
 #include <stdexcept>
 #include <algorithm>
 
+#include <traph/tensor/tensor.h>
 #include <traph/tensor/arithmetic.h>
 
 #include <eigen3/Eigen/Dense>
@@ -17,142 +18,108 @@
 
 namespace traph
 {
-	/*
-	Tensor<f32> add(const Tensor<f32> &t, f32 v)
-	{
-		Tensor<f32> result(t.size());
-#ifdef TRAPH_BUILD_EIGEN
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = 0; i < flat_size_end; ++i)
-		{
-			result.data()[i] = t.data()[i] + v;
-		}
-#elif defined TRAPH_BUILD_MKL
-		result.fill_(v);
-		cblas_saxpy(t.size().flat_size(), 1.f, t.data(), 1, result.data(), 1);
-#endif
-		return result;
-	}
-
-	Tensor<f64> add(const Tensor<f64> &t, f64 v)
-	{
-		Tensor<f64> result(t.size());
-#ifdef TRAPH_BUILD_EIGEN
-		idx_type flat_size_end = t.size().flat_size();
-		for (idx_type i = 0; i < flat_size_end; ++i)
-		{
-			result.data()[i] = t.data()[i] + v;
-		}
-#elif defined TRAPH_BUILD_MKL
-		result.fill_(v);
-		cblas_daxpy(t.size().flat_size(), 1.f, t.data(), 1, result.data(), 1);
-#endif
-		return result;
-	}
-	*/
-
-	Tensor<u8> matmul(const Tensor<u8> &a, const Tensor<u8> &b)
+	std::shared_ptr<Tensor<u8>> matmul_impl(const Tensor<u8>& a, const Tensor<u8>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<u8> result = zeros<u8>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<u8>> result(new Tensor<u8>(a.size()[0], b.size()[1]));
 
 		// copy data
-		Eigen::Map<const Eigen::Matrix<u8, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<u8, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<u8, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<u8, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<u8, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 		return result;
 	}
 
-	Tensor<i8> matmul(const Tensor<i8> &a, const Tensor<i8> &b)
+	std::shared_ptr<Tensor<i8>> matmul_impl(const Tensor<i8>& a, const Tensor<i8>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<i8> result = zeros<i8>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<i8>> result(new Tensor<i8>(a.size()[0], b.size()[1]));
 
 		// copy data
-		Eigen::Map<const Eigen::Matrix<i8, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<i8, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i8, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i8, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<i8, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 		return result;
 	}
 
-	Tensor<i16> matmul(const Tensor<i16> &a, const Tensor<i16> &b)
+	std::shared_ptr<Tensor<i16>> matmul_impl(const Tensor<i16>& a, const Tensor<i16>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<i16> result = zeros<i16>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<i16>> result(new Tensor<i16>(a.size()[0], b.size()[1]));
 
 		// copy data
-		Eigen::Map<const Eigen::Matrix<i16, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<i16, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i16, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i16, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<i16, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 		return result;
 	}
 
-	Tensor<i32> matmul(const Tensor<i32> &a, const Tensor<i32> &b)
+	std::shared_ptr<Tensor<i32>> matmul_impl(const Tensor<i32>& a, const Tensor<i32>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<i32> result = zeros<i32>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<i32>> result(new Tensor<i32>(a.size()[0], b.size()[1]));
 
 		// copy data
-		Eigen::Map<const Eigen::Matrix<i32, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<i32, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i32, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i32, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<i32, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 		return result;
 	}
 
-	Tensor<i64> matmul(const Tensor<i64> &a, const Tensor<i64> &b)
+	std::shared_ptr<Tensor<i64>> matmul_impl(const Tensor<i64>& a, const Tensor<i64>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<i64> result = zeros<i64>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<i64>> result(new Tensor<i64>(a.size()[0], b.size()[1]));
 
 		// copy data
-		Eigen::Map<const Eigen::Matrix<i64, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<i64, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i64, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<i64, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<i64, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 		return result;
 	}
 
-	Tensor<f32> matmul(const Tensor<f32> &a, const Tensor<f32> &b)
+	std::shared_ptr<Tensor<f32>> matmul_impl(const Tensor<f32>& a, const Tensor<f32>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<f32> result = zeros<f32>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<f32>> result(new Tensor<f32>(a.size()[0], b.size()[1]));
 
 #ifdef TRAPH_BUILD_EIGEN
 		// copy data
-		Eigen::Map<const Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result->data_ptr());
 #elif defined TRAPH_BUILD_MKL
-		CBLAS_LAYOUT a_layout = a.layout() == layout_type::column_major ? CBLAS_LAYOUT::CblasColMajor : CBLAS_LAYOUT::CblasRowMajor;
+		CBLAS_LAYOUT a_layout = a.order() == layout_type::column_major ? CBLAS_LAYOUT::CblasColMajor : CBLAS_LAYOUT::CblasRowMajor;
 
 		cblas_sgemm(a_layout,
 			CBLAS_TRANSPOSE::CblasNoTrans,
@@ -161,34 +128,34 @@ namespace traph
 			b.size()[1],
 			a.size()[1],
 			1.f,
-			a.data(),
+			a.data_ptr(),
 			a.size()[0],
-			b.data(),
+			b.data_ptr(),
 			b.size()[0],
 			0.f,
-			result.data(),
-			result.size()[0]);
+			result->data_ptr(),
+			result->size()[0]);
 #endif
 		return result;
 	}
 
-	Tensor<f64> matmul(const Tensor<f64> &a, const Tensor<f64> &b)
+	std::shared_ptr<Tensor<f64>> matmul_impl(const Tensor<f64>& a, const Tensor<f64>& b)
 	{
 		// check
 		matmul_check(a, b);
 		// result
-		Tensor<f64> result = zeros<f64>({ a.size()[0], b.size()[1] });
+		std::shared_ptr<Tensor<f64>> result(new Tensor<f64>(a.size()[0], b.size()[1]));
 
 #ifdef TRAPH_BUILD_EIGEN
 		// copy data
-		Eigen::Map<const Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data() + a.offset(), a.size()[0], a.size()[1]);
-		Eigen::Map<const Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data() + b.offset(), b.size()[0], b.size()[1]);
+		Eigen::Map<const Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+		Eigen::Map<const Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic>> eigen_b(b.data_ptr() + b.offset(), b.size()[0], b.size()[1]);
 
 		Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a * eigen_b;
 		// copy to result
-		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data());
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * b.size()[1], result.data_ptr());
 #elif defined TRAPH_BUILD_MKL
-		CBLAS_LAYOUT a_layout = a.layout() == layout_type::column_major ? CBLAS_LAYOUT::CblasColMajor : CBLAS_LAYOUT::CblasRowMajor;
+		CBLAS_LAYOUT a_layout = a.order() == layout_type::column_major ? CBLAS_LAYOUT::CblasColMajor : CBLAS_LAYOUT::CblasRowMajor;
 
 		cblas_dgemm(a_layout,
 			CBLAS_TRANSPOSE::CblasNoTrans,
@@ -197,15 +164,43 @@ namespace traph
 			b.size()[1],
 			a.size()[1],
 			1.f,
-			a.data(),
+			a.data_ptr(),
 			a.size()[0],
-			b.data(),
+			b.data_ptr(),
 			b.size()[0],
 			0.f,
-			result.data(),
-			result.size()[0]);
+			result->data_ptr(),
+			result->size()[0]);
 #endif
 
 		return result;
 	}
+
+	std::shared_ptr<Tensor<f32>> inverse_impl(const Tensor<f32>& a)
+	{
+		// result
+		std::shared_ptr<Tensor<f32>> result(new Tensor<f32>(a.size()[0], a.size()[1]));
+
+		// copy data
+		Eigen::Map<const Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+
+		Eigen::Matrix<f32, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a.inverse();
+		// copy to result
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * a.size()[1], result->data_ptr());
+		return result;
+	}
+
+	std::shared_ptr<Tensor<f64>> inverse_impl(const Tensor<f64>& a)
+	{
+		// result
+		std::shared_ptr<Tensor<f64>> result(new Tensor<f64>(a.size()[0], a.size()[1]));
+
+		// copy data
+		Eigen::Map<const Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic>> eigen_a(a.data_ptr() + a.offset(), a.size()[0], a.size()[1]);
+
+		Eigen::Matrix<f64, Eigen::Dynamic, Eigen::Dynamic> eigen_c = eigen_a.inverse();
+		// copy to result
+		std::copy(eigen_c.data(), eigen_c.data() + a.size()[0] * a.size()[1], result->data_ptr());
+		return result;
+	}
 }