diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
index 61da9b957e..ef3746024a 100644
--- a/interface/ceed-basis.c
+++ b/interface/ceed-basis.c
@@ -551,10 +551,9 @@ int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp,
       c1 = c2;
     }
   }
-  //  // Pass to CeedBasisCreateTensorH1
+  // Pass to CeedBasisCreateTensorH1
   ierr = CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d,
-                                 q_ref_1d,
-                                 q_weight_1d, basis); CeedChk(ierr);
+                                 q_ref_1d, q_weight_1d, basis); CeedChk(ierr);
 cleanup:
   ierr2 = CeedFree(&interp_1d); CeedChk(ierr2);
   ierr2 = CeedFree(&grad_1d); CeedChk(ierr2);
diff --git a/interface/ceed-fortran.c b/interface/ceed-fortran.c
index 28366df9d9..941f9e71d0 100644
--- a/interface/ceed-fortran.c
+++ b/interface/ceed-fortran.c
@@ -552,6 +552,15 @@ void fCeedQRFactorization(int *ceed, CeedScalar *mat, CeedScalar *tau, int *m,
   *err = CeedQRFactorization(Ceed_dict[*ceed], mat, tau, *m, *n);
 }
 
+#define fCeedHouseholderApplyQ \
+    FORTRAN_NAME(ceedhouseholderapplyq, CEEDHOUSEHOLDERAPPLYQ)
+void fCeedHouseholderApplyQ(CeedScalar *A, CeedScalar *Q, CeedScalar *tau,
+                            int *t_mode,
+                            int *m, int *n, int *k, int *row, int *col, int *err) {
+  *err = CeedHouseholderApplyQ(A, Q, tau, (CeedTransposeMode)*t_mode, *m, *n, *k,
+                               *row, *col);
+}
+
 #define fCeedSymmetricSchurDecomposition \
     FORTRAN_NAME(ceedsymmetricschurdecomposition, CEEDSYMMETRICSCHURDECOMPOSITION)
 void fCeedSymmetricSchurDecomposition(int *ceed, CeedScalar *mat,
@@ -607,6 +616,26 @@ void fCeedBasisGetInterp1D(int *basis, CeedScalar *interp_1d, int64_t *offset,
   *offset = interp1d_ - interp_1d;
 }
 
+#define fCeedBasisGetGrad1D \
+    FORTRAN_NAME(ceedbasisgetgrad1d, CEEDBASISGETGRAD1D)
+void fCeedBasisGetGrad1D(int *basis, CeedScalar *grad_1d, int64_t *offset,
+                         int *err) {
+  const CeedScalar *grad1d_;
+  CeedBasis basis_ = CeedBasis_dict[*basis];
+  *err = CeedBasisGetGrad1D(basis_, &grad1d_);
+  *offset = grad1d_ - grad_1d;
+}
+
+#define fCeedBasisGetQRef \
+    FORTRAN_NAME(ceedbasisgetqref, CEEDBASISGETQREF)
+void fCeedBasisGetQRef(int *basis, CeedScalar *q_ref, int64_t *offset,
+                       int *err) {
+  const CeedScalar *qref_;
+  CeedBasis basis_ = CeedBasis_dict[*basis];
+  *err = CeedBasisGetQRef(basis_, &qref_);
+  *offset = qref_ - q_ref;
+}
+
 #define fCeedBasisDestroy FORTRAN_NAME(ceedbasisdestroy,CEEDBASISDESTROY)
 void fCeedBasisDestroy(int *basis, int *err) {
   if (*basis == FORTRAN_NULL) return;
diff --git a/python/build_ceed_cffi.py b/python/build_ceed_cffi.py
index 6ffcfe0137..e4c6b99122 100644
--- a/python/build_ceed_cffi.py
+++ b/python/build_ceed_cffi.py
@@ -36,6 +36,7 @@
     header = '\n'.join(lines)
     header = header.split("static inline CeedInt CeedIntPow", 1)[0]
     header += '\nextern int CeedVectorGetState(CeedVector, uint64_t*);'
+    header += '\nextern int CeedElemRestrictionGetELayout(CeedElemRestriction, CeedInt *layout);'
     # Note: cffi cannot handle vargs
     header = re.sub("va_list", "const char *", header)
 ffibuilder.cdef(header)
diff --git a/python/ceed_elemrestriction.py b/python/ceed_elemrestriction.py
index d1ae83fbe4..a08f461394 100644
--- a/python/ceed_elemrestriction.py
+++ b/python/ceed_elemrestriction.py
@@ -131,6 +131,29 @@ def get_multiplicity(self):
         # Return
         return mult
 
+    # Get ElemRestrition Layout
+    def get_layout(self):
+        """Get the element vector layout of an ElemRestriction.
+
+           Returns:
+             layout: Vector containing layout array, stored as [nodes, components, elements].
+                     The data for node i, component j, element k in the element
+                     vector is given by i*layout[0] + j*layout[1] + k*layout[2]."""
+
+        # Create output array
+        layout = np.zeros(3, dtype="int32")
+        array_pointer = ffi.cast(
+            "CeedInt *",
+            layout.__array_interface__['data'][0])
+
+        # libCEED call
+        err_code = lib.CeedElemRestrictionGetELayout(
+            self._pointer[0], array_pointer)
+        self._ceed._check_error(err_code)
+
+        # Return
+        return layout
+
 # ------------------------------------------------------------------------------
 
 
diff --git a/python/tests/output/test_202.out b/python/tests/output/test_202.out
deleted file mode 100644
index 7e9ed6855d..0000000000
--- a/python/tests/output/test_202.out
+++ /dev/null
@@ -1,33 +0,0 @@
-CeedVector length 20
-  10.000000
-  11.000000
-  12.000000
-  13.000000
-  14.000000
-  11.000000
-  12.000000
-  13.000000
-  14.000000
-  15.000000
-  15.000000
-  16.000000
-  17.000000
-  17.000000
-  17.000000
-  16.000000
-  17.000000
-  18.000000
-  18.000000
-  18.000000
-
-CeedVector length 9
-  10.000000
-  22.000000
-  24.000000
-  26.000000
-  28.000000
-  30.000000
-  32.000000
-  34.000000
-  18.000000
-
diff --git a/python/tests/output/test_208.out b/python/tests/output/test_208.out
deleted file mode 100644
index e40943e09b..0000000000
--- a/python/tests/output/test_208.out
+++ /dev/null
@@ -1,23 +0,0 @@
-CeedVector length 10
-  15.000000
-  16.000000
-  17.000000
-  17.000000
-  17.000000
-  16.000000
-  17.000000
-  18.000000
-  18.000000
-  18.000000
-
-CeedVector length 9
-  0.000000
-  0.000000
-  0.000000
-  0.000000
-  0.000000
-  15.000000
-  32.000000
-  34.000000
-  18.000000
-
diff --git a/python/tests/output/test_301.out b/python/tests/output/test_301.out
deleted file mode 100644
index 74a87f6bc9..0000000000
--- a/python/tests/output/test_301.out
+++ /dev/null
@@ -1,15 +0,0 @@
- -2.00000000
- -3.00000000
- -2.00000000
-  0.33333333
- -5.00000000
-  2.00000000
-  0.33333333
-  0.40000000
- -4.00000000
-  0.33333333
- -0.20000000
- -0.50000000
-  1.50000000
-  1.66666667
-  1.60000000
diff --git a/python/tests/output/test_304.out b/python/tests/output/test_304.out
deleted file mode 100644
index 6361fddd1e..0000000000
--- a/python/tests/output/test_304.out
+++ /dev/null
@@ -1,22 +0,0 @@
-Q: 
-  0.69153918
- -0.70710678
-  0.14755868
-  0.00004347
- -0.14721060
-  0.00047835
-  0.69240821
- -0.70632831
-  0.14790715
-  0.00047882
- -0.69066919
- -0.70788369
- -0.69153892
- -0.70710646
- -0.14755804
- -0.00100064
-lambda: 
-  0.13487964
-  0.23325338
-  0.86513545
-  1.16666509
diff --git a/python/tests/test-2-elemrestriction.py b/python/tests/test-2-elemrestriction.py
index 95de43d5bf..d3ccbe1dc6 100644
--- a/python/tests/test-2-elemrestriction.py
+++ b/python/tests/test-2-elemrestriction.py
@@ -30,26 +30,26 @@
 def test_200(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
-    x = ceed.Vector(ne + 1)
-    a = np.arange(10, 10 + ne + 1, dtype="float64")
+    x = ceed.Vector(num_elem + 1)
+    a = np.arange(10, 10 + num_elem + 1, dtype="float64")
     x.set_array(a, cmode=libceed.USE_POINTER)
 
-    ind = np.zeros(2 * ne, dtype="int32")
-    for i in range(ne):
+    ind = np.zeros(2 * num_elem, dtype="int32")
+    for i in range(num_elem):
         ind[2 * i + 0] = i
         ind[2 * i + 1] = i + 1
-    r = ceed.ElemRestriction(ne, 2, 1, 1, ne + 1, ind,
+    r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind,
                              cmode=libceed.USE_POINTER)
 
-    y = ceed.Vector(2 * ne)
+    y = ceed.Vector(2 * num_elem)
     y.set_value(0)
 
     r.apply(x, y)
 
     with y.array_read() as y_array:
-        for i in range(2 * ne):
+        for i in range(2 * num_elem):
             assert 10 + (i + 1) // 2 == y_array[i]
 
 # -------------------------------------------------------------------------------
@@ -60,22 +60,22 @@ def test_200(ceed_resource):
 def test_201(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
-    x = ceed.Vector(2 * ne)
-    a = np.arange(10, 10 + 2 * ne, dtype="float64")
+    x = ceed.Vector(2 * num_elem)
+    a = np.arange(10, 10 + 2 * num_elem, dtype="float64")
     x.set_array(a, cmode=libceed.USE_POINTER)
 
     strides = np.array([1, 2, 2], dtype="int32")
-    r = ceed.StridedElemRestriction(ne, 2, 1, 2 * ne, strides)
+    r = ceed.StridedElemRestriction(num_elem, 2, 1, 2 * num_elem, strides)
 
-    y = ceed.Vector(2 * ne)
+    y = ceed.Vector(2 * num_elem)
     y.set_value(0)
 
     r.apply(x, y)
 
     with y.array_read() as y_array:
-        for i in range(2 * ne):
+        for i in range(2 * num_elem):
             assert 10 + i == y_array[i]
 
 # -------------------------------------------------------------------------------
@@ -86,71 +86,93 @@ def test_201(ceed_resource):
 def test_202(ceed_resource, capsys):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 8
-    blksize = 5
+    num_elem = 8
+    elem_size = 2
+    num_blk = 2
+    blk_size = 5
 
-    x = ceed.Vector(ne + 1)
-    a = np.arange(10, 10 + ne + 1, dtype="float64")
-    x.set_array(a, cmode=libceed.USE_POINTER)
+    x = ceed.Vector(num_elem + 1)
+    a = np.arange(10, 10 + num_elem + 1, dtype="float64")
+    x.set_array(a, cmode=libceed.COPY_VALUES)
 
-    ind = np.zeros(2 * ne, dtype="int32")
-    for i in range(ne):
+    ind = np.zeros(2 * num_elem, dtype="int32")
+    for i in range(num_elem):
         ind[2 * i + 0] = i
         ind[2 * i + 1] = i + 1
-    r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne + 1, ind,
-                                    cmode=libceed.USE_POINTER)
+    r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1,
+                                    num_elem + 1, ind, cmode=libceed.COPY_VALUES)
 
-    y = ceed.Vector(2 * blksize * 2)
+    y = ceed.Vector(num_blk * blk_size * elem_size)
     y.set_value(0)
 
+    # NoTranspose
     r.apply(x, y)
-
-    print(y)
+    layout = r.get_layout()
+    with y.array_read() as y_array:
+        for i in range(elem_size):
+            for j in range(1):
+                for k in range(num_elem):
+                    block = int(k / blk_size)
+                    elem = k % blk_size
+                    indx = (i * blk_size + elem) * \
+                        layout[0] + j * blk_size * layout[1] + \
+                        block * blk_size * layout[2]
+                assert y_array[indx] == a[ind[k * elem_size + i]]
 
     x.set_value(0)
     r.T.apply(y, x)
-    print(x)
-
-    stdout, stderr, ref_stdout = check.output(capsys)
-    assert not stderr
-    assert stdout == ref_stdout
+    with x.array_read() as x_array:
+        for i in range(num_elem + 1):
+            assert x_array[i] == (10 + i) * (2.0 if i >
+                                             0 and i < num_elem else 1.0)
 
 # -------------------------------------------------------------------------------
 # Test creation, use, and destruction of a blocked element restriction
 # -------------------------------------------------------------------------------
 
 
-def test_208(ceed_resource, capsys):
+def test_208(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 8
-    blksize = 5
+    num_elem = 8
+    elem_size = 2
+    num_blk = 2
+    blk_size = 5
 
-    x = ceed.Vector(ne + 1)
-    a = np.arange(10, 10 + ne + 1, dtype="float64")
-    x.set_array(a, cmode=libceed.USE_POINTER)
+    x = ceed.Vector(num_elem + 1)
+    a = np.arange(10, 10 + num_elem + 1, dtype="float64")
+    x.set_array(a, cmode=libceed.COPY_VALUES)
 
-    ind = np.zeros(2 * ne, dtype="int32")
-    for i in range(ne):
+    ind = np.zeros(2 * num_elem, dtype="int32")
+    for i in range(num_elem):
         ind[2 * i + 0] = i
         ind[2 * i + 1] = i + 1
-    r = ceed.BlockedElemRestriction(ne, 2, blksize, 1, 1, ne + 1, ind,
-                                    cmode=libceed.USE_POINTER)
+    r = ceed.BlockedElemRestriction(num_elem, elem_size, blk_size, 1, 1,
+                                    num_elem + 1, ind, cmode=libceed.COPY_VALUES)
 
-    y = ceed.Vector(blksize * 2)
+    y = ceed.Vector(blk_size * elem_size)
     y.set_value(0)
 
+    # NoTranspose
     r.apply_block(1, x, y)
-
-    print(y)
+    layout = r.get_layout()
+    with y.array_read() as y_array:
+        for i in range(elem_size):
+            for j in range(1):
+                for k in range(blk_size, num_elem):
+                    block = int(k / blk_size)
+                    elem = k % blk_size
+                    indx = (i * blk_size + elem) * layout[0] + j * blk_size * \
+                        layout[1] + block * blk_size * \
+                        layout[2] - blk_size * elem_size
+                assert y_array[indx] == a[ind[k * elem_size + i]]
 
     x.set_value(0)
     r.T.apply_block(1, y, x)
-    print(x)
-
-    stdout, stderr, ref_stdout = check.output(capsys)
-    assert not stderr
-    assert stdout == ref_stdout
+    with x.array_read() as x_array:
+        for i in range(blk_size, num_elem + 1):
+            assert x_array[i] == (10 + i) * (2.0 if i >
+                                             blk_size and i < num_elem else 1.0)
 
 # -------------------------------------------------------------------------------
 # Test getting the multiplicity of the indices in an element restriction
@@ -160,22 +182,22 @@ def test_208(ceed_resource, capsys):
 def test_209(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
-    ind = np.zeros(4 * ne, dtype="int32")
-    for i in range(ne):
+    ind = np.zeros(4 * num_elem, dtype="int32")
+    for i in range(num_elem):
         ind[4 * i + 0] = i * 3 + 0
         ind[4 * i + 1] = i * 3 + 1
         ind[4 * i + 2] = i * 3 + 2
         ind[4 * i + 3] = i * 3 + 3
-    r = ceed.ElemRestriction(ne, 4, 1, 1, 3 * ne + 1, ind,
+    r = ceed.ElemRestriction(num_elem, 4, 1, 1, 3 * num_elem + 1, ind,
                              cmode=libceed.USE_POINTER)
 
     mult = r.get_multiplicity()
 
     with mult.array_read() as mult_array:
-        for i in range(3 * ne + 1):
-            val = 1 + (1 if (i > 0 and i < 3 * ne and i % 3 == 0) else 0)
+        for i in range(3 * num_elem + 1):
+            val = 1 + (1 if (i > 0 and i < 3 * num_elem and i % 3 == 0) else 0)
             assert val == mult_array[i]
 
 # -------------------------------------------------------------------------------
@@ -186,13 +208,13 @@ def test_209(ceed_resource):
 def test_210(ceed_resource, capsys):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
-    ind = np.zeros(2 * ne, dtype="int32")
-    for i in range(ne):
+    ind = np.zeros(2 * num_elem, dtype="int32")
+    for i in range(num_elem):
         ind[2 * i + 0] = i + 0
         ind[2 * i + 1] = i + 1
-    r = ceed.ElemRestriction(ne, 2, 1, 1, ne + 1, ind,
+    r = ceed.ElemRestriction(num_elem, 2, 1, 1, num_elem + 1, ind,
                              cmode=libceed.USE_POINTER)
 
     print(r)
@@ -209,10 +231,10 @@ def test_210(ceed_resource, capsys):
 def test_211(ceed_resource, capsys):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
     strides = np.array([1, 2, 2], dtype="int32")
-    r = ceed.StridedElemRestriction(ne, 2, 1, ne + 1, strides)
+    r = ceed.StridedElemRestriction(num_elem, 2, 1, num_elem + 1, strides)
 
     print(r)
 
@@ -228,10 +250,11 @@ def test_211(ceed_resource, capsys):
 def test_212(ceed_resource, capsys):
     ceed = libceed.Ceed(ceed_resource)
 
-    ne = 3
+    num_elem = 3
 
     strides = np.array([1, 2, 2], dtype="int32")
-    r = ceed.BlockedStridedElemRestriction(ne, 2, 2, 1, ne + 1, strides)
+    r = ceed.BlockedStridedElemRestriction(
+        num_elem, 2, 2, 1, num_elem + 1, strides)
 
     print(r)
 
diff --git a/python/tests/test-3-basis.py b/python/tests/test-3-basis.py
index a9366f74ce..0b83aeefe0 100644
--- a/python/tests/test-3-basis.py
+++ b/python/tests/test-3-basis.py
@@ -71,34 +71,32 @@ def test_300(ceed_resource, capsys):
 # -------------------------------------------------------------------------------
 
 
-def test_301(ceed_resource, capsys):
+def test_301(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
+    m = 4
+    n = 3
+    a = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype="float64")
     qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype="float64")
     tau = np.empty(3, dtype="float64")
 
-    qr, tau = ceed.qr_factorization(qr, tau, 4, 3)
+    qr, tau = ceed.qr_factorization(qr, tau, m, n)
+    np_qr, np_tau = np.linalg.qr(a.reshape(m, n), mode="raw")
 
-    for i in range(len(qr)):
-        if qr[i] <= TOL and qr[i] >= -TOL:
-            qr[i] = 0
-        print("%12.8f" % qr[i])
+    for i in range(n):
+        assert tau[i] == np_tau[i]
 
-    for i in range(len(tau)):
-        if tau[i] <= TOL and tau[i] >= -TOL:
-            tau[i] = 0
-        print("%12.8f" % tau[i])
-
-    stdout, stderr, ref_stdout = check.output(capsys)
-    assert not stderr
-    assert stdout == ref_stdout
+    qr = qr.reshape(m, n)
+    for i in range(m):
+        for j in range(n):
+            assert round(qr[i, j] - np_qr[j, i], 10) == 0
 
 # -------------------------------------------------------------------------------
 # Test Symmetric Schur Decomposition
 # -------------------------------------------------------------------------------
 
 
-def test_304(ceed_resource, capsys):
+def test_304(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
     A = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333,
@@ -120,7 +118,7 @@ def test_304(ceed_resource, capsys):
 # -------------------------------------------------------------------------------
 
 
-def test_305(ceed_resource, capsys):
+def test_305(ceed_resource):
     ceed = libceed.Ceed(ceed_resource)
 
     M = np.array([0.2, 0.0745355993, -0.0745355993, 0.0333333333,
diff --git a/tests/output/t202-elemrestriction-f.out b/tests/output/t202-elemrestriction-f.out
deleted file mode 100644
index 55e23727e1..0000000000
--- a/tests/output/t202-elemrestriction-f.out
+++ /dev/null
@@ -1,31 +0,0 @@
-CeedVector length 20
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   18.00000000
-   18.00000000
-CeedVector length 9
-   10.00000000
-   22.00000000
-   24.00000000
-   26.00000000
-   28.00000000
-   30.00000000
-   32.00000000
-   34.00000000
-   18.00000000
diff --git a/tests/output/t202-elemrestriction.out b/tests/output/t202-elemrestriction.out
deleted file mode 100644
index 55e23727e1..0000000000
--- a/tests/output/t202-elemrestriction.out
+++ /dev/null
@@ -1,31 +0,0 @@
-CeedVector length 20
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   18.00000000
-   18.00000000
-CeedVector length 9
-   10.00000000
-   22.00000000
-   24.00000000
-   26.00000000
-   28.00000000
-   30.00000000
-   32.00000000
-   34.00000000
-   18.00000000
diff --git a/tests/output/t203-elemrestriction.out b/tests/output/t203-elemrestriction.out
deleted file mode 100644
index cfea41922d..0000000000
--- a/tests/output/t203-elemrestriction.out
+++ /dev/null
@@ -1,117 +0,0 @@
-CeedVector length 27
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   20.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   25.00000000
-   26.00000000
-   27.00000000
-   28.00000000
-   30.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   35.00000000
-   36.00000000
-   37.00000000
-   38.00000000
-CeedVector length 60
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   20.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   25.00000000
-   30.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   35.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   18.00000000
-   18.00000000
-   25.00000000
-   26.00000000
-   27.00000000
-   27.00000000
-   27.00000000
-   26.00000000
-   27.00000000
-   28.00000000
-   28.00000000
-   28.00000000
-   35.00000000
-   36.00000000
-   37.00000000
-   37.00000000
-   37.00000000
-   36.00000000
-   37.00000000
-   38.00000000
-   38.00000000
-   38.00000000
-CeedVector length 27
-   10.00000000
-   22.00000000
-   24.00000000
-   26.00000000
-   28.00000000
-   30.00000000
-   32.00000000
-   34.00000000
-   18.00000000
-   20.00000000
-   42.00000000
-   44.00000000
-   46.00000000
-   48.00000000
-   50.00000000
-   52.00000000
-   54.00000000
-   28.00000000
-   30.00000000
-   62.00000000
-   64.00000000
-   66.00000000
-   68.00000000
-   70.00000000
-   72.00000000
-   74.00000000
-   38.00000000
diff --git a/tests/output/t208-elemrestriction-f.out b/tests/output/t208-elemrestriction-f.out
deleted file mode 100644
index 7337918474..0000000000
--- a/tests/output/t208-elemrestriction-f.out
+++ /dev/null
@@ -1,21 +0,0 @@
-CeedVector length 10
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   18.00000000
-   18.00000000
-CeedVector length 9
-    0.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-   15.00000000
-   32.00000000
-   34.00000000
-   18.00000000
diff --git a/tests/output/t208-elemrestriction.out b/tests/output/t208-elemrestriction.out
deleted file mode 100644
index a4871fd5db..0000000000
--- a/tests/output/t208-elemrestriction.out
+++ /dev/null
@@ -1,21 +0,0 @@
-CeedVector length 10
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-CeedVector length 9
-    0.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-    0.00000000
-   15.00000000
-   32.00000000
-   17.00000000
-    0.00000000
diff --git a/tests/output/t213-elemrestriction.out b/tests/output/t213-elemrestriction.out
deleted file mode 100644
index cfea41922d..0000000000
--- a/tests/output/t213-elemrestriction.out
+++ /dev/null
@@ -1,117 +0,0 @@
-CeedVector length 27
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   20.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   25.00000000
-   26.00000000
-   27.00000000
-   28.00000000
-   30.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   35.00000000
-   36.00000000
-   37.00000000
-   38.00000000
-CeedVector length 60
-   10.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   11.00000000
-   12.00000000
-   13.00000000
-   14.00000000
-   15.00000000
-   20.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   21.00000000
-   22.00000000
-   23.00000000
-   24.00000000
-   25.00000000
-   30.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   31.00000000
-   32.00000000
-   33.00000000
-   34.00000000
-   35.00000000
-   15.00000000
-   16.00000000
-   17.00000000
-   17.00000000
-   17.00000000
-   16.00000000
-   17.00000000
-   18.00000000
-   18.00000000
-   18.00000000
-   25.00000000
-   26.00000000
-   27.00000000
-   27.00000000
-   27.00000000
-   26.00000000
-   27.00000000
-   28.00000000
-   28.00000000
-   28.00000000
-   35.00000000
-   36.00000000
-   37.00000000
-   37.00000000
-   37.00000000
-   36.00000000
-   37.00000000
-   38.00000000
-   38.00000000
-   38.00000000
-CeedVector length 27
-   10.00000000
-   22.00000000
-   24.00000000
-   26.00000000
-   28.00000000
-   30.00000000
-   32.00000000
-   34.00000000
-   18.00000000
-   20.00000000
-   42.00000000
-   44.00000000
-   46.00000000
-   48.00000000
-   50.00000000
-   52.00000000
-   54.00000000
-   28.00000000
-   30.00000000
-   62.00000000
-   64.00000000
-   66.00000000
-   68.00000000
-   70.00000000
-   72.00000000
-   74.00000000
-   38.00000000
diff --git a/tests/output/t301-basis-f.out b/tests/output/t301-basis-f.out
deleted file mode 100644
index 74a87f6bc9..0000000000
--- a/tests/output/t301-basis-f.out
+++ /dev/null
@@ -1,15 +0,0 @@
- -2.00000000
- -3.00000000
- -2.00000000
-  0.33333333
- -5.00000000
-  2.00000000
-  0.33333333
-  0.40000000
- -4.00000000
-  0.33333333
- -0.20000000
- -0.50000000
-  1.50000000
-  1.66666667
-  1.60000000
diff --git a/tests/output/t301-basis.out b/tests/output/t301-basis.out
deleted file mode 100644
index 74a87f6bc9..0000000000
--- a/tests/output/t301-basis.out
+++ /dev/null
@@ -1,15 +0,0 @@
- -2.00000000
- -3.00000000
- -2.00000000
-  0.33333333
- -5.00000000
-  2.00000000
-  0.33333333
-  0.40000000
- -4.00000000
-  0.33333333
- -0.20000000
- -0.50000000
-  1.50000000
-  1.66666667
-  1.60000000
diff --git a/tests/output/t302-basis-f.out b/tests/output/t302-basis-f.out
deleted file mode 100644
index cd6dbc37ad..0000000000
--- a/tests/output/t302-basis-f.out
+++ /dev/null
@@ -1,14 +0,0 @@
-collograd[0]: -3.00000000  4.04508497 -1.54508497  0.50000000
-collograd[1]: -0.80901699  0.00000000  1.11803399 -0.30901699
-collograd[2]:  0.30901699 -1.11803399  0.00000000  0.80901699
-collograd[3]: -0.50000000  1.54508497 -4.04508497  3.00000000
-collograd[0]: -3.33200024  4.86015442 -2.10878235  0.58062817
-collograd[1]: -0.75755761 -0.38441439  1.47067023 -0.32869822
-collograd[2]:  0.32869822 -1.47067023  0.38441439  0.75755761
-collograd[3]: -0.58062817  2.10878235 -4.86015442  3.33200024
-collograd[0]: -3.46592188  1.89221184  2.88856261 -0.64549009 -1.79434637  1.12498389
-collograd[1]: -1.80001352  0.46288791  1.30407559  0.28647572 -0.39349201  0.14006630
-collograd[2]: -0.10371066 -0.80402118 -0.31095294  0.88495577  0.82906498 -0.49533597
-collograd[3]:  0.49533597 -0.82906498 -0.88495577  0.31095294  0.80402118  0.10371066
-collograd[4]: -0.14006630  0.39349201 -0.28647572 -1.30407559 -0.46288791  1.80001352
-collograd[5]: -1.12498389  1.79434637  0.64549009 -2.88856261 -1.89221184  3.46592188
diff --git a/tests/output/t302-basis.out b/tests/output/t302-basis.out
deleted file mode 100644
index 00337185df..0000000000
--- a/tests/output/t302-basis.out
+++ /dev/null
@@ -1,14 +0,0 @@
- collograd1d[0]:	 -3.00000000	  4.04508497	 -1.54508497	  0.50000000
- collograd1d[1]:	 -0.80901699	  0.00000000	  1.11803399	 -0.30901699
- collograd1d[2]:	  0.30901699	 -1.11803399	  0.00000000	  0.80901699
- collograd1d[3]:	 -0.50000000	  1.54508497	 -4.04508497	  3.00000000
- collograd1d[0]:	 -3.33200024	  4.86015442	 -2.10878235	  0.58062817
- collograd1d[1]:	 -0.75755761	 -0.38441439	  1.47067023	 -0.32869822
- collograd1d[2]:	  0.32869822	 -1.47067023	  0.38441439	  0.75755761
- collograd1d[3]:	 -0.58062817	  2.10878235	 -4.86015442	  3.33200024
- collograd1d[0]:	 -3.46592188	  1.89221184	  2.88856261	 -0.64549009	 -1.79434637	  1.12498389
- collograd1d[1]:	 -1.80001352	  0.46288791	  1.30407559	  0.28647572	 -0.39349201	  0.14006630
- collograd1d[2]:	 -0.10371066	 -0.80402118	 -0.31095294	  0.88495577	  0.82906498	 -0.49533597
- collograd1d[3]:	  0.49533597	 -0.82906498	 -0.88495577	  0.31095294	  0.80402118	  0.10371066
- collograd1d[4]:	 -0.14006630	  0.39349201	 -0.28647572	 -1.30407559	 -0.46288791	  1.80001352
- collograd1d[5]:	 -1.12498389	  1.79434637	  0.64549009	 -2.88856261	 -1.89221184	  3.46592188
diff --git a/tests/t201-elemrestriction-f.f90 b/tests/t201-elemrestriction-f.f90
index fbf7571519..306427d797 100644
--- a/tests/t201-elemrestriction-f.f90
+++ b/tests/t201-elemrestriction-f.f90
@@ -35,8 +35,8 @@ program test
       strides=[1,2,2]
       call ceedelemrestrictioncreatestrided(ceed,ne,2,1,2*ne,strides,r,err)
 
-      call ceedvectorcreate(ceed,2*ne,y,err);
-      call ceedvectorsetvalue(y,0.d0,err);
+      call ceedvectorcreate(ceed,2*ne,y,err)
+      call ceedvectorsetvalue(y,0.d0,err)
       call ceedelemrestrictionapply(r,ceed_notranspose,x,y,&
      & ceed_request_immediate,err)
 
diff --git a/tests/t201-elemrestriction.c b/tests/t201-elemrestriction.c
index 30a882743c..05fcb27649 100644
--- a/tests/t201-elemrestriction.c
+++ b/tests/t201-elemrestriction.c
@@ -28,13 +28,13 @@ int main(int argc, char **argv) {
 
   CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
   CeedElemRestrictionGetELayout(r, &layout);
-  for (CeedInt i=0; i<2; i++)      // Node
-    for (CeedInt j=0; j<1; j++)    // Component
+  for (CeedInt i=0; i<2; i++)            // Node
+    for (CeedInt j=0; j<1; j++)          // Component
       for (CeedInt k=0; k<num_elem; k++) // Element
         if (yy[i*layout[0] + j*layout[1] + k*layout[2]] !=
             a[i*strides[0] + j*strides[1] + k*strides[2]])
           // LCOV_EXCL_START
-          printf("Error in restricted array y[%d][%d][%d] = %f",
+          printf("Error in restricted array y[%d][%d][%d] = %f\n",
                  i, j, k, (double)yy[i*strides[0] + j*strides[1] + j*strides[2]]);
   // LCOV_EXCL_STOP
   CeedVectorRestoreArrayRead(y, &yy);
diff --git a/tests/t202-elemrestriction-f.f90 b/tests/t202-elemrestriction-f.f90
index 36e1852c5f..ea199b5016 100644
--- a/tests/t202-elemrestriction-f.f90
+++ b/tests/t202-elemrestriction-f.f90
@@ -5,16 +5,20 @@ program test
 
       integer ceed,err
       integer x,y
-      integer i,r
+      integer r
+      integer i,j,k
 
-      integer ne
-      parameter(ne=8)
-      integer blksize
-      parameter(blksize=5)
+      integer ne,elemsize,nb,blksize
+      parameter(ne=8,elemsize=2,nb=2,blksize=5)
+      integer ind(elemsize*ne)
+      integer layout(3)
+      integer blk,elem,indx
 
-      integer*4 ind(2*ne)
       real*8 a(ne+1)
-      integer*8 aoffset
+      real*8 yy(nb*blksize*elemsize)
+      real*8 xx(ne+1)
+      real*8 diff
+      integer*8 aoffset,xoffset,yoffset
 
       character arg*32
 
@@ -34,29 +38,56 @@ program test
         ind(2*i-1)=i-1
         ind(2*i  )=i
       enddo
-
-      call ceedelemrestrictioncreateblocked(ceed,ne,2,blksize,1,1,ne+1,&
+      call ceedelemrestrictioncreateblocked(ceed,ne,elemsize,blksize,1,1,ne+1,&
      & ceed_mem_host,ceed_use_pointer,ind,r,err)
 
-      call ceedvectorcreate(ceed,2*blksize*2,y,err);
-      call ceedvectorsetvalue(y,0.d0,err);
+      call ceedvectorcreate(ceed,nb*blksize*elemsize,y,err);
+      call ceedvectorsetvalue(y,0.d0,err)
 
-!    No Transpose
+!     NoTranspose
       call ceedelemrestrictionapply(r,ceed_notranspose,x,y,&
      & ceed_request_immediate,err)
-      call ceedvectorview(y,err)
-
-!    Transpose
-      call ceedvectorgetarray(x,ceed_mem_host,a,aoffset,err)
-      do i=1,ne+1
-        a(aoffset+i)=0.0
+      call ceedvectorgetarrayread(y,ceed_mem_host,yy,yoffset,err)
+      call ceedelemrestrictiongetelayout(r,layout,err)
+      do i=0,elemsize-1
+        do j=0,0
+          do k=0,ne-1
+            blk = k/blksize
+            elem = mod(k,blksize)
+            indx = (i*blksize+elem)*layout(1)+j*layout(2)*blksize+blk*layout(3)*blksize
+            diff=yy(yoffset+indx+1)
+            diff=diff-a(ind(k*elemsize+i+1)+1)
+            if (abs(diff) > 1.0D-15) then
+! LCOV_EXCL_START
+             write(*,*) 'Error in restricted array y(',i,')(',j,')(',k,')=',&
+     &         yy(yoffset+indx+1)
+! LCOV_EXCL_STOP
+            endif
+          enddo
+        enddo
       enddo
-      call ceedvectorrestorearray(x,a,aoffset,err)
-      
+      call ceedvectorrestorearrayread(y,yy,yoffset,err)
+
+!     Transpose
+      call ceedvectorsetvalue(x,0.d0,err)
       call ceedelemrestrictionapply(r,ceed_transpose,y,x,&
      & ceed_request_immediate,err)
-
-      call ceedvectorview(x,err)
+      call ceedvectorgetarrayread(x,ceed_mem_host,xx,xoffset,err)
+      do i=1,ne+1
+        diff=xx(xoffset+i)
+        if (i > 1 .and. i < ne+1) then
+          diff=diff-2*(10+i-1)
+        else
+          diff=diff-(10+i-1)
+        endif
+        if (abs(diff) > 1.0D-15) then
+! LCOV_EXCL_START
+         write(*,*) 'Error in restricted array x(',i,')=',&
+     &         xx(xoffset+i)
+! LCOV_EXCL_STOP
+        endif
+      enddo
+      call ceedvectorrestorearrayread(x,xx,xoffset,err)
 
       call ceedvectordestroy(x,err)
       call ceedvectordestroy(y,err)
diff --git a/tests/t202-elemrestriction.c b/tests/t202-elemrestriction.c
index 8fea5331fb..0a11b126cb 100644
--- a/tests/t202-elemrestriction.c
+++ b/tests/t202-elemrestriction.c
@@ -2,14 +2,19 @@
 /// Test creation, use, and destruction of a blocked element restriction
 /// \test Test creation, use, and destruction of a blocked element restriction
 #include <ceed.h>
+#include <ceed/backend.h>
 
 int main(int argc, char **argv) {
   Ceed ceed;
   CeedVector x, y;
   CeedInt num_elem = 8;
+  CeedInt elem_size = 2;
+  CeedInt num_blk = 2;
   CeedInt blk_size = 5;
-  CeedInt ind[2*num_elem];
-  CeedScalar a[num_elem+1];
+  CeedInt ind[elem_size*num_elem];
+  CeedScalar a[num_elem + 1];
+  const CeedScalar *xx, *yy;
+  CeedInt layout[3];
   CeedElemRestriction r;
 
   CeedInit(argv[1], &ceed);
@@ -23,22 +28,42 @@ int main(int argc, char **argv) {
     ind[2*i+0] = i;
     ind[2*i+1] = i+1;
   }
-  CeedElemRestrictionCreateBlocked(ceed, num_elem, 2, blk_size, 1, 1, num_elem+1,
-                                   CEED_MEM_HOST, CEED_USE_POINTER, ind, &r);
-  CeedVectorCreate(ceed, 2*blk_size*2, &y);
+  CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, 1, 1,
+                                   num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER,
+                                   ind, &r);
+  CeedVectorCreate(ceed, num_blk*blk_size*elem_size, &y);
   CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(y, "%12.8f", stdout);
+  CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
+  CeedElemRestrictionGetELayout(r, &layout);
+  for (CeedInt i=0; i<elem_size; i++)      // Node
+    for (CeedInt j=0; j<1; j++)            // Component
+      for (CeedInt k=0; k<num_elem; k++) { // Element
+        CeedInt block = k / blk_size;
+        CeedInt elem = k % blk_size;
+        CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
+                        block*layout[2]*blk_size;
+        if (yy[index] != a[ind[k*elem_size + i]])
+          // LCOV_EXCL_START
+          printf("Error in restricted array y[%d][%d][%d] = %f\n",
+                 i, j, k, (double)yy[index]);
+        // LCOV_EXCL_STOP
+      }
+  CeedVectorRestoreArrayRead(y, &yy);
 
   // Transpose
-  CeedVectorGetArray(x, CEED_MEM_HOST, (CeedScalar **)&a);
-  for (CeedInt i=0; i<num_elem+1; i++)
-    a[i] = 0;
-  CeedVectorRestoreArray(x, (CeedScalar **)&a);
+  CeedVectorSetValue(x, 0);
   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(x, "%12.8f", stdout);
+  CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
+  for (CeedInt i=0; i<num_elem+1; i++)
+    if (xx[i] != (10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
+      // LCOV_EXCL_START
+      printf("Error in restricted array x[%d] = %f\n",
+             i, (double)xx[i]);
+  // LCOV_EXCL_STOP
+  CeedVectorRestoreArrayRead(x, &xx);
 
   CeedVectorDestroy(&x);
   CeedVectorDestroy(&y);
diff --git a/tests/t203-elemrestriction.c b/tests/t203-elemrestriction.c
index c9e0609666..848ea8d2a1 100644
--- a/tests/t203-elemrestriction.c
+++ b/tests/t203-elemrestriction.c
@@ -2,46 +2,75 @@
 /// Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector
 /// \test Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector
 #include <ceed.h>
+#include <ceed/backend.h>
 
 int main(int argc, char **argv) {
   Ceed ceed;
   CeedVector x, y;
   CeedInt num_elem = 8;
+  CeedInt elem_size = 2;
+  CeedInt num_blk = 2;
   CeedInt blk_size = 5;
   CeedInt num_comp = 3;
-  CeedInt ind[2*num_elem];
-  CeedScalar a[num_comp*(num_elem+1)];
+  CeedInt ind[elem_size*num_elem];
+  CeedScalar a[num_comp*(num_elem + 1)];
+  const CeedScalar *xx, *yy;
+  CeedInt layout[3];
   CeedElemRestriction r;
 
   CeedInit(argv[1], &ceed);
 
-  CeedVectorCreate(ceed, (num_elem+1)*num_comp, &x);
-  for (CeedInt i=0; i<(num_elem+1); i++) {
+  CeedVectorCreate(ceed, num_comp*(num_elem+1), &x);
+  for (CeedInt i=0; i<num_elem+1; i++) {
     a[i+0*(num_elem+1)] = 10 + i;
     a[i+1*(num_elem+1)] = 20 + i;
     a[i+2*(num_elem+1)] = 30 + i;
   }
   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
-  CeedVectorView(x, "%12.8f", stdout);
+
   for (CeedInt i=0; i<num_elem; i++) {
     ind[2*i+0] = i;
     ind[2*i+1] = i+1;
   }
-  CeedElemRestrictionCreateBlocked(ceed, num_elem, 2, blk_size, num_comp,
-                                   num_elem+1,
-                                   num_comp*(num_elem+1), CEED_MEM_HOST,
+  CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp,
+                                   num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
                                    CEED_USE_POINTER, ind, &r);
-  CeedVectorCreate(ceed, 2*blk_size*2*num_comp, &y);
+  CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
   CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(y, "%12.8f", stdout);
+  CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
+  CeedElemRestrictionGetELayout(r, &layout);
+  for (CeedInt i=0; i<elem_size; i++)      // Node
+    for (CeedInt j=0; j<num_comp; j++)     // Component
+      for (CeedInt k=0; k<num_elem; k++) { // Element
+        CeedInt block = k / blk_size;
+        CeedInt elem = k % blk_size;
+        CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
+                        block*layout[2]*blk_size;
+        if (yy[index] != a[ind[k*elem_size + i]+j*(num_elem+1)])
+          // LCOV_EXCL_START
+          printf("Error in restricted array y[%d][%d][%d] = %f\n",
+                 i, j, k, (double)yy[index]);
+        // LCOV_EXCL_STOP
+      }
+  CeedVectorRestoreArrayRead(y, &yy);
 
   // Transpose
-  CeedVectorSetValue(x, 0.0);
+  CeedVectorSetValue(x, 0);
   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(x, "%12.8f", stdout);
+  CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
+  for (CeedInt i=0; i<num_elem+1; i++) {
+    for (CeedInt j=0; j<num_comp; j++) {
+      if (xx[i+j*(num_elem+1)] != ((j+1)*10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
+        // LCOV_EXCL_START
+        printf("Error in restricted array x[%d][%d] = %f\n",
+               j, i, (double)xx[i+j*(num_elem+1)]);
+      // LCOV_EXCL_STOP
+    }
+  }
+  CeedVectorRestoreArrayRead(x, &xx);
 
   CeedVectorDestroy(&x);
   CeedVectorDestroy(&y);
diff --git a/tests/t208-elemrestriction-f.f90 b/tests/t208-elemrestriction-f.f90
index eb0a2e1347..a703e4b3c0 100644
--- a/tests/t208-elemrestriction-f.f90
+++ b/tests/t208-elemrestriction-f.f90
@@ -5,16 +5,20 @@ program test
 
       integer ceed,err
       integer x,y
-      integer i,r
+      integer r
+      integer i,j,k
 
-      integer ne
-      parameter(ne=8)
-      integer blksize
-      parameter(blksize=5)
+      integer ne,elemsize,nb,blksize
+      parameter(ne=8,elemsize=2,nb=2,blksize=5)
+      integer ind(elemsize*ne)
+      integer layout(3)
+      integer blk,elem,indx
 
-      integer*4 ind(2*ne)
       real*8 a(ne+1)
-      integer*8 aoffset
+      real*8 yy(nb*blksize*elemsize)
+      real*8 xx(ne+1)
+      real*8 diff
+      integer*8 aoffset,xoffset,yoffset
 
       character arg*32
 
@@ -34,29 +38,57 @@ program test
         ind(2*i-1)=i-1
         ind(2*i  )=i
       enddo
-
-      call ceedelemrestrictioncreateblocked(ceed,ne,2,blksize,1,1,ne+1,&
+      call ceedelemrestrictioncreateblocked(ceed,ne,elemsize,blksize,1,1,ne+1,&
      & ceed_mem_host,ceed_use_pointer,ind,r,err)
-
-      call ceedvectorcreate(ceed,blksize*2,y,err);
-      call ceedvectorsetvalue(y,0.d0,err);
-
-!    No Transpose
-      call ceedelemrestrictionapplyblock(r,1,ceed_notranspose,&
-     & x,y,ceed_request_immediate,err)
-      call ceedvectorview(y,err)
-
-!    Transpose
-      call ceedvectorgetarray(x,ceed_mem_host,a,aoffset,err)
-      do i=1,ne+1
-        a(aoffset+i)=0.0
+  
+      call ceedvectorcreate(ceed,blksize*elemsize,y,err);
+      call ceedvectorsetvalue(y,0.d0,err)
+
+!     NoTranspose
+      call ceedelemrestrictionapplyblock(r,1,ceed_notranspose,x,y,&
+     & ceed_request_immediate,err)
+      call ceedvectorgetarrayread(y,ceed_mem_host,yy,yoffset,err)
+      call ceedelemrestrictiongetelayout(r,layout,err)
+      do i=0,elemsize-1
+        do j=0,0
+          do k=blksize,ne-1
+            blk = k/blksize
+            elem = mod(k,blksize)
+            indx = (i*blksize+elem)*layout(1)+j*layout(2)*blksize+blk*layout(3)*blksize-&
+     & blksize*elemsize
+            diff=yy(yoffset+indx+1)
+            diff=diff-a(ind(k*elemsize+i+1)+1)
+            if (abs(diff) > 1.0D-15) then
+! LCOV_EXCL_START
+             write(*,*) 'Error in restricted array y(',i,')(',j,')(',k,')=',&
+     &         yy(yoffset+indx+1)
+! LCOV_EXCL_STOP
+            endif
+          enddo
+        enddo
       enddo
-      call ceedvectorrestorearray(x,a,aoffset,err)
-      
-      call ceedelemrestrictionapplyblock(r,1,ceed_transpose,&
-     & y,x,ceed_request_immediate,err)
-
-      call ceedvectorview(x,err)
+      call ceedvectorrestorearrayread(y,yy,yoffset,err)
+
+!     Transpose
+      call ceedvectorsetvalue(x,0.d0,err)
+      call ceedelemrestrictionapplyblock(r,1,ceed_transpose,y,x,&
+     & ceed_request_immediate,err)
+      call ceedvectorgetarrayread(x,ceed_mem_host,xx,xoffset,err)
+      do i=blksize+1,ne+1
+        diff=xx(xoffset+i)
+        if (i > blksize+1 .and. i < ne+1) then
+          diff=diff-2*(10+i-1)
+        else
+          diff=diff-(10+i-1)
+        endif
+        if (abs(diff) > 1.0D-15) then
+! LCOV_EXCL_START
+         write(*,*) 'Error in restricted array x(',i,')=',&
+     &         xx(xoffset+i)
+! LCOV_EXCL_STOP
+        endif
+      enddo
+      call ceedvectorrestorearrayread(x,xx,xoffset,err)
 
       call ceedvectordestroy(x,err)
       call ceedvectordestroy(y,err)
diff --git a/tests/t208-elemrestriction.c b/tests/t208-elemrestriction.c
index cd9dee718d..343516decd 100644
--- a/tests/t208-elemrestriction.c
+++ b/tests/t208-elemrestriction.c
@@ -2,17 +2,20 @@
 /// Test creation, use, and destruction of a blocked element restriction
 /// \test Test creation, use, and destruction of a blocked element restriction
 #include <ceed.h>
+#include <ceed/backend.h>
 
 int main(int argc, char **argv) {
   Ceed ceed;
   CeedVector x, y;
   CeedInt num_elem = 8;
-  CeedInt blk_size = 5;
   CeedInt elem_size = 2;
+  CeedInt num_blk = 2;
+  CeedInt blk_size = 5;
   CeedInt ind[elem_size*num_elem];
-  CeedScalar a[num_elem+1];
+  CeedScalar a[num_elem + 1];
+  const CeedScalar *xx, *yy;
+  CeedInt layout[3];
   CeedElemRestriction r;
-  CeedScalar *y_array;
 
   CeedInit(argv[1], &ceed);
 
@@ -22,40 +25,47 @@ int main(int argc, char **argv) {
   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
 
   for (CeedInt i=0; i<num_elem; i++) {
-    for (CeedInt k=0; k<elem_size; k++) {
-      ind[elem_size*i+k] = i+k;
-    }
+    ind[2*i+0] = i;
+    ind[2*i+1] = i+1;
   }
-
   CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, 1, 1,
-                                   num_elem+1,
-                                   CEED_MEM_HOST, CEED_USE_POINTER, ind, &r);
-
+                                   num_elem+1, CEED_MEM_HOST, CEED_USE_POINTER,
+                                   ind, &r);
   CeedVectorCreate(ceed, blk_size*elem_size, &y);
   CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApplyBlock(r, 1, CEED_NOTRANSPOSE, x, y,
                                 CEED_REQUEST_IMMEDIATE);
-
-  // Zero padded entries
-  CeedVectorGetArray(y, CEED_MEM_HOST, &y_array);
-  for (CeedInt i = (elem_size*num_elem - blk_size*elem_size);
-       i < blk_size*elem_size;
-       ++i) {
-    y_array[i] = 0;
-  }
-  CeedVectorRestoreArray(y, &y_array);
-  CeedVectorView(y, "%12.8f", stdout);
+  CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
+  CeedElemRestrictionGetELayout(r, &layout);
+  for (CeedInt i=0; i<elem_size; i++)             // Node
+    for (CeedInt j=0; j<1; j++)                   // Component
+      for (CeedInt k=blk_size; k<num_elem; k++) { // Element
+        CeedInt block = k / blk_size;
+        CeedInt elem = k % blk_size;
+        CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
+                        block*layout[2]*blk_size - blk_size*elem_size;
+        if (yy[index] != a[ind[k*elem_size + i]])
+          // LCOV_EXCL_START
+          printf("Error in restricted array y[%d][%d][%d] = %f\n",
+                 i, j, k, (double)yy[index]);
+        // LCOV_EXCL_STOP
+      }
+  CeedVectorRestoreArrayRead(y, &yy);
 
   // Transpose
-  CeedVectorGetArray(x, CEED_MEM_HOST, (CeedScalar **)&a);
-  for (CeedInt i=0; i<num_elem+1; i++)
-    a[i] = 0;
-  CeedVectorRestoreArray(x, (CeedScalar **)&a);
+  CeedVectorSetValue(x, 0);
   CeedElemRestrictionApplyBlock(r, 1, CEED_TRANSPOSE, y, x,
                                 CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(x, "%12.8f", stdout);
+  CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
+  for (CeedInt i=blk_size; i<num_elem+1; i++)
+    if (xx[i] != (10+i)*(i > blk_size && i < num_elem ? 2.0 : 1.0))
+      // LCOV_EXCL_START
+      printf("Error in restricted array x[%d] = %f\n",
+             i, (double)xx[i]);
+  // LCOV_EXCL_STOP
+  CeedVectorRestoreArrayRead(x, &xx);
 
   CeedVectorDestroy(&x);
   CeedVectorDestroy(&y);
diff --git a/tests/t213-elemrestriction.c b/tests/t213-elemrestriction.c
index 9b32734301..17bf192c9a 100644
--- a/tests/t213-elemrestriction.c
+++ b/tests/t213-elemrestriction.c
@@ -1,51 +1,80 @@
 /// @file
 /// Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer
 /// \test Test creation, use, and destruction of a blocked element restriction with multiple components in the lvector, and libCEED owns ind pointer
-#include <stdlib.h>
 #include <ceed.h>
+#include <ceed/backend.h>
+#include <stdlib.h>
+#include <string.h>
 
 int main(int argc, char **argv) {
   Ceed ceed;
   CeedVector x, y;
   CeedInt num_elem = 8;
+  CeedInt elem_size = 2;
+  CeedInt num_blk = 2;
   CeedInt blk_size = 5;
   CeedInt num_comp = 3;
-  CeedInt *ind = malloc(sizeof(CeedInt)*2*num_elem);
-  CeedScalar a[num_comp*(num_elem+1)];
+  CeedInt ind[elem_size*num_elem];
+  CeedInt *ceed_ind = malloc(sizeof(CeedInt)*elem_size*num_elem);
+  CeedScalar a[num_comp*(num_elem + 1)];
+  const CeedScalar *xx, *yy;
+  CeedInt layout[3];
   CeedElemRestriction r;
 
   CeedInit(argv[1], &ceed);
 
-  CeedVectorCreate(ceed, (num_elem+1)*num_comp, &x);
-  for (CeedInt i=0; i<(num_elem+1); i++) {
+  CeedVectorCreate(ceed, num_comp*(num_elem+1), &x);
+  for (CeedInt i=0; i<num_elem+1; i++) {
     a[i+0*(num_elem+1)] = 10 + i;
     a[i+1*(num_elem+1)] = 20 + i;
     a[i+2*(num_elem+1)] = 30 + i;
   }
   CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, a);
-  CeedVectorView(x, "%12.8f", stdout);
+
   for (CeedInt i=0; i<num_elem; i++) {
     ind[2*i+0] = i;
     ind[2*i+1] = i+1;
   }
-  CeedElemRestrictionCreateBlocked(ceed, num_elem, 2, blk_size, num_comp,
-                                   num_elem+1,
-                                   num_comp*(num_elem+1), CEED_MEM_HOST,
-                                   CEED_OWN_POINTER, ind, &r);
-  CeedVectorCreate(ceed, 2*blk_size*2*num_comp, &y);
+  memcpy(ceed_ind, ind, sizeof(CeedInt)*elem_size*num_elem);
+  CeedElemRestrictionCreateBlocked(ceed, num_elem, elem_size, blk_size, num_comp,
+                                   num_elem+1, num_comp*(num_elem+1), CEED_MEM_HOST,
+                                   CEED_OWN_POINTER, ceed_ind, &r);
+  CeedVectorCreate(ceed, num_comp*num_blk*blk_size*elem_size, &y);
   CeedVectorSetValue(y, 0); // Allocates array
 
   // NoTranspose
   CeedElemRestrictionApply(r, CEED_NOTRANSPOSE, x, y, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(y, "%12.8f", stdout);
+  CeedVectorGetArrayRead(y, CEED_MEM_HOST, &yy);
+  CeedElemRestrictionGetELayout(r, &layout);
+  for (CeedInt i=0; i<elem_size; i++)      // Node
+    for (CeedInt j=0; j<num_comp; j++)     // Component
+      for (CeedInt k=0; k<num_elem; k++) { // Element
+        CeedInt block = k / blk_size;
+        CeedInt elem = k % blk_size;
+        CeedInt index = (i*blk_size+elem)*layout[0] + j*layout[1]*blk_size +
+                        block*layout[2]*blk_size;
+        if (yy[index] != a[ind[k*elem_size + i]+j*(num_elem+1)])
+          // LCOV_EXCL_START
+          printf("Error in restricted array y[%d][%d][%d] = %f\n",
+                 i, j, k, (double)yy[index]);
+        // LCOV_EXCL_STOP
+      }
+  CeedVectorRestoreArrayRead(y, &yy);
 
   // Transpose
-  CeedVectorGetArray(x, CEED_MEM_HOST, (CeedScalar **)&a);
-  for (CeedInt i=0; i<(num_elem+1)*num_comp; i++)
-    a[i] = 0;
-  CeedVectorRestoreArray(x, (CeedScalar **)&a);
+  CeedVectorSetValue(x, 0);
   CeedElemRestrictionApply(r, CEED_TRANSPOSE, y, x, CEED_REQUEST_IMMEDIATE);
-  CeedVectorView(x, "%12.8f", stdout);
+  CeedVectorGetArrayRead(x, CEED_MEM_HOST, &xx);
+  for (CeedInt i=0; i<num_elem+1; i++) {
+    for (CeedInt j=0; j<num_comp; j++) {
+      if (xx[i+j*(num_elem+1)] != ((j+1)*10+i)*(i > 0 && i < num_elem ? 2.0 : 1.0))
+        // LCOV_EXCL_START
+        printf("Error in restricted array x[%d][%d] = %f\n",
+               j, i, (double)xx[i+j*(num_elem+1)]);
+      // LCOV_EXCL_STOP
+    }
+  }
+  CeedVectorRestoreArrayRead(x, &xx);
 
   CeedVectorDestroy(&x);
   CeedVectorDestroy(&y);
diff --git a/tests/t301-basis-f.f90 b/tests/t301-basis-f.f90
index d010ef3789..3496fbfa20 100644
--- a/tests/t301-basis-f.f90
+++ b/tests/t301-basis-f.f90
@@ -3,32 +3,33 @@ program test
       implicit none
       include 'ceed/fortran.h'
 
-      integer ceed,err,i
-      real*8 qr(12), tau(3)
+      integer ceed,err,i,j
+      real*8 a(12),qr(12),a_qr(12),tau(3)
 
       character arg*32
 
+      A = (/ 1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0 /)
       qr = (/ 1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0 /)
+      a_qr = (/ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /)
 
       call getarg(1,arg)
 
       call ceedinit(trim(arg)//char(0),ceed,err)
-      call ceedqrfactorization(ceed,qr,tau,4,3,err);
-      do i=1,12
-        if (abs(qr(i))<1.0D-14) then
-! LCOV_EXCL_START
-          qr(i) = 0
-! LCOV_EXCL_STOP
-        endif
-        write(*,'(A,F12.8)') '',qr(i)
-      enddo
+      call ceedqrfactorization(ceed,qr,tau,4,3,err)
       do i=1,3
-        if (abs(tau(i))<1.0D-14) then
+        do j=i,3
+          a_qr((i-1)*3+j)=qr((i-1)*3+j)
+        enddo
+      enddo
+      call ceedhouseholderapplyq(a_qr,qr,tau,ceed_notranspose,4,3,3,3,1,err)
+
+      do i=1,12
+        if (abs(a(i)-a_qr(i))>1.0D-14) then
 ! LCOV_EXCL_START
-          tau(i) = 0
+        write(*,*) 'Error in QR factorization a_qr(',i,') = ',a_qr(i),&
+    &    ' != a(',i,') = ',a(i)
 ! LCOV_EXCL_STOP
-        endif
-        write(*,'(A,F12.8)') '',tau(i)
+      endif
       enddo
 
       call ceeddestroy(ceed,err)
diff --git a/tests/t301-basis.c b/tests/t301-basis.c
index b8de71ee55..53334d4522 100644
--- a/tests/t301-basis.c
+++ b/tests/t301-basis.c
@@ -2,23 +2,31 @@
 /// Test QR Factorization
 /// \test Test QR Factorization
 #include <ceed.h>
+#include <ceed/backend.h>
+#include <math.h>
 
 int main(int argc, char **argv) {
   Ceed ceed;
-  CeedScalar qr[12] = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
-  CeedScalar tau[3];
+  CeedScalar A[12]    = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
+  CeedScalar qr[12]   = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
+  CeedScalar A_qr[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  CeedScalar tau[4];
 
   CeedInit(argv[1], &ceed);
 
   CeedQRFactorization(ceed, qr, tau, 4, 3);
-  for (int i=0; i<12; i++) {
-    if (qr[i] <= 1E-14 && qr[i] >= -1E-14) qr[i] = 0;
-    fprintf(stdout, "%12.8f\n", qr[i]);
-  }
-  for (int i=0; i<3; i++) {
-    if (tau[i] <= 1E-14 && tau[i] >= -1E-14) tau[i] = 0;
-    fprintf(stdout, "%12.8f\n", tau[i]);
-  }
+  for (CeedInt i=0; i<3; i++)
+    for (CeedInt j=i; j<3; j++)
+      A_qr[i*3+j] = qr[i*3+j];
+  CeedHouseholderApplyQ(A_qr, qr, tau, CEED_NOTRANSPOSE, 4, 3, 3, 3, 1);
+
+  for (CeedInt i=0; i<12; i++)
+    if (fabs(A_qr[i] - A[i]) > 10*CEED_EPSILON)
+      // LCOV_EXCL_START
+      printf("Error in QR factorization A_qr[%d] = %f != A[%d] = %f\n",
+             i, A_qr[i], i, A[i]);
+  // LCOV_EXCL_STOP
+
   CeedDestroy(&ceed);
   return 0;
 }
diff --git a/tests/t302-basis-f.f90 b/tests/t302-basis-f.f90
index c114753445..4206afcb7f 100644
--- a/tests/t302-basis-f.f90
+++ b/tests/t302-basis-f.f90
@@ -4,8 +4,12 @@ program test
       include 'ceed/fortran.h'
 
       integer ceed,err,i,j
-      integer b
-      real*8 collograd1d(16), collograd1d2(36)
+      integer b,p
+      parameter(p=4)
+      real*8 collograd1d(36),x2(6)
+      real*8 grad1d(16),qref(6)
+      integer*8 gradoffset,qoffset
+      real*8 sum
 
       character arg*32
 
@@ -14,56 +18,69 @@ program test
       call ceedinit(trim(arg)//char(0),ceed,err)
 
 !     Already collocated, GetCollocatedGrad will return grad1d
-      call ceedbasiscreatetensorh1lagrange(ceed,1,1,4,4,ceed_gauss_lobatto,b,&
+      call ceedbasiscreatetensorh1lagrange(ceed,1,1,p,p,ceed_gauss_lobatto,b,&
      & err)
       call ceedbasisgetcollocatedgrad(b,collograd1d,err)
-      do i=1,16
-        if (abs(collograd1d(i))<1.0D-14) then
-          collograd1d(i) = 0
-        endif
-      enddo
-      do i=0,3
-        write(*,'(A,I1,A,F12.8,F12.8,F12.8,F12.8,F12.8,F12.8)')&
-     &   'collograd[',i,']:',(collograd1d(j+4*i),j=1,4)
-      call flush(6)
+      call ceedbasisgetgrad1d(b,grad1d,gradoffset,err)
+      do i=0,p-1
+        do j=1,p
+          if (abs(collograd1d(j+p*i)-grad1d(j+p*i+gradoffset))>1.0D-13) then
+! LCOV_EXCL_START
+            write(*,*) 'Error in collocated gradient ',collograd1d(j+p*i),' != ',&
+     &       grad1d(j+p*i+gradoffset)
+! LCOV_EXCL_STOP
+          endif
+        enddo
       enddo
       call ceedbasisdestroy(b,err)
 
 !     Q = P, not already collocated
-      call ceedbasiscreatetensorh1lagrange(ceed,1,1,4,4,ceed_gauss,b,err)
+      call ceedbasiscreatetensorh1lagrange(ceed,1,1,p,p,ceed_gauss,b,err)
       call ceedbasisgetcollocatedgrad(b,collograd1d,err)
-      do i=1,16
-        if (abs(collograd1d(i))<1.0D-14) then
+
+      call ceedbasisgetqref(b,qref,qoffset,err)
+      do i=1,p
+        x2(i)=qref(i+qoffset)*qref(i+qoffset)
+      enddo
+
+      do i=0,p-1
+        sum=0
+        do j=1,p
+            sum=sum+collograd1d(j+p*i)*x2(j)
+        enddo
+        if (abs(sum-2*qref(i+1+qoffset))>1.0D-13) then
 ! LCOV_EXCL_START
-          collograd1d(i) = 0
+            write(*,*) 'Error in collocated gradient ',sum,' != ',&
+            &       2*qref(i+1+qoffset)
 ! LCOV_EXCL_STOP
         endif
       enddo
-      do i=0,3
-        write(*,'(A,I1,A,F12.8,F12.8,F12.8,F12.8,F12.8,F12.8)')&
-     &   'collograd[',i,']:',(collograd1d(j+4*i),j=1,4)
-      call flush(6)
-      enddo
       call ceedbasisdestroy(b,err)
 
 !     Q = P + 2, not already collocated
-      call ceedbasiscreatetensorh1lagrange(ceed,1,1,4,6,ceed_gauss,b,err)
-      call ceedbasisgetcollocatedgrad(b,collograd1d2,err)
-      do i=1,36
-        if (abs(collograd1d2(i))<1.0D-14) then
+      call ceedbasiscreatetensorh1lagrange(ceed,1,1,p,p+2,ceed_gauss,b,err)
+      call ceedbasisgetcollocatedgrad(b,collograd1d,err)
+
+      call ceedbasisgetqref(b,qref,qoffset,err)
+      do i=1,p+2
+        x2(i)=qref(i+qoffset)*qref(i+qoffset)
+      enddo
+
+      do i=0,p+1
+        sum=0
+        do j=1,p+2
+            sum=sum+collograd1d(j+(p+2)*i)*x2(j)
+        enddo
+        if (abs(sum-2*qref(i+1+qoffset))>1.0D-13) then
 ! LCOV_EXCL_START
-          collograd1d2(i) = 0
+            write(*,*) 'Error in collocated gradient ',sum,' != ',&
+     &       2*qref(i+1+qoffset)
 ! LCOV_EXCL_STOP
         endif
       enddo
-      do i=0,5
-        write(*,'(A,I1,A,F12.8,F12.8,F12.8,F12.8,F12.8,F12.8)')&
-     &   'collograd[',i,']:',(collograd1d2(j+6*i),j=1,6)
-      call flush(6)
-      enddo
       call ceedbasisdestroy(b,err)
 
       call ceeddestroy(ceed,err)
 
       end
-!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
\ No newline at end of file
diff --git a/tests/t302-basis.c b/tests/t302-basis.c
index c52a30e6c8..7463dbacf2 100644
--- a/tests/t302-basis.c
+++ b/tests/t302-basis.c
@@ -7,51 +7,65 @@
 
 int main(int argc, char **argv) {
   Ceed ceed;
-  CeedInt P=4, Q=4;
-  CeedScalar collo_grad_1d[Q*Q], collo_grad_1d_2[(P+2)*(P+2)];
+  CeedInt P = 4;
+  CeedScalar collo_grad_1d[(P+2)*(P+2)], x_2[P+2];
+  const CeedScalar *grad_1d, *q_ref;
+  CeedScalar sum = 0.0;
   CeedBasis b;
 
   CeedInit(argv[1], &ceed);
 
   // Already collocated, GetCollocatedGrad will return grad_1d
-  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, Q, CEED_GAUSS_LOBATTO, &b);
+  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS_LOBATTO, &b);
   CeedBasisGetCollocatedGrad(b, collo_grad_1d);
+  CeedBasisGetGrad(b, &grad_1d);
 
-  for (int i=0; i<Q; i++) {
-    fprintf(stdout, "%12s[%d]:", "collograd1d", i);
-    for (int j=0; j<Q; j++) {
-      if (fabs(collo_grad_1d[j+Q*i]) <= 1E-14) collo_grad_1d[j+Q*i] = 0;
-      fprintf(stdout, "\t% 12.8f", collo_grad_1d[j+Q*i]);
-    }
-    fprintf(stdout, "\n");
-  }
+  for (CeedInt i=0; i<P; i++)
+    for (CeedInt j=0; j<P; j++)
+      if (fabs(collo_grad_1d[j+P*i] - grad_1d[j+P*i]) > 100*CEED_EPSILON)
+        // LCOV_EXCL_START
+        printf("Error in collocated gradient %f != %f\n", collo_grad_1d[j+P*i],
+               grad_1d[j+P*i]);
+  // LCOV_EXCL_START
   CeedBasisDestroy(&b);
 
   // Q = P, not already collocated
-  CeedBasisCreateTensorH1Lagrange(ceed, 1,  1, P, Q, CEED_GAUSS, &b);
+  CeedBasisCreateTensorH1Lagrange(ceed, 1,  1, P, P, CEED_GAUSS, &b);
   CeedBasisGetCollocatedGrad(b, collo_grad_1d);
 
-  for (int i=0; i<Q; i++) {
-    fprintf(stdout, "%12s[%d]:", "collograd1d", i);
-    for (int j=0; j<Q; j++) {
-      if (fabs(collo_grad_1d[j+Q*i]) <= 1E-14) collo_grad_1d[j+Q*i] = 0;
-      fprintf(stdout, "\t% 12.8f", collo_grad_1d[j+Q*i]);
-    }
-    fprintf(stdout, "\n");
+  CeedBasisGetQRef(b, &q_ref);
+  for (CeedInt i=0; i<P; i++)
+    x_2[i] = q_ref[i]*q_ref[i];
+
+  // Verify collo_grad * x^2 = 2x for quadrature points
+  for (CeedInt i=0; i<P; i++) {
+    sum = 0.0;
+    for (CeedInt j=0; j<P; j++)
+      sum += collo_grad_1d[j+P*i]*x_2[j];
+    if (fabs(sum - 2*q_ref[i]) > 100*CEED_EPSILON)
+      // LCOV_EXCL_START
+      printf("Error in collocated gradient %f != %f\n", sum, 2*q_ref[i]);
+    // LCOV_EXCL_START
   }
   CeedBasisDestroy(&b);
 
   // Q = P + 2, not already collocated
   CeedBasisCreateTensorH1Lagrange(ceed, 1,  1, P, P+2, CEED_GAUSS, &b);
-  CeedBasisGetCollocatedGrad(b, collo_grad_1d_2);
-
-  for (int i=0; i<P+2; i++) {
-    fprintf(stdout, "%12s[%d]:", "collograd1d", i);
-    for (int j=0; j<P+2; j++) {
-      if (fabs(collo_grad_1d_2[j+(P+2)*i]) <= 1E-14) collo_grad_1d_2[j+(P+2)*i] = 0;
-      fprintf(stdout, "\t% 12.8f", collo_grad_1d_2[j+(P+2)*i]);
-    }
-    fprintf(stdout, "\n");
+  CeedBasisGetCollocatedGrad(b, collo_grad_1d);
+
+  CeedBasisGetQRef(b, &q_ref);
+  for (CeedInt i=0; i<P+2; i++)
+    x_2[i] = q_ref[i]*q_ref[i];
+
+  // Verify collo_grad * x^2 = 2x for quadrature points
+  for (CeedInt i=0; i<P+2; i++) {
+    sum = 0.0;
+    for (CeedInt j=0; j<P+2; j++)
+      sum += collo_grad_1d[j+(P+2)*i]*x_2[j];
+    if (fabs(sum - 2*q_ref[i]) > 100*CEED_EPSILON)
+      // LCOV_EXCL_START
+      printf("Error in collocated gradient %f != %f\n", sum, 2*q_ref[i]);
+    // LCOV_EXCL_START
   }
   CeedBasisDestroy(&b);