Repository: sandialabs/pyttb

Issue:
Display tensors by frontal slice rather than lateral slice
The default display in pyttb is to display the 2D slices w.r.t. the last 2 modes. This is in contrast to MATLAB which displays them w.r.t. the first two modes. For a 3-way tensor, pyttb displays the lateral slices rather than the frontal slices. I think it should be changed or at least there should be a function that displays it the MATLAB TTB way. 

Example in MATLAB:
```
X = tensor([1:8],[2 2 2])
```
Output:
```
X is a tensor of size 2 x 2 x 2
	X(:,:,1) = 
	     1     3
	     2     4
	X(:,:,2) = 
	     5     7
	     6     8
```

Same example in Python:
```
Y = ttb.tensor(np.arange(1,9), (2, 2, 2));
Y
```
Output:
```
tensor of shape (2, 2, 2)
data[0, :, :] =
[[1 5]
 [3 7]]
data[1, :, :] =
[[2 6]
 [4 8]]
```

I expect these display orders are artifacts of the underlying default storage. By default, MATLAB prints in such a way that it's linear in memory. Python uses the opposite ordering by default, so it's perhaps not surprising that its printing is in the opposite order. But since pyttb is storing things in the MATLAB order internally, it should also print them this way.

Discussion:
@ntjohnson1: I think this was just a bug in translated the MATLAB. The indices and slices are consistent but we weren't actually printing it in a way that was compatible with MATLAB Tensortoolbox or with the underlying numpy printout. @ntjohnson1: ```python
>>>import numpy as np; import pyttb as ttb
>>>a = np.arange(16)
>>>b = a.reshape((2,2,2,2), order="F")
>>>b
array([[[[ 0,  8],
         [ 4, 12]],
        [[ 2, 10],
         [ 6, 14]]],
       [[[ 1,  9],
         [ 5, 13]],
        [[ 3, 11],
         [ 7, 15]]]])
>>>c = ttb.tensor(b)
>>>c
tensor of shape (2, 2, 2, 2)
data[0, 0, :, :] =
[[ 0  8]
 [ 4 12]]
data[1, 0, :, :] =
[[ 1  9]
 [ 5 13]]
data[0, 1, :, :] =
[[ 2 10]
 [ 6 14]]
data[1, 1, :, :] =
[[ 3 11]
 [ 7 15]]
