diff --git a/dymos/transcriptions/common/control_group.py b/dymos/transcriptions/common/control_group.py
index 34db4658b..7addfe93f 100644
--- a/dymos/transcriptions/common/control_group.py
+++ b/dymos/transcriptions/common/control_group.py
@@ -113,11 +113,12 @@ def _configure_controls(self):
 
             size = np.prod(shape)
             self.sizes[name] = size
+            sp_eye = sp.eye(size, format='csr')
 
             # The partial of interpolated value wrt the control input values is linear
             # and can be computed as the kronecker product of the interpolation matrix (L)
             # and eye(size).
-            J_val = sp.kron(self.L, sp.eye(size), format='csr')
+            J_val = sp.kron(self.L, sp_eye, format='csr')
             rs, cs, data = sp.find(J_val)
             self.declare_partials(of=self._output_val_names[name],
                                   wrt=self._input_names[name],
@@ -138,14 +139,14 @@ def _configure_controls(self):
             # The partials of the rates and second derivatives are nonlinear but the sparsity
             # pattern is obtained from the kronecker product of the 1st and 2nd differentiation
             # matrices (D and D2) and eye(size).
-            self.rate_jacs[name] = sp.kron(sp.csr_matrix(self.D), sp.eye(size), format='csr')
+            self.rate_jacs[name] = sp.kron(self.D, sp_eye, format='csr')
             rs, cs = self.rate_jacs[name].nonzero()
 
             self.declare_partials(of=self._output_rate_names[name],
                                   wrt=self._input_names[name],
                                   rows=rs, cols=cs)
 
-            self.rate2_jacs[name] = sp.kron(sp.csr_matrix(self.D2), sp.eye(size), format='csr')
+            self.rate2_jacs[name] = sp.kron(self.D2, sp_eye, format='csr')
             rs, cs = self.rate2_jacs[name].nonzero()
 
             self.declare_partials(of=self._output_rate2_names[name],
diff --git a/dymos/transcriptions/pseudospectral/components/state_interp_comp.py b/dymos/transcriptions/pseudospectral/components/state_interp_comp.py
index 4704a3b59..1af8e0e16 100644
--- a/dymos/transcriptions/pseudospectral/components/state_interp_comp.py
+++ b/dymos/transcriptions/pseudospectral/components/state_interp_comp.py
@@ -131,11 +131,10 @@ def configure_io(self):
 
             for key in self.jacs:
                 # Each jacobian matrix has a form that is defined by the Kronecker product
-                # of the interpolation matrix and np.eye(size). Make sure to specify csc format
-                # here to avoid spurious zeros.
+                # of the interpolation matrix eye(size).
                 self.jacs[key][name] = sp.kron(sp.csr_matrix(self.matrices[key]),
-                                               sp.eye(size),
-                                               format='csc')
+                                               sp.eye(size, format='csr'),
+                                               format='csr')
 
             self.sizes[name] = size
 
@@ -157,7 +156,7 @@ def configure_io(self):
                 self.declare_partials(of=self.xc_str[name], wrt=self.xd_str[name],
                                       rows=Ai_rows, cols=Ai_cols, val=data)
 
-                Bi_rows, Bi_cols, _ = sp.find(self.jacs['Bi'][name])
+                Bi_rows, Bi_cols = self.jacs['Bi'][name].nonzero()
                 self.declare_partials(of=self.xc_str[name], wrt=self.fd_str[name],
                                       rows=Bi_rows, cols=Bi_cols)
 
@@ -165,7 +164,7 @@ def configure_io(self):
                 self.declare_partials(of=self.xdotc_str[name], wrt=self.fd_str[name],
                                       rows=Bd_rows, cols=Bd_cols, val=data)
 
-            Ad_rows, Ad_cols, _ = sp.find(self.jacs['Ad'][name])
+            Ad_rows, Ad_cols = self.jacs['Ad'][name].nonzero()
             self.declare_partials(of=self.xdotc_str[name], wrt=self.xd_str[name],
                                   rows=Ad_rows, cols=Ad_cols)
 
@@ -242,7 +241,8 @@ def _compute_partials_radau(self, inputs, partials):
 
             dstau_dt_x_size = np.repeat(dstau_dt, size)[:, np.newaxis]
 
-            partials[xdotc_name, xd_name] = self.jacs['Ad'][name].multiply(dstau_dt_x_size).data
+            dxdotc_dxd = self.jacs['Ad'][name].multiply(dstau_dt_x_size)
+            partials[xdotc_name, xd_name] = dxdotc_dxd.data
 
     def _compute_partials_gauss_lobatto(self, inputs, partials):
         ndn = self.options['grid_data'].subset_num_nodes['state_disc']