     return ret;
   }
 
-  TC_FORCE_INLINE __device__ TVector operator-(const TVector &o) {
+  TC_FORCE_INLINE __host__ __device__ TVector operator-(const TVector &o) {
     TVector ret;
     for (int i = 0; i < dim; i++) {
       ret[i] = d[i] - o[i];

  }

  template <int dim__ = dim_>
  TC_FORCE_INLINE __device__ TVector(T x, T y) {
    static_assert(dim__ == 2, "");
    d[0] = x;
    d[1] = y;
  }

  template <int dim__ = dim_>
  TC_FORCE_INLINE __device__ TVector(T x, T y, T z) {
    static_assert(dim__ == 3, "");
    d[0] = x;
    d[1] = y;
    d[2] = z;
  }

  TC_FORCE_INLINE __device__ TVector(T x = 0) {
    for (int i = 0; i < dim; i++) {
      d[i] = x;
    }

  real det = determinant(mat);
  return 1.0f / det *
         TMatrix<real, 3>(mat[1][1] * mat[2][2] - mat[2][1] * mat[1][2],
                          mat[2][1] * mat[0][2] - mat[0][1] * mat[2][2],
                          mat[0][1] * mat[1][2] - mat[1][1] * mat[0][2],
                          mat[2][0] * mat[1][2] - mat[1][0] * mat[2][2],
                          mat[0][0] * mat[2][2] - mat[2][0] * mat[0][2],
                          mat[1][0] * mat[0][2] - mat[0][0] * mat[1][2],
                          mat[1][0] * mat[2][1] - mat[2][0] * mat[1][1],
                          mat[2][0] * mat[0][1] - mat[0][0] * mat[2][1],
                          mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]);
}

TC_FORCE_INLINE __device__ real sgn(real x) {