```

Solution patch:
diff --git a/pyttb/cp_als.py b/pyttb/cp_als.py
index 86d5585..98bc801 100644
--- a/pyttb/cp_als.py
+++ b/pyttb/cp_als.py
@@ -92,7 +92,7 @@ def cp_als(  # noqa: PLR0912,PLR0913,PLR0915
      Iter 1: f = ... f-delta = ...
      Final f = ...
     >>> print(M) # doctest: +ELLIPSIS
-    ktensor of shape (2, 2)
+    ktensor of shape (2, 2) with order F
     weights=[108.4715... 8.6114...]
     factor_matrices[0] =
     [[0.4187... 0.3989...]
@@ -101,7 +101,7 @@ def cp_als(  # noqa: PLR0912,PLR0913,PLR0915
     [[0.6188... 0.2581...]
      [0.7854... 0.9661...]]
     >>> print(Minit) # doctest: +ELLIPSIS
-    ktensor of shape (2, 2)
+    ktensor of shape (2, 2) with order F
     weights=[1. 1.]
     factor_matrices[0] =
     [[4.1702...e-01 7.2032...e-01]
@@ -185,7 +185,7 @@ def cp_als(  # noqa: PLR0912,PLR0913,PLR0915
 
     # Reduce dimorder to only those modes we will optimize
     dimorder_in = dimorder  # save for output
-    dimorder = [d for d in dimorder if d in optdims]
+    dimorder = [int(d) for d in dimorder if d in optdims]
 
     # Store the last MTTKRP result to accelerate fitness computation
     U_mttkrp = np.zeros((input_tensor.shape[dimorder[-1]], rank))
diff --git a/pyttb/cp_apr.py b/pyttb/cp_apr.py
index 293b60a..c5b1e1d 100644
--- a/pyttb/cp_apr.py
+++ b/pyttb/cp_apr.py
@@ -917,8 +917,8 @@ def tt_cp_apr_pqnr(  # noqa: PLR0912,PLR0913,PLR0915
                 delg = np.zeros((rank, lbfgsMem))
                 rho = np.zeros((lbfgsMem,))
                 lbfgsPos = 0
-                m_rowOLD = np.empty(())
-                gradOLD = np.empty(())
+                m_rowOLD = np.empty((), dtype=m_row.dtype)
+                gradOLD = np.empty((), dtype=m_row.dtype)
 
                 # Iteratively solve the row subproblem with projected quasi-Newton steps
                 for i in range(maxinneriters):
@@ -983,8 +983,8 @@ def tt_cp_apr_pqnr(  # noqa: PLR0912,PLR0913,PLR0915
                     isRowNOTconverged[jj] = 1
 
                     # Update the L-BFGS approximation.
-                    tmp_delm = m_row - m_rowOLD
-                    tmp_delg = gradM - gradOLD
+                    tmp_delm: np.ndarray = m_row - m_rowOLD
+                    tmp_delg: np.ndarray = gradM - gradOLD
                     tmp_delm_dot = tmp_delm.dot(tmp_delg.transpose())
                     if not np.any(tmp_delm_dot == 0):
                         tmp_rho = 1 / tmp_delm_dot
diff --git a/pyttb/gcp/fg_est.py b/pyttb/gcp/fg_est.py
index 32ccd00..7dc7a51 100644
--- a/pyttb/gcp/fg_est.py
+++ b/pyttb/gcp/fg_est.py
@@ -162,7 +162,7 @@ def estimate_helper(
     ndim = subs.shape[1]
 
     # Create exploded U's from the model factor matrices
-    Uexp = [np.empty(())] * ndim
+    Uexp = [np.empty((), dtype=factors[0].dtype)] * ndim
     for k in range(ndim):
         Uexp[k] = factors[k][subs[:, k], :]
 
diff --git a/pyttb/gcp/samplers.py b/pyttb/gcp/samplers.py
index 8af8e05..bd9e29a 100644
--- a/pyttb/gcp/samplers.py
+++ b/pyttb/gcp/samplers.py
@@ -272,7 +272,7 @@ def nonzeros(
 
     # Select nonzeros
     if samples == nnz:
-        nidx = np.arange(0, nnz)
+        nidx: np.ndarray = np.arange(0, nnz, dtype=int)
     elif with_replacement or samples < nnz:
         nidx = np.random.choice(nnz, size=samples, replace=with_replacement)
     else:
diff --git a/pyttb/hosvd.py b/pyttb/hosvd.py
index 4243e22..7936e33 100644
--- a/pyttb/hosvd.py
+++ b/pyttb/hosvd.py
@@ -129,7 +129,7 @@ def hosvd(  # noqa: PLR0912,PLR0913,PLR0915
 
         # Shrink!
         if sequential:
-            Y = Y.ttm(factor_matrices[k].transpose(), k)
+            Y = Y.ttm(factor_matrices[k].transpose(), int(k))
     # Extract final core
     if sequential:
         G = Y
diff --git a/pyttb/ktensor.py b/pyttb/ktensor.py
index 94b52a7..e32dfd6 100644
--- a/pyttb/ktensor.py
+++ b/pyttb/ktensor.py
@@ -102,7 +102,7 @@ class ktensor:
 
         >>> K = ttb.ktensor()
         >>> print(K)
-        ktensor of shape ()
+        ktensor of shape () with order F
         weights=[]
         factor_matrices=[]
 
@@ -114,7 +114,7 @@ class ktensor:
         >>> fm1 = np.array([[5.0, 6.0], [7.0, 8.0]])
         >>> K = ttb.ktensor([fm0, fm1], weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -131,7 +131,7 @@ class ktensor:
         >>> factor_matrices = [fm0, fm1]
         >>> K = ttb.ktensor([fm0, fm1])
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[1. 2.]
@@ -235,7 +235,7 @@ class ktensor:
         >>> np.random.seed(1)
         >>> K = ttb.ktensor.from_function(np.random.random_sample, (2, 3, 4), 2)
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[4.1702...e-01 7.2032...e-01]
@@ -254,7 +254,7 @@ class ktensor:
 
         >>> K = ttb.ktensor.from_function(np.ones, (2, 3, 4), 2)
         >>> print(K)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[1. 1.]
@@ -273,7 +273,7 @@ class ktensor:
 
         >>> K = ttb.ktensor.from_function(np.zeros, (2, 3, 4), 2)
         >>> print(K)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[0. 0.]
@@ -334,7 +334,7 @@ class ktensor:
         >>> data = np.arange(1, rank * sum(shape) + 1).astype(float)
         >>> K = ttb.ktensor.from_vector(data[:], shape, False)
         >>> print(K)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[1. 3.]
@@ -356,7 +356,7 @@ class ktensor:
         >>> weights_and_data = np.concatenate((weights, data), axis=0)
         >>> K = ttb.ktensor.from_vector(weights_and_data[:], shape, True)
         >>> print(K)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[2. 2.]
         factor_matrices[0] =
         [[1. 3.]
@@ -412,6 +412,11 @@ class ktensor:
 
         return cls(factor_matrices, weights, copy=False)
 
+    @property
+    def order(self) -> Literal["F"]:
+        """Return the data layout of the underlying storage."""
+        return "F"
+
     def arrange(
         self,
         weight_factor: Optional[int] = None,
@@ -449,7 +454,7 @@ class ktensor:
         >>> fm1 = np.array([[5.0, 6.0], [7.0, 8.0]])
         >>> K = ttb.ktensor([fm0, fm1], weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -463,7 +468,7 @@ class ktensor:
         >>> p = [1, 0]
         >>> K.arrange(permutation=p)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[2. 1.]
         factor_matrices[0] =
         [[2. 1.]
@@ -477,7 +482,7 @@ class ktensor:
 
         >>> K.arrange()
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[89.4427... 27.2029...]
         factor_matrices[0] =
         [[0.4472... 0.3162...]
@@ -490,7 +495,7 @@ class ktensor:
 
         >>> K.arrange(weight_factor=1)
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[0.4472... 0.3162...]
@@ -551,7 +556,7 @@ class ktensor:
         >>> np.random.seed(1)
         >>> K = ttb.ktensor.from_function(np.random.random_sample, (2, 3, 4), 2)
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[4.1702...e-01 7.2032...e-01]
@@ -571,7 +576,7 @@ class ktensor:
         >>> K2 = K.copy()
         >>> K2.weights = np.array([2.0, 3.0])
         >>> print(K2)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[2. 3.]
         factor_matrices[0] =
         [[4.1702...e-01 7.2032...e-01]
@@ -589,7 +594,7 @@ class ktensor:
         Show that the original :class:`pyttb.ktensor` is unchanged:
 
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[4.1702...e-01 7.2032...e-01]
@@ -660,7 +665,7 @@ class ktensor:
         >>> fm1 = np.array([[5.0, 6.0], [7.0, 8.0]])
         >>> K = ttb.ktensor([fm0, fm1], weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -673,7 +678,7 @@ class ktensor:
         component from each factor of the original :class:`pyttb.ktensor`:
 
         >>> K.extract([1])
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[2.]
         factor_matrices[0] =
         [[2.]
@@ -747,7 +752,7 @@ class ktensor:
         >>> K.factor_matrices[0][1, 1] = -K.factor_matrices[0][1, 1]
         >>> K.factor_matrices[1][1, 1] = -K.factor_matrices[1][1, 1]
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[ 1.  2.]
@@ -759,7 +764,7 @@ class ktensor:
         Fix the signs of the largest magnitude entries:
 
         >>> print(K.fixsigns())
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[ 1. -2.]
@@ -776,7 +781,7 @@ class ktensor:
         >>> K2.factor_matrices[1][1, 1] = -K2.factor_matrices[1][1, 1]
         >>> K = K.fixsigns(K2)
         >>> print(K)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[27.2029... 89.4427...]
         factor_matrices[0] =
         [[ 0.3162... -0.4472...]
@@ -883,7 +888,7 @@ class ktensor:
         >>> fm1 = np.array([[5.0, 6.0], [7.0, 8.0]])
         >>> K = ttb.ktensor([fm0, fm1], weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -892,7 +897,7 @@ class ktensor:
         [[5. 6.]
          [7. 8.]]
         >>> print(K.full())  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[29. 39.]
          [63. 85.]]
@@ -965,7 +970,7 @@ class ktensor:
         >>> fm1 = np.array([[5.0, 6.0], [7.0, 8.0]])
         >>> K = ttb.ktensor([fm0, fm1], weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -974,12 +979,12 @@ class ktensor:
         [[5. 6.]
          [7. 8.]]
         >>> K.full()  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[29. 39.]
          [63. 85.]]
         >>> K.to_tenmat(np.array([0]))  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -1187,7 +1192,9 @@ class ktensor:
             vals = vals + tmpvals
         return vals
 
-    def mttkrp(self, U: Union[ktensor, Sequence[np.ndarray]], n: int) -> np.ndarray:
+    def mttkrp(
+        self, U: Union[ktensor, Sequence[np.ndarray]], n: Union[int, np.integer]
+    ) -> np.ndarray:
         """
         Matricized tensor times Khatri-Rao product for :class:`pyttb.ktensor`.
 
@@ -1311,7 +1318,7 @@ class ktensor:
         --------
         >>> K = ttb.ktensor.from_function(np.ones, (2, 3, 4), 2)
         >>> print(K.normalize())  # doctest: +ELLIPSIS
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[4.898... 4.898...]
         factor_matrices[0] =
         [[0.7071... 0.7071...]
@@ -1475,7 +1482,7 @@ class ktensor:
         >>> factor_matrices = [fm0, fm1]
         >>> K = ttb.ktensor(factor_matrices, weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -1488,7 +1495,7 @@ class ktensor:
 
         >>> K1 = K.permute(np.array([1, 0]))
         >>> print(K1)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[5. 6.]
@@ -1528,7 +1535,7 @@ class ktensor:
         >>> factor_matrices = [fm0, fm1]
         >>> K = ttb.ktensor(factor_matrices, weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -1540,7 +1547,7 @@ class ktensor:
         Distribute weights of that :class:`pyttb.ktensor` to mode 0:
 
         >>> K.redistribute(0)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[1. 4.]
@@ -1709,7 +1716,7 @@ class ktensor:
             best_perm = -1 * np.ones((RA), dtype=int)
             best_score = 0.0
             for _ in range(RB):
-                idx = np.argmax(C.reshape(prod(C.shape), order="F"))
+                idx = np.argmax(C.reshape(prod(C.shape), order=self.order))
                 ij = tt_ind2sub((RA, RB), np.array(idx))
                 best_score = best_score + C[ij[0], ij[1]]
                 C[ij[0], :] = -10
@@ -1748,7 +1755,7 @@ class ktensor:
         >>> factor_matrices = [fm0, fm1]
         >>> K = ttb.ktensor(factor_matrices, weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -1762,7 +1769,7 @@ class ktensor:
 
         >>> K1 = K.symmetrize()
         >>> print(K1)  # doctest: +ELLIPSIS
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[2.3404... 4.9519...]
@@ -1828,7 +1835,7 @@ class ktensor:
         >>> factor_matrices = [fm0, fm1]
         >>> K = ttb.ktensor(factor_matrices, weights)
         >>> print(K)
-        ktensor of shape (2, 2)
+        ktensor of shape (2, 2) with order F
         weights=[1. 2.]
         factor_matrices[0] =
         [[1. 2.]
@@ -1906,7 +1913,7 @@ class ktensor:
         >>> weights_and_data = np.concatenate((weights, data), axis=0)
         >>> K = ttb.ktensor.from_vector(weights_and_data[:], shape, True)
         >>> print(K)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[2. 2.]
         factor_matrices[0] =
         [[1. 3.]
@@ -1926,7 +1933,7 @@ class ktensor:
 
         >>> K2 = ttb.ktensor.from_vector(K.tovec(), shape, True)
         >>> print(K2)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[2. 2.]
         factor_matrices[0] =
         [[1. 3.]
@@ -2019,7 +2026,7 @@ class ktensor:
         >>> K = ttb.ktensor.from_vector(weights_and_data[:], shape, True)
         >>> K0 = K.ttv(np.array([1, 1, 1]), dims=1)  # compute along a single dimension
         >>> print(K0)
-        ktensor of shape (2, 4)
+        ktensor of shape (2, 4) with order F
         weights=[36. 54.]
         factor_matrices[0] =
         [[1. 3.]
@@ -2045,7 +2052,7 @@ class ktensor:
 
         >>> K2 = K.ttv([vec4, vec3], np.array([2, 1]))
         >>> print(K2)
-        ktensor of shape (2,)
+        ktensor of shape (2,) with order F
         weights=[1800. 3564.]
         factor_matrices[0] =
         [[1. 3.]
@@ -2131,7 +2138,7 @@ class ktensor:
         >>> K1 = K.copy()
         >>> K1 = K1.update(0, vec0)
         >>> print(K1)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[2. 2.]
@@ -2152,7 +2159,7 @@ class ktensor:
         >>> vec_all = np.concatenate((vec0, vec1, vec2))
         >>> K2 = K2.update([0, 1, 2], vec_all)
         >>> print(K2)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[2. 2.]
@@ -2173,7 +2180,7 @@ class ktensor:
         >>> vec_some = np.concatenate((vec0, vec2))
         >>> K3 = K3.update([0, 2], vec_some)
         >>> print(K3)
-        ktensor of shape (2, 3, 4)
+        ktensor of shape (2, 3, 4) with order F
         weights=[1. 1.]
         factor_matrices[0] =
         [[2. 2.]
@@ -2211,7 +2218,7 @@ class ktensor:
                 self.factor_matrices[k] = np.reshape(
                     data[loc:endloc].copy(),
                     (self.shape[k], self.ncomponents),
-                    order="F",
+                    order=self.order,
                 )
                 loc = endloc
             else:
@@ -2533,7 +2540,7 @@ class ktensor:
         -------
         str:
         """
-        s = f"ktensor of shape {self.shape}\n"
+        s = f"ktensor of shape {self.shape} with order {self.order}\n"
         s += f"weights={self.weights}"
         if len(self.shape) == 0:
             s += "\nfactor_matrices=[]"
diff --git a/pyttb/pyttb_utils.py b/pyttb/pyttb_utils.py
index d96d9d2..870e589 100644
--- a/pyttb/pyttb_utils.py
+++ b/pyttb/pyttb_utils.py
@@ -510,7 +510,11 @@ def tt_ismember_rows(
     return matched, results.astype(int)
 
 
-def tt_ind2sub(shape: Tuple[int, ...], idx: np.ndarray) -> np.ndarray:
+def tt_ind2sub(
+    shape: Tuple[int, ...],
+    idx: np.ndarray,
+    order: Union[Literal["F"], Literal["C"]] = "F",
+) -> np.ndarray:
     """
     Multiple subscripts from linear indices.
 
@@ -526,7 +530,7 @@ def tt_ind2sub(shape: Tuple[int, ...], idx: np.ndarray) -> np.ndarray:
     if idx.size == 0:
         return np.empty(shape=(0, len(shape)), dtype=int)
     idx[idx < 0] += prod(shape)  # Handle negative indexing as simply as possible
-    return np.array(np.unravel_index(idx, shape, order="F")).transpose()
+    return np.array(np.unravel_index(idx, shape, order=order)).transpose()
 
 
 def tt_subsubsref(obj, s):
@@ -572,7 +576,11 @@ def tt_intvec2str(v: np.ndarray) -> str:
     return np.array2string(v)
 
 
-def tt_sub2ind(shape: Tuple[int, ...], subs: np.ndarray) -> np.ndarray:
+def tt_sub2ind(
+    shape: Tuple[int, ...],
+    subs: np.ndarray,
+    order: Union[Literal["F"], Literal["C"]] = "F",
+) -> np.ndarray:
     """
     Converts multidimensional subscripts to linear indices.
 
@@ -582,6 +590,8 @@ def tt_sub2ind(shape: Tuple[int, ...], subs: np.ndarray) -> np.ndarray:
         Shape of tensor
     subs:
         Subscripts for tensor
+    order:
+        Memory layout
 
     Returns
     -------
@@ -594,7 +604,7 @@ def tt_sub2ind(shape: Tuple[int, ...], subs: np.ndarray) -> np.ndarray:
     """
     if subs.size == 0:
         return np.array([])
-    idx = np.ravel_multi_index(tuple(subs.transpose()), shape, order="F")
+    idx = np.ravel_multi_index(tuple(subs.transpose()), shape, order=order)
     return idx
 
 
@@ -806,7 +816,7 @@ def get_index_variant(indices: IndexType) -> IndexVariant:
 
 
 def get_mttkrp_factors(
-    U: Union[ttb.ktensor, Sequence[np.ndarray]], n: int, ndims: int
+    U: Union[ttb.ktensor, Sequence[np.ndarray]], n: Union[int, np.integer], ndims: int
 ) -> Sequence[np.ndarray]:
     """Apply standard checks and type conversions for mttkrp factors"""
     if isinstance(U, ttb.ktensor):
diff --git a/pyttb/sptenmat.py b/pyttb/sptenmat.py
index 196f06b..c97ad50 100644
--- a/pyttb/sptenmat.py
+++ b/pyttb/sptenmat.py
@@ -343,7 +343,7 @@ class sptenmat:
         >>> S1[0, 0, 0] = 1
         >>> ST1 = S1.to_sptenmat(np.array([0]))
         >>> ST1.full()  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
diff --git a/pyttb/sptensor.py b/pyttb/sptensor.py
index a75ab7d..c111475 100644
--- a/pyttb/sptensor.py
+++ b/pyttb/sptensor.py
@@ -415,7 +415,7 @@ class sptensor:
 
         # Generate each column of the subscripts in turn
         for n in range(0, self.ndims):
-            i = o.copy()
+            i: list[np.ndarray] = o.copy()
             i[n] = np.expand_dims(np.arange(0, self.shape[n]), axis=1)
             s[:, n] = np.squeeze(ttb.khatrirao(*i))
 
@@ -516,7 +516,7 @@ class sptensor:
         >>> T = ttb.tensor(np.ones((2, 2, 2)))
         >>> S = T.to_sptensor()
         >>> S.contract(0, 1)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [2. 2.]
 
@@ -693,7 +693,7 @@ class sptensor:
         >>> S = ttb.sptensor()
         >>> S[1, 1] = 1
         >>> S.to_tensor()
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0. 0.]
          [0. 1.]]
@@ -807,7 +807,9 @@ class sptensor:
         csize = np.array(self.shape)[cdims]
 
         if rsize.size == 0:
-            ridx = np.zeros((self.nnz, 1))
+            ridx: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.zeros(
+                (self.nnz, 1), dtype=int
+            )
         elif self.subs.size == 0:
             ridx = np.array([], dtype=int)
         else:
@@ -815,7 +817,9 @@ class sptensor:
         ridx = ridx.reshape((ridx.size, 1)).astype(int)
 
         if csize.size == 0:
-            cidx = np.zeros((self.nnz, 1))
+            cidx: np.ndarray[tuple[int, ...], np.dtype[Any]] = np.zeros(
+                (self.nnz, 1), dtype=int
+            )
         elif self.subs.size == 0:
             cidx = np.array([], dtype=int)
         else:
@@ -1079,7 +1083,7 @@ class sptensor:
 
         >>> T = S.to_tensor()
         >>> S.logical_or(T)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1. 0.]
          [0. 1.]]
@@ -1087,7 +1091,7 @@ class sptensor:
         Compute logical OR with a scalar value:
 
         >>> S.logical_or(1)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1. 1.]
          [1. 1.]]
@@ -1149,7 +1153,7 @@ class sptensor:
         >>> T = S.to_tensor()
         >>> T[1, 0] = 1.0
         >>> S.logical_xor(T)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0. 0.]
          [1. 0.]]
@@ -1157,7 +1161,7 @@ class sptensor:
         Compute logical XOR with a scalar value:
 
         >>> S.logical_xor(1)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0. 1.]
          [1. 0.]]
@@ -1243,7 +1247,9 @@ class sptensor:
         vals[matching_indices] = self.vals[matching_indices]
         return vals
 
-    def mttkrp(self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: int) -> np.ndarray:
+    def mttkrp(
+        self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: Union[int, np.integer]
+    ) -> np.ndarray:
         """
         Matricized tensor times Khatri-Rao product using the
         :class:`pyttb.sptensor`. This is an efficient form of the matrix
@@ -1312,7 +1318,7 @@ class sptensor:
                 else:
                     Z.append(np.array([]))
             # Perform ttv multiplication
-            ttv = self.ttv(Z, exclude_dims=n)
+            ttv = self.ttv(Z, exclude_dims=int(n))
             # TODO is is possible to hit the float condition here?
             if isinstance(ttv, float):  # pragma: no cover
                 V[:, r] = ttv
@@ -1917,7 +1923,7 @@ class sptensor:
         result is a :class:`pyttb.tensor`:
 
         >>> S.ttv(np.ones(2), 0)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [4. 2.]
 
@@ -2799,7 +2805,7 @@ class sptensor:
         Subtract a scalar value, returning a dense tensor:
 
         >>> S - 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[-1. -1.]
          [-1.  0.]]
@@ -2848,7 +2854,7 @@ class sptensor:
         Add a scalar value, returning a dense tensor:
 
         >>> S + 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1. 1.]
          [1. 2.]]
@@ -3363,7 +3369,7 @@ class sptensor:
         >>> S = ttb.sptensor(shape=(2, 2))
         >>> S[:, :] = 2.0
         >>> 1 / S
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0.5 0.5]
          [0.5 0.5]]
@@ -3462,24 +3468,24 @@ class sptensor:
 
         >>> A = 2 * np.ones((2, 1))
         >>> S.ttm([A, A], dims=[0, 1], transpose=True)
-        tensor of shape (1, 1, 2, 2)
-        data[0, 0, :, :] =
-        [[8. 0.]
-         [8. 0.]]
+        tensor of shape (1, 1, 2, 2) with order F
+        data[:, :, 0, 0] =
+        [[8.]]
+        data[:, :, 1, 0] =
+        [[8.]]
+        data[:, :, 0, 1] =
+        [[0.]]
+        data[:, :, 1, 1] =
+        [[0.]]
 
         Compute sparse tensor matrix product specifying which two tensor
         dimensions to exclude in the multiplication:
 
         >>> S.ttm([A, A], exclude_dims=[0, 1], transpose=True)
-        tensor of shape (2, 2, 1, 1)
-        data[0, 0, :, :] =
-        [[8.]]
-        data[1, 0, :, :] =
-        [[8.]]
-        data[0, 1, :, :] =
-        [[0.]]
-        data[1, 1, :, :] =
-        [[0.]]
+        tensor of shape (2, 2, 1, 1) with order F
+        data[:, :, 0, 0] =
+        [[8. 0.]
+         [8. 0.]]
         """
         # Handle list of matrices
         if isinstance(matrices, Sequence):
diff --git a/pyttb/sumtensor.py b/pyttb/sumtensor.py
index 93e6834..771b028 100644
--- a/pyttb/sumtensor.py
+++ b/pyttb/sumtensor.py
@@ -108,7 +108,7 @@ class sumtensor:
         >>> ttb.sumtensor([T1, T2])  # doctest: +NORMALIZE_WHITESPACE
         sumtensor of shape (2, 2) with 2 parts:
         Part 0:
-            tensor of shape (2, 2)
+            tensor of shape (2, 2) with order F
             data[:, :] =
             [[1. 1.]
              [1. 1.]]
@@ -260,7 +260,7 @@ class sumtensor:
         >>> T = ttb.tenones((2, 2))
         >>> S = ttb.sumtensor([T, T])
         >>> print(S.full())  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2. 2.]
          [2. 2.]]
@@ -318,7 +318,9 @@ class sumtensor:
             result += part.innerprod(other)
         return result
 
-    def mttkrp(self, U: Union[ttb.ktensor, List[np.ndarray]], n: int) -> np.ndarray:
+    def mttkrp(
+        self, U: Union[ttb.ktensor, List[np.ndarray]], n: Union[int, np.integer]
+    ) -> np.ndarray:
         """
         Matricized tensor times Khatri-Rao product. The matrices used in the
         Khatri-Rao product are passed as a :class:`pyttb.ktensor` (where the
@@ -396,17 +398,17 @@ class sumtensor:
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> S = ttb.sumtensor([T, T])
         >>> T.ttv(np.ones(2), 0)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [4. 6.]
         >>> S.ttv(np.ones(2), 0)  # doctest: +NORMALIZE_WHITESPACE
         sumtensor of shape (2,) with 2 parts:
         Part 0:
-             tensor of shape (2,)
+             tensor of shape (2,) with order F
              data[:] =
              [4. 6.]
         Part 1:
-             tensor of shape (2,)
+             tensor of shape (2,) with order F
              data[:] =
              [4. 6.]
         >>> T.ttv([np.ones(2), np.ones(2)])
diff --git a/pyttb/tenmat.py b/pyttb/tenmat.py
index d9f294b..531be95 100644
--- a/pyttb/tenmat.py
+++ b/pyttb/tenmat.py
@@ -7,7 +7,7 @@
 from __future__ import annotations
 
 from math import prod
-from typing import Optional, Tuple, Union
+from typing import Literal, Optional, Tuple, Union
 
 import numpy as np
 
@@ -53,7 +53,7 @@ class tenmat:
         Create an empty :class:`pyttb.tenmat`.
 
         >>> ttb.tenmat()  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape ()
+        matrix corresponding to a tensor of shape () with order F
         rindices = [  ] (modes of tensor corresponding to rows)
         cindices = [  ] (modes of tensor corresponding to columns)
         data = []
@@ -88,7 +88,7 @@ class tenmat:
                 [6., 7.]]])
         """
 
-        # Case 0a: Empty Contructor
+        # Case 0a: Empty Constructor
         # data is empty, return empty tenmat unless rdims, cdims, or tshape are
         # not empty
         if data is None or (isinstance(data, np.ndarray) and data.size == 0):
@@ -102,7 +102,7 @@ class tenmat:
             self.tshape: Union[Tuple[()], Tuple[int, ...]] = ()
             self.rindices = np.array([])
             self.cindices = np.array([])
-            self.data = np.array([], ndmin=2, order="F")
+            self.data = np.array([], ndmin=2, order=self.order)
             return
 
         # Verify that data is a numeric numpy.ndarray
@@ -116,7 +116,7 @@ class tenmat:
                 assert False, "tshape must be specified when data is 1d array."
             else:
                 # make data a 2d array with shape (1, data.shape[0]), i.e., a row vector
-                data = np.reshape(data.copy(), (1, data.shape[0]), order="F")
+                data = np.reshape(data.copy(), (1, data.shape[0]), order=self.order)
 
         if len(data.shape) != 2:
             raise ValueError(
@@ -169,6 +169,11 @@ class tenmat:
             self.data = data
         return
 
+    @property
+    def order(self) -> Literal["F"]:
+        """Return the data layout of the underlying storage."""
+        return "F"
+
     def copy(self) -> tenmat:
         """
         Return a deep copy of the :class:`pyttb.tenmat`.
@@ -234,13 +239,13 @@ class tenmat:
         Extract original tensor shaped data.
 
         >>> tm.to_tensor()  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2, 2)
-        data[0, :, :] =
-        [[0. 1.]
-         [2. 3.]]
-        data[1, :, :] =
-        [[4. 5.]
-         [6. 7.]]
+        tensor of shape (2, 2, 2) with order F
+        data[:, :, 0] =
+        [[0. 2.]
+         [4. 6.]]
+        data[:, :, 1] =
+        [[1. 3.]
+         [5. 7.]]
         """
         # RESHAPE TENSOR-AS-MATRIX
         # Here we just reverse what was done in the tenmat constructor.
@@ -251,7 +256,7 @@ class tenmat:
         data = self.data
         if copy:
             data = self.data.copy()
-        data = np.reshape(data, np.array(shape)[order], order="F")
+        data = np.reshape(data, np.array(shape)[order], order=self.order)
         if order.size > 1:
             data = np.transpose(data, np.argsort(order))
         return ttb.tensor(data, shape, copy=False)
@@ -267,14 +272,14 @@ class tenmat:
         >>> T = ttb.tenones((2, 2, 2))
         >>> TM = T.to_tenmat(rdims=np.array([0]))
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
         [[1. 1. 1. 1.]
          [1. 1. 1. 1.]]
         >>> TM.ctranspose()  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 1, 2 ] (modes of tensor corresponding to rows)
         cindices = [ 0 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -300,7 +305,7 @@ class tenmat:
         >>> T = ttb.tenones((2, 2, 2))
         >>> TM = T.to_tenmat(rdims=np.array([0]))
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -341,7 +346,7 @@ class tenmat:
         >>> T = ttb.tenones((2, 2, 2))
         >>> TM = T.to_tenmat(rdims=np.array([0]))
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -406,7 +411,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2, 2)).to_tenmat(np.array([0]))
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -414,7 +419,7 @@ class tenmat:
          [1. 1. 1. 1.]]
         >>> TM[0, 0] = 2.0
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -451,7 +456,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> TM * TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -508,7 +513,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> TM * TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -533,14 +538,14 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> TM + TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
         [[2. 2.]
          [2. 2.]]
         >>> TM + 1.0  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2)  with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -579,7 +584,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> 1.0 + TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -604,14 +609,14 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> TM - TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
         [[0. 0.]
          [0. 0.]]
         >>> TM - 1.0  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -650,7 +655,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> 1.0 - TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -685,7 +690,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> +TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -710,7 +715,7 @@ class tenmat:
         --------
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> -TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -737,7 +742,7 @@ class tenmat:
         Print an empty :class:`pyttb.tenmat`.
 
         >>> ttb.tenmat()  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape ()
+        matrix corresponding to a tensor of shape () with order F
         rindices = [  ] (modes of tensor corresponding to rows)
         cindices = [  ] (modes of tensor corresponding to columns)
         data = []
@@ -746,7 +751,7 @@ class tenmat:
 
         >>> TM = ttb.tenones((2, 2)).to_tenmat(np.array([0]))
         >>> TM  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2)
+        matrix corresponding to a tensor of shape (2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -762,6 +767,7 @@ class tenmat:
         s = ""
         s += "matrix corresponding to a tensor of shape "
         s += str(np_to_python(self.tshape))
+        s += f" with order {self.order}"
         s += "\n"
 
         s += "rindices = "
diff --git a/pyttb/tensor.py b/pyttb/tensor.py
index 80d1045..b527f94 100644
--- a/pyttb/tensor.py
+++ b/pyttb/tensor.py
@@ -19,6 +19,7 @@ from typing import (
     Sequence,
     Tuple,
     Union,
+    cast,
     overload,
 )
 
@@ -104,14 +105,14 @@ class tensor:
 
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> print(T)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 2]
          [3 4]]
         """
         if data is None:
             # EMPTY / DEFAULT CONSTRUCTOR
-            self.data: np.ndarray = np.array([])
+            self.data: np.ndarray = np.array([], order=self.order)
             self.shape: Tuple = ()
             return
 
@@ -139,16 +140,33 @@ class tensor:
         # Make sure the data is indeed the right shape
         if data.size > 0 and len(shape) > 0:
             # reshaping using Fortran ordering to match Matlab conventions
-            data = np.reshape(data, np.array(shape), order="F")
+            data = np.reshape(data, np.array(shape), order=self.order)
 
         # Create the tensor
         if copy:
-            self.data = data.copy()
+            self.data = data.copy(self.order)
         else:
+            if not self.matches_order(data):
+                logging.warning(
+                    f"Selected no copy, but input data isn't {self.order} ordered "
+                    "so must copy."
+                )
             self.data = data
         self.shape = shape
         return
 
+    @property
+    def order(self) -> Literal["F"]:
+        """Return the data layout of the underlying storage."""
+        return "F"
+
+    def matches_order(self, array: np.ndarray) -> bool:
+        if array.flags["C_CONTIGUOUS"] and self.order == "C":
+            return True
+        if array.flags["F_CONTIGUOUS"] and self.order == "F":
+            return True
+        return False
+
     @classmethod
     def from_function(
         cls,
@@ -178,15 +196,19 @@ class tensor:
 
         >>> T = ttb.tensor.from_function(np.ones, (2, 3, 4))
         >>> print(T)
-        tensor of shape (2, 3, 4)
-        data[0, :, :] =
-        [[1. 1. 1. 1.]
-         [1. 1. 1. 1.]
-         [1. 1. 1. 1.]]
-        data[1, :, :] =
-        [[1. 1. 1. 1.]
-         [1. 1. 1. 1.]
-         [1. 1. 1. 1.]]
+        tensor of shape (2, 3, 4) with order F
+        data[:, :, 0] =
+        [[1. 1. 1.]
+         [1. 1. 1.]]
+        data[:, :, 1] =
+        [[1. 1. 1.]
+         [1. 1. 1.]]
+        data[:, :, 2] =
+        [[1. 1. 1.]
+         [1. 1. 1.]]
+        data[:, :, 3] =
+        [[1. 1. 1.]
+         [1. 1. 1.]]
         """
         # Check size
         shape = parse_shape(shape)
@@ -245,7 +267,7 @@ class tensor:
         >>> T.collapse()
         4.0
         >>> T.collapse(np.array([0]))
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [2. 2.]
         >>> T.collapse(np.arange(T.ndims), sum)
@@ -254,7 +276,7 @@ class tensor:
         1.0
         """
         if self.data.size == 0:
-            return np.array([])
+            return np.array([], order=self.order)
 
         if dims is None:
             dims = np.arange(0, self.ndims)
@@ -267,7 +289,7 @@ class tensor:
 
         # Check for the case where we accumulate over *all* dimensions
         if remdims.size == 0:
-            result = fun(self.data.flatten("F"))
+            result = fun(self.data.flatten(self.order))
             if isinstance(result, np.generic):
                 result = result.item()
             return result
@@ -279,7 +301,7 @@ class tensor:
         A = self.to_tenmat(remdims, dims).double()
 
         ## Apply the collapse function
-        B = np.zeros((A.shape[0], 1))
+        B = np.zeros((A.shape[0], 1), order=self.order)
         for i in range(0, A.shape[0]):
             B[i] = fun(A[i, :])
 
@@ -308,23 +330,23 @@ class tensor:
         2.0
         >>> T = ttb.tensor(np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]]))
         >>> print(T)
-        tensor of shape (2, 2, 2)
-        data[0, :, :] =
-        [[1 2]
-         [3 4]]
-        data[1, :, :] =
-        [[5 6]
-         [7 8]]
+        tensor of shape (2, 2, 2) with order F
+        data[:, :, 0] =
+        [[1 3]
+         [5 7]]
+        data[:, :, 1] =
+        [[2 4]
+         [6 8]]
         >>> T.contract(0, 1)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [ 8. 10.]
         >>> T.contract(0, 2)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [ 7. 11.]
         >>> T.contract(1, 2)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [ 5. 13.]
         """
@@ -354,16 +376,16 @@ class tensor:
         x = self.permute(np.concatenate((remdims, np.array([i1, i2]))))
 
         # Reshape data to be 3D
-        data = np.reshape(x.data, (m, n, n), order="F")
+        data = np.reshape(x.data, (m, n, n), order=self.order)
 
         # Add diagonal entries for each slice
-        newdata = np.zeros((m, 1))
+        newdata = np.zeros((m, 1), order=self.order)
         for idx in range(0, n):
             newdata += data[:, idx, idx][:, None]
 
         # Reshape result
         if prod(newsize) > 1:
-            newdata = np.reshape(newdata, newsize, order="F")
+            newdata = np.reshape(newdata, newsize, order=self.order)
 
         return ttb.tensor(newdata, newsize, copy=False)
 
@@ -415,7 +437,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1,2],[3,4]]))
         >>> print(T)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 2]
          [3 4]]
@@ -428,7 +450,7 @@ class tensor:
         array([[ True],
                [ True]])
         """
-        idx = np.nonzero(np.ravel(self.data, order="F"))[0]
+        idx = np.nonzero(np.ravel(self.data, order=self.order))[0]
         subs = tt_ind2sub(self.shape, idx)
         vals = self.data[tuple(subs.T)][:, None]
         return subs, vals
@@ -445,7 +467,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[0, 2], [3, 0]]))
         >>> print(T)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0 2]
          [3 0]]
@@ -518,13 +540,13 @@ class tensor:
         >>> data = np.reshape(np.arange(prod(tshape)), tshape)
         >>> T = ttb.tensor(data)
         >>> T  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2, 2)
-        data[0, :, :] =
-        [[0 1]
-         [2 3]]
-        data[1, :, :] =
-        [[4 5]
-         [6 7]]
+        tensor of shape (2, 2, 2) with order F
+        data[:, :, 0] =
+        [[0 2]
+         [4 6]]
+        data[:, :, 1] =
+        [[1 3]
+         [5 7]]
 
         Convert to a :class:`pyttb.tenmat` unwrapping around the first dimension.
             Either allow for implicit column or explicit column dimension
@@ -540,7 +562,7 @@ class tensor:
 
         >>> TM3 = T.to_tenmat(rdims=np.array([0]), cdims_cyclic="fc")
         >>> TM3  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2)  with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 1, 2 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -551,7 +573,7 @@ class tensor:
 
         >>> TM4 = T.to_tenmat(rdims=np.array([0]), cdims_cyclic="bc")
         >>> TM4  # doctest: +NORMALIZE_WHITESPACE
-        matrix corresponding to a tensor of shape (2, 2, 2)
+        matrix corresponding to a tensor of shape (2, 2, 2) with order F
         rindices = [ 0 ] (modes of tensor corresponding to rows)
         cindices = [ 2, 1 ] (modes of tensor corresponding to columns)
         data[:, :] =
@@ -588,7 +610,7 @@ class tensor:
         data = np.reshape(
             self.permute(dims).data,
             (rprod, cprod),
-            order="F",
+            order=self.order,
         )
         return ttb.tenmat(data, rdims, cdims, tshape=tshape, copy=copy)
 
@@ -616,8 +638,8 @@ class tensor:
         if isinstance(other, ttb.tensor):
             if self.shape != other.shape:
                 assert False, "Inner product must be between tensors of the same size"
-            x = np.reshape(self.data, (self.data.size,), order="F")
-            y = np.reshape(other.data, (other.data.size,), order="F")
+            x = np.reshape(self.data, (self.data.size,), order=self.order)
+            y = np.reshape(other.data, (other.data.size,), order=self.order)
             return x.dot(y).item()
         if isinstance(other, (ttb.ktensor, ttb.sptensor, ttb.ttensor)):
             # Reverse arguments and call specializer code
@@ -712,9 +734,8 @@ class tensor:
         sz = np.array(self.shape)
 
         if grps is None:
-            grps = np.arange(0, n)
-
-        if len(grps.shape) == 1:
+            grps = np.arange(0, n)[None, :]
+        elif len(grps.shape) == 1:
             grps = np.array([grps])
 
         # Substantially different routines are called depending on whether the user
@@ -881,7 +902,9 @@ class tensor:
         # Extract those non-zero values
         return self.data[tuple(wsubs.transpose())]
 
-    def mttkrp(self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: int) -> np.ndarray:
+    def mttkrp(
+        self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: Union[int, np.integer]
+    ) -> np.ndarray:
         """
         Matricized tensor times Khatri-Rao product. The matrices used in the
         Khatri-Rao product are passed as a :class:`pyttb.ktensor` (where the
@@ -931,20 +954,20 @@ class tensor:
 
         if n == 0:
             Ur = ttb.khatrirao(*U[1 : self.ndims], reverse=True)
-            Y = np.reshape(self.data, (szn, szr), order="F")
+            Y = np.reshape(self.data, (szn, szr), order=self.order)
             return Y @ Ur
         if n == self.ndims - 1:
             Ul = ttb.khatrirao(*U[0 : self.ndims - 1], reverse=True)
-            Y = np.reshape(self.data, (szl, szn), order="F")
+            Y = np.reshape(self.data, (szl, szn), order=self.order)
             return Y.T @ Ul
         else:
             Ul = ttb.khatrirao(*U[n + 1 :], reverse=True)
             Ur = np.reshape(
-                ttb.khatrirao(*U[0:n], reverse=True), (szl, 1, R), order="F"
+                ttb.khatrirao(*U[0:n], reverse=True), (szl, 1, R), order=self.order
             )
-            Y = np.reshape(self.data, (-1, szr), order="F")
+            Y = np.reshape(self.data, (-1, szr), order=self.order)
             Y = Y @ Ul
-            Y = np.reshape(Y, (szl, szn, R), order="F")
+            Y = np.reshape(Y, (szl, szn, R), order=self.order)
             V = np.zeros((szn, R))
             for r in range(R):
                 V[:, [r]] = Y[:, :, r].T @ Ur[:, :, r]
@@ -980,14 +1003,14 @@ class tensor:
         split_idx = min_split(self.shape)
         V = [np.empty_like(self.data, shape=())] * self.ndims
         K = ttb.khatrirao(*U[split_idx + 1 :], reverse=True)
-        W = np.reshape(self.data, (-1, K.shape[0]), order="F").dot(K)
+        W = np.reshape(self.data, (-1, K.shape[0]), order=self.order).dot(K)
         for k in range(split_idx):
             # Loop entry invariant: W has modes (mk x ... x ms, C)
             V[k] = mttv_mid(W, U[k + 1 : split_idx + 1])
             W = mttv_left(W, U[k])
         V[split_idx] = W
         K = ttb.khatrirao(*U[0 : split_idx + 1], reverse=True)
-        W = np.reshape(self.data, (K.shape[0], -1), order="F").transpose().dot(K)
+        W = np.reshape(self.data, (K.shape[0], -1), order=self.order).transpose().dot(K)
         for k in range(split_idx + 1, self.ndims - 1):
             # Loop invariant: W has modes (mk x .. x md, C)
             V[k] = mttv_mid(W, U[k + 1 :])
@@ -1116,12 +1139,12 @@ class tensor:
         --------
         >>> T1 = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 2]
          [3 4]]
         >>> T1.permute(np.array((1, 0)))
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 3]
          [2 4]]
@@ -1163,7 +1186,7 @@ class tensor:
         if prod(self.shape) != prod(shape):
             assert False, "Reshaping a tensor cannot change number of elements"
 
-        return ttb.tensor(np.reshape(self.data, shape, order="F"), shape)
+        return ttb.tensor(np.reshape(self.data, shape, order=self.order), shape)
 
     def scale(
         self,
@@ -1250,7 +1273,9 @@ class tensor:
         else:
             idx = np.where(shapeArray > 1)
             if idx[0].size == 0:
-                return self.data.item()
+                # Why is item annotated as str?
+                single_item: float = cast(float, self.data.item())
+                return single_item
             return ttb.tensor(np.squeeze(self.data))
 
     def symmetrize(  # noqa: PLR0912,PLR0915
@@ -1277,11 +1302,11 @@ class tensor:
         --------
         >>> T = ttb.tenones((2, 2, 2))
         >>> T.symmetrize(np.array([0, 2]))
-        tensor of shape (2, 2, 2)
-        data[0, :, :] =
+        tensor of shape (2, 2, 2) with order F
+        data[:, :, 0] =
         [[1. 1.]
          [1. 1.]]
-        data[1, :, :] =
+        data[:, :, 1] =
         [[1. 1.]
          [1. 1.]]
         """
@@ -1443,20 +1468,20 @@ class tensor:
         >>> T = ttb.tenones((2, 2, 2, 2))
         >>> A = 2 * np.ones((2, 1))
         >>> T.ttm([A, A], dims=[0, 1], transpose=True)
-        tensor of shape (1, 1, 2, 2)
-        data[0, 0, :, :] =
-        [[16. 16.]
-         [16. 16.]]
-        >>> T.ttm([A, A], exclude_dims=[0, 1], transpose=True)
-        tensor of shape (2, 2, 1, 1)
-        data[0, 0, :, :] =
+        tensor of shape (1, 1, 2, 2) with order F
+        data[:, :, 0, 0] =
         [[16.]]
-        data[1, 0, :, :] =
+        data[:, :, 1, 0] =
         [[16.]]
-        data[0, 1, :, :] =
+        data[:, :, 0, 1] =
         [[16.]]
-        data[1, 1, :, :] =
+        data[:, :, 1, 1] =
         [[16.]]
+        >>> T.ttm([A, A], exclude_dims=[0, 1], transpose=True)
+        tensor of shape (2, 2, 1, 1) with order F
+        data[:, :, 0, 0] =
+        [[16. 16.]
+         [16. 16.]]
         """
         if isinstance(matrix, Sequence):
             # Check that the dimensions are valid
@@ -1485,7 +1510,7 @@ class tensor:
         second_dim = 1
         if len(ids) > 0:
             second_dim = np.prod(shape[ids])
-        newdata = np.reshape(newdata, (shape[n], second_dim), order="F")
+        newdata = np.reshape(newdata, (shape[n], second_dim), order=self.order)
         if transpose:
             newdata = matrix.T @ newdata
             p = matrix.shape[1]
@@ -1496,7 +1521,7 @@ class tensor:
         newshape = np.array(
             [p, *list(shape[range(0, n)]), *list(shape[range(n + 1, self.ndims)])]
         )
-        Y_data = np.reshape(newdata, newshape, order="F")
+        Y_data: np.ndarray = np.reshape(newdata, newshape, order=self.order)
         Y_data = np.transpose(Y_data, np.argsort(order))
         return ttb.tensor(Y_data, copy=False)
 
@@ -1533,26 +1558,26 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T.ttt(T)
-        tensor of shape (2, 2, 2, 2)
-        data[0, 0, :, :] =
+        tensor of shape (2, 2, 2, 2) with order F
+        data[:, :, 0, 0] =
         [[1 2]
          [3 4]]
-        data[1, 0, :, :] =
+        data[:, :, 1, 0] =
         [[ 3  6]
          [ 9 12]]
-        data[0, 1, :, :] =
+        data[:, :, 0, 1] =
         [[2 4]
          [6 8]]
-        data[1, 1, :, :] =
+        data[:, :, 1, 1] =
         [[ 4  8]
          [12 16]]
         >>> T.ttt(T, 0)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[10 14]
          [14 20]]
         >>> T.ttt(T, selfdims=0, otherdims=1)
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ 7 15]
          [10 22]]
@@ -1630,11 +1655,11 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T.ttv(np.ones(2), 0)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [4. 6.]
         >>> T.ttv(np.ones(2), 1)
-        tensor of shape (2,)
+        tensor of shape (2,) with order F
         data[:] =
         [3. 7.]
         >>> T.ttv([np.ones(2), np.ones(2)])
@@ -1667,7 +1692,9 @@ class tensor:
         sz = np.array(self.shape)[np.concatenate((remdims, dims))]
 
         for i in range(dims.size - 1, -1, -1):
-            c = np.reshape(c, tuple([np.prod(sz[0 : n - 1]), sz[n - 1]]), order="F")
+            c = np.reshape(
+                c, tuple([np.prod(sz[0 : n - 1]), sz[n - 1]]), order=self.order
+            )
             c = c.dot(vector[vidx[i]])
             n -= 1
         # If needed, convert the final result back to tensor
@@ -1734,21 +1761,23 @@ class tensor:
 
             y = self.data.copy()
             for i in range(drem, 0, -1):
-                yy = np.reshape(y, (sz ** (dnew + i - 1), sz), order="F")
+                yy = np.reshape(y, (sz ** (dnew + i - 1), sz), order=self.order)
                 y = yy.dot(vector)
 
             # Convert to matrix if 2-way or convert back to tensor if result is >= 3-way
             if dnew == 2:
-                return np.reshape(y, [sz, sz], order="F")
+                return np.reshape(y, [sz, sz], order=self.order)
             if dnew > 2:
                 return ttb.tensor(
-                    np.reshape(y, newshape=sz * np.ones(dnew, dtype=int), order="F"),
+                    np.reshape(
+                        y, newshape=sz * np.ones(dnew, dtype=int), order=self.order
+                    ),
                     copy=False,
                 )
 
             # extract scalar if needed
             if len(y) == 1:
-                y = y.item()
+                return cast(float, y.item())
 
             return y
         assert False, "Invalid value for version; should be None, 1, or 2"
@@ -1923,18 +1952,15 @@ class tensor:
         1.0
         >>> # produces a tensor of order 1 and size 1
         >>> T[1, 1, 1, :]  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (1,)
+        tensor of shape (1,) with order F
         data[:] =
         [1.]
         >>> # produces a tensor of size 2 x 2 x 1
         >>> T[0:2, [2, 3], 1, :]  # doctest: +NORMALIZE_WHITESPACE
-        tensor of shape (2, 2, 1)
-        data[0, :, :] =
-        [[1.]
-         [1.]]
-        data[1, :, :] =
-        [[1.]
-         [1.]]
+        tensor of shape (2, 2, 1) with order F
+        data[:, :, 0] =
+        [[1. 1.]
+         [1. 1.]]
         >>> # returns a vector of length 2
         >>> # Equivalent to selecting [0,0,0,0] and [1,1,1,0] separately
         >>> T[np.array([[0, 0, 0, 0], [1, 1, 1, 0]])]
@@ -2031,12 +2057,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T == T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True  True]
          [ True  True]]
         >>> T == 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True False]
          [False False]]
@@ -2063,12 +2089,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T != T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False False]
          [False False]]
         >>> T != 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False  True]
          [ True  True]]
@@ -2095,12 +2121,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T >= T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True  True]
          [ True  True]]
         >>> T >= 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True  True]
          [ True  True]]
@@ -2127,12 +2153,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T <= T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True  True]
          [ True  True]]
         >>> T <= 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ True False]
          [False False]]
@@ -2159,12 +2185,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T > T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False False]
          [False False]]
         >>> T > 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False  True]
          [ True  True]]
@@ -2191,12 +2217,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T < T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False False]
          [False False]]
         >>> T < 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[False False]
          [False False]]
@@ -2223,12 +2249,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T - T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0 0]
          [0 0]]
         >>> T - 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0 1]
          [2 3]]
@@ -2255,12 +2281,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T + T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2 4]
          [6 8]]
         >>> T + 1
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2 3]
          [4 5]]
@@ -2290,7 +2316,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> 1 + T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2 3]
          [4 5]]
@@ -2313,7 +2339,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T**2
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ 1  4]
          [ 9 16]]
@@ -2340,12 +2366,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T * T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[ 1  4]
          [ 9 16]]
         >>> T * 2
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2 4]
          [6 8]]
@@ -2375,7 +2401,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> 2 * T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2 4]
          [6 8]]
@@ -2398,12 +2424,12 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T / T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1. 1.]
          [1. 1.]]
         >>> T / 2
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[0.5 1. ]
          [1.5 2. ]]
@@ -2434,7 +2460,7 @@ class tensor:
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> np.set_printoptions(precision=8)
         >>> 2 / T  # doctest: +ELLIPSIS
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[2.         1.        ]
          [0.66666... 0.5       ]]
@@ -2460,7 +2486,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> +T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 2]
          [3 4]]
@@ -2479,7 +2505,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> -T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[-1 -2]
          [-3 -4]]
@@ -2499,7 +2525,7 @@ class tensor:
         --------
         >>> T = ttb.tensor(np.array([[1, 2], [3, 4]]))
         >>> T
-        tensor of shape (2, 2)
+        tensor of shape (2, 2) with order F
         data[:, :] =
         [[1 2]
          [3 4]]
@@ -2512,7 +2538,7 @@ class tensor:
             return s
 
         s = ""
-        s += f"tensor of shape {np_to_python(self.shape)}"
+        s += f"tensor of shape {np_to_python(self.shape)} with order {self.order}"
 
         if self.ndims == 1:
             s += "\ndata"
@@ -2521,22 +2547,22 @@ class tensor:
                 s += " =\n"
                 s += str(self.data)
                 return s
-        for i in np.arange(np.prod(self.shape[:-2])):
+        for i in np.arange(np.prod(self.shape[2:])):
             s += "\ndata"
             if self.ndims == 2:
                 s += "[:, :]"
                 s += " =\n"
                 s += str(self.data)
             elif self.ndims > 2:
-                idx = tt_ind2sub(self.shape[:-2], np.array([i]))
-                s += str(idx[0].tolist())[0:-1]
-                s += ", :, :]"
+                idx = tt_ind2sub(self.shape[2:], np.array([i]))
+                s += "[:, :, "
+                s += str(idx[0].tolist())[1:]
                 s += " =\n"
                 s += str(
                     self.data[
                         tuple(
                             np.concatenate(
-                                (idx[0], np.array([slice(None), slice(None)]))
+                                (np.array([slice(None), slice(None)]), idx[0])
                             )
                         )
                     ]
@@ -2547,7 +2573,7 @@ class tensor:
     __str__ = __repr__
 
 
-def tenones(shape: Shape) -> tensor:
+def tenones(shape: Shape, order: Union[Literal["F"], Literal["C"]] = "F") -> tensor:
     """
     Creates a tensor of all ones.
 
@@ -2555,6 +2581,8 @@ def tenones(shape: Shape) -> tensor:
     ----------
     shape:
         Shape of resulting tensor.
+    order:
+        Memory layout for resulting tensor.
 
     Returns
     -------
@@ -2564,21 +2592,25 @@ def tenones(shape: Shape) -> tensor:
     --------
     >>> T = ttb.tenones((3,))
     >>> T
-    tensor of shape (3,)
+    tensor of shape (3,) with order F
     data[:] =
     [1. 1. 1.]
     >>> T = ttb.tenones((3, 3))
     >>> T
-    tensor of shape (3, 3)
+    tensor of shape (3, 3) with order F
     data[:, :] =
     [[1. 1. 1.]
      [1. 1. 1.]
      [1. 1. 1.]]
     """
-    return tensor.from_function(np.ones, shape)
 
+    def ones(shape: Tuple[int, ...]) -> np.ndarray:
+        return np.ones(shape, order=order)
 
-def tenzeros(shape: Shape) -> tensor:
+    return tensor.from_function(ones, shape)
+
+
+def tenzeros(shape: Shape, order: Union[Literal["F"], Literal["C"]] = "F") -> tensor:
     """
     Creates a tensor of all zeros.
 
@@ -2586,6 +2618,8 @@ def tenzeros(shape: Shape) -> tensor:
     ----------
     shape:
         Shape of resulting tensor.
+    order:
+        Memory layout for resulting tensor.
 
     Returns
     -------
@@ -2595,21 +2629,25 @@ def tenzeros(shape: Shape) -> tensor:
     --------
     >>> T = ttb.tenzeros((3,))
     >>> T
-    tensor of shape (3,)
+    tensor of shape (3,) with order F
     data[:] =
     [0. 0. 0.]
     >>> T = ttb.tenzeros((3, 3))
     >>> T
-    tensor of shape (3, 3)
+    tensor of shape (3, 3) with order F
     data[:, :] =
     [[0. 0. 0.]
      [0. 0. 0.]
      [0. 0. 0.]]
     """
-    return tensor.from_function(np.zeros, shape)
+
+    def zeros(shape: Tuple[int, ...]) -> np.ndarray:
+        return np.zeros(shape, order=order)
+
+    return tensor.from_function(zeros, shape)
 
 
-def tenrand(shape: Shape) -> tensor:
+def tenrand(shape: Shape, order: Union[Literal["F"], Literal["C"]] = "F") -> tensor:
     """
     Creates a tensor with entries drawn from a uniform
     distribution on the unit interval.
@@ -2618,6 +2656,8 @@ def tenrand(shape: Shape) -> tensor:
     ----------
     shape:
         Shape of resulting tensor.
+    order:
+        Memory layout for resulting tensor.
 
     Returns
     -------
@@ -2628,7 +2668,7 @@ def tenrand(shape: Shape) -> tensor:
     >>> np.random.seed(1)
     >>> T = ttb.tenrand((3,))
     >>> T
-    tensor of shape (3,)
+    tensor of shape (3,) with order F
     data[:] =
     [4.170...e-01 7.203...e-01 1.143...e-04]
     """
@@ -2636,12 +2676,19 @@ def tenrand(shape: Shape) -> tensor:
     # Typing doesn't play nice with partial
     # mypy issue: 1484
     def unit_uniform(pass_through_shape: Tuple[int, ...]) -> np.ndarray:
-        return np.random.uniform(low=0, high=1, size=pass_through_shape)
+        data = np.random.uniform(low=0, high=1, size=pass_through_shape)
+        if order == "F":
+            return np.asfortranarray(data)
+        return data
 
     return tensor.from_function(unit_uniform, shape)
 
 
-def tendiag(elements: OneDArray, shape: Optional[Shape] = None) -> tensor:
+def tendiag(
+    elements: OneDArray,
+    shape: Optional[Shape] = None,
+    order: Union[Literal["F"], Literal["C"]] = "F",
+) -> tensor:
     """
     Creates a tensor with elements along super diagonal. If provided shape is too
     small the tensor will be enlarged to accomodate.
@@ -2652,6 +2699,8 @@ def tendiag(elements: OneDArray, shape: Optional[Shape] = None) -> tensor:
         Elements to set along the diagonal.
     shape:
         Shape of resulting tensor.
+    order:
+        Memory layout for resulting tensor.
 
     Returns
     -------
@@ -2674,13 +2723,15 @@ def tendiag(elements: OneDArray, shape: Optional[Shape] = None) -> tensor:
     else:
         shape = parse_shape(shape)
         constructed_shape = tuple(max(N, dim) for dim in shape)
-    X = tenzeros(constructed_shape)
+    X = tenzeros(constructed_shape, order=order)
     subs = np.tile(np.arange(0, N)[:, None], (len(constructed_shape),))
     X[subs] = elements
     return X
 
 
-def teneye(order: int, size: int) -> tensor:
+def teneye(
+    ndims: int, size: int, order: Union[Literal["F"], Literal["C"]] = "F"
+) -> tensor:
     """Create identity tensor of specified shape.
 
     T is an "identity tensor if T.ttsv(x, skip_dim=0) = x for all x such that
@@ -2692,13 +2743,15 @@ def teneye(order: int, size: int) -> tensor:
 
     Parameters
     ----------
-    order: Number of dimensions of tensor.
+    ndims: Number of dimensions of tensor.
     size: Number of elements in any dimension of the tensor.
+    order:
+        Memory layout for resulting tensor.
 
     Examples
     --------
     >>> ttb.teneye(2, 3)
-    tensor of shape (3, 3)
+    tensor of shape (3, 3) with order F
     data[:, :] =
     [[1. 0. 0.]
      [0. 1. 0.]
@@ -2713,17 +2766,17 @@ def teneye(order: int, size: int) -> tensor:
     -------
     Identity tensor.
     """
-    if order % 2 != 0:
-        raise ValueError(f"Order must be even but received {order}")
-    idx_iterator = combinations_with_replacement(range(size), order)
-    A = tenzeros((size,) * order)
-    s = np.zeros((factorial(order), order // 2))
+    if ndims % 2 != 0:
+        raise ValueError(f"Order must be even but received {ndims}")
+    idx_iterator = combinations_with_replacement(range(size), ndims)
+    A = tenzeros((size,) * ndims, order=order)
+    s = np.zeros((factorial(ndims), ndims // 2), order=order)
     for _i, indices in enumerate(idx_iterator):
         p = np.array(list(permutations(indices)))
-        for j in range(order // 2):
+        for j in range(ndims // 2):
             s[:, j] = p[:, 2 * j - 1] == p[:, 2 * j]
-        v = np.sum(np.sum(s, axis=1) == order // 2)
-        A[tuple(zip(*p))] = v / factorial(order)
+        v = np.sum(np.sum(s, axis=1) == ndims // 2)
+        A[tuple(zip(*p))] = v / factorial(ndims)
     return A
 
 
diff --git a/pyttb/ttensor.py b/pyttb/ttensor.py
index c33c8ae..cab3665 100644
--- a/pyttb/ttensor.py
+++ b/pyttb/ttensor.py
@@ -401,7 +401,9 @@ class ttensor:
         assert not isinstance(newcore, float)
         return ttensor(newcore, [self.factor_matrices[dim] for dim in remdims])
 
-    def mttkrp(self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: int) -> np.ndarray:
+    def mttkrp(
+        self, U: Union[ttb.ktensor, Sequence[np.ndarray]], n: Union[int, np.integer]
+    ) -> np.ndarray:
         """
         Matricized tensor times Khatri-Rao product for ttensors.
 


Test patch:
diff --git a/tests/test_ktensor.py b/tests/test_ktensor.py
index 30c0af8..a31f30a 100644
--- a/tests/test_ktensor.py
+++ b/tests/test_ktensor.py
@@ -1228,5 +1228,5 @@ def test_ktensor__mul__(sample_ktensor_2way, sample_ktensor_3way):
 
 def test_ktensor__str__(sample_ktensor_2way):
     (data0, K0) = sample_ktensor_2way
-    s = """ktensor of shape (2, 2)\nweights=[1. 2.]\nfactor_matrices[0] =\n[[1. 2.]\n [3. 4.]]\nfactor_matrices[1] =\n[[5. 6.]\n [7. 8.]]"""
+    s = """ktensor of shape (2, 2) with order F\nweights=[1. 2.]\nfactor_matrices[0] =\n[[1. 2.]\n [3. 4.]]\nfactor_matrices[1] =\n[[5. 6.]\n [7. 8.]]"""
     assert K0.__str__() == s
diff --git a/tests/test_tenmat.py b/tests/test_tenmat.py
index 987b3eb..150b644 100644
--- a/tests/test_tenmat.py
+++ b/tests/test_tenmat.py
@@ -643,7 +643,7 @@ def test_tenmat__str__(
     # Empty
     tenmatInstance = ttb.tenmat()
     s = ""
-    s += "matrix corresponding to a tensor of shape ()\n"
+    s += "matrix corresponding to a tensor of shape () with order F\n"
     s += "rindices = [  ] (modes of tensor corresponding to rows)\n"
     s += "cindices = [  ] (modes of tensor corresponding to columns)\n"
     s += "data = []\n"
@@ -655,7 +655,7 @@ def test_tenmat__str__(
     s = ""
     s += "matrix corresponding to a tensor of shape "
     s += str(tenmatInstance.tshape)
-    s += "\n"
+    s += " with order F\n"
     s += "rindices = "
     s += "[ " + (", ").join([str(int(d)) for d in tenmatInstance.rindices]) + " ] "
     s += "(modes of tensor corresponding to rows)\n"
@@ -672,7 +672,7 @@ def test_tenmat__str__(
     s = ""
     s += "matrix corresponding to a tensor of shape "
     s += str(tenmatInstance.tshape)
-    s += "\n"
+    s += " with order F\n"
     s += "rindices = "
     s += "[ " + (", ").join([str(int(d)) for d in tenmatInstance.rindices]) + " ] "
     s += "(modes of tensor corresponding to rows)\n"
@@ -690,7 +690,7 @@ def test_tenmat__str__(
     s = ""
     s += "matrix corresponding to a tensor of shape "
     s += str(tenmatInstance.tshape)
-    s += "\n"
+    s += " with order F\n"
     s += "rindices = "
     s += "[ " + (", ").join([str(int(d)) for d in tenmatInstance.rindices]) + " ] "
     s += "(modes of tensor corresponding to rows)\n"
diff --git a/tests/test_tensor.py b/tests/test_tensor.py
index 638db35..762a67c 100644
--- a/tests/test_tensor.py
+++ b/tests/test_tensor.py
@@ -1492,7 +1492,7 @@ def test_tensor__str__(sample_tensor_2way):
     data = np.random.normal(size=(4,))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
+    s += f"tensor of shape {tensorInstance.shape} with order F"
     s += "\ndata"
     s += "[:] =\n"
     s += data.__str__()
@@ -1502,7 +1502,7 @@ def test_tensor__str__(sample_tensor_2way):
     data = np.random.normal(size=(4, 3))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
+    s += f"tensor of shape {tensorInstance.shape} with order F"
     s += "\ndata"
     s += "[:, :] =\n"
     s += data.__str__()
@@ -1512,46 +1512,46 @@ def test_tensor__str__(sample_tensor_2way):
     data = np.random.normal(size=(4, 3, 2))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
-    for i in range(data.shape[0]):
+    s += f"tensor of shape {tensorInstance.shape} with order F"
+    for i in range(data.shape[-1]):
         s += "\ndata"
-        s += "[{}, :, :] =\n".format(i)
-        s += data[i, :, :].__str__()
+        s += "[:, :, {}] =\n".format(i)
+        s += data[:, :, i].__str__()
     assert s == tensorInstance.__str__()
 
     data = np.random.normal(size=(2, 3, 4))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
-    for i in range(data.shape[0]):
+    s += f"tensor of shape {tensorInstance.shape} with order F"
+    for i in range(data.shape[-1]):
         s += "\ndata"
-        s += "[{}, :, :] =\n".format(i)
-        s += data[i, :, :].__str__()
+        s += "[:, :, {}] =\n".format(i)
+        s += data[:, :, i].__str__()
     assert s == tensorInstance.__str__()
 
     # Test 4D
     data = np.random.normal(size=(4, 4, 3, 2))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
-    for i in range(data.shape[0]):
-        for j in range(data.shape[1]):
+    s += f"tensor of shape {tensorInstance.shape} with order F"
+    for i in range(data.shape[-1]):
+        for j in range(data.shape[-2]):
             s += "\ndata"
-            s += "[{}, {}, :, :] =\n".format(j, i)
-            s += data[j, i, :, :].__str__()
+            s += "[:, :, {}, {}] =\n".format(j, i)
+            s += data[:, :, j, i].__str__()
     assert s == tensorInstance.__str__()
 
     # Test 5D
     data = np.random.normal(size=(2, 2, 2, 2, 2))
     tensorInstance = ttb.tensor(data)
     s = ""
-    s += f"tensor of shape {tensorInstance.shape}"
-    for i in range(data.shape[0]):
-        for j in range(data.shape[1]):
-            for k in range(data.shape[2]):
+    s += f"tensor of shape {tensorInstance.shape} with order F"
+    for i in range(data.shape[-1]):
+        for j in range(data.shape[-2]):
+            for k in range(data.shape[-3]):
                 s += "\ndata"
-                s += "[{}, {}, {}, :, :] =\n".format(k, j, i)
-                s += data[k, j, i, :, :].__str__()
+                s += "[:, :, {}, {}, {}] =\n".format(k, j, i)
+                s += data[:, :, k, j, i].__str__()
     assert s == tensorInstance.__str__()