diff --git a/arraycontext/context.py b/arraycontext/context.py
index aa3054d94c320d0252d1e2705c00e68558698f1b..0c9232d8c40de87500ef19a55004ae5b06cce6cc 100644
--- a/arraycontext/context.py
+++ b/arraycontext/context.py
@@ -151,13 +151,15 @@ class ArrayContext(ABC):
     .. attribute:: np
 
          Provides access to a namespace that serves as a work-alike to
-         :mod:`numpy`.  The actual level of functionality provided is up to the
+         :mod:`numpy`. The actual level of functionality provided is up to the
          individual array context implementation, however the functions and
          objects available under this namespace must not behave differently
          from :mod:`numpy`.
 
          As a baseline, special functions available through :mod:`loopy`
          (e.g. ``sin``, ``exp``) are accessible through this interface.
+         A full list of implemented functionality is given in
+         :ref:`numpy-coverage`.
 
          Callables accessible through this namespace vectorize over object
          arrays, including :class:`arraycontext.ArrayContainer`\ s.
diff --git a/arraycontext/impl/pyopencl/fake_numpy.py b/arraycontext/impl/pyopencl/fake_numpy.py
index 514193ff4e448dde3c441899411f4f7ee6e34124..bae5b34fae1bac2939e2396150686cfaf66e80b0 100644
--- a/arraycontext/impl/pyopencl/fake_numpy.py
+++ b/arraycontext/impl/pyopencl/fake_numpy.py
@@ -58,137 +58,37 @@ class PyOpenCLFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
     def _get_fake_numpy_linalg_namespace(self):
         return _PyOpenCLFakeNumpyLinalgNamespace(self._array_context)
 
-    # {{{ comparisons
+    # NOTE: the order of these follows the order in numpy docs
+    # NOTE: when adding a function here, also add it to `array_context.rst` docs!
 
-    # FIXME: This should be documentation, not a comment.
-    # These are here mainly because some arrays may choose to interpret
-    # equality comparison as a binary predicate of structural identity,
-    # i.e. more like "are you two equal", and not like numpy semantics.
-    # These operations provide access to numpy-style comparisons in that
-    # case.
-
-    def equal(self, x, y):
-        return rec_multimap_array_container(operator.eq, x, y)
-
-    def not_equal(self, x, y):
-        return rec_multimap_array_container(operator.ne, x, y)
-
-    def greater(self, x, y):
-        return rec_multimap_array_container(operator.gt, x, y)
-
-    def greater_equal(self, x, y):
-        return rec_multimap_array_container(operator.ge, x, y)
-
-    def less(self, x, y):
-        return rec_multimap_array_container(operator.lt, x, y)
-
-    def less_equal(self, x, y):
-        return rec_multimap_array_container(operator.le, x, y)
-
-    # }}}
+    # {{{ array creation routines
 
     def ones_like(self, ary):
-        def _ones_like(subary):
+        return self.full_like(ary, 1)
+
+    def full_like(self, ary, fill_value):
+        def _full_like(subary):
             ones = self._array_context.empty_like(subary)
-            ones.fill(1)
+            ones.fill(fill_value)
             return ones
 
-        return self._new_like(ary, _ones_like)
-
-    def maximum(self, x, y):
-        return rec_multimap_array_container(
-                partial(cl_array.maximum, queue=self._array_context.queue),
-                x, y)
-
-    def minimum(self, x, y):
-        return rec_multimap_array_container(
-                partial(cl_array.minimum, queue=self._array_context.queue),
-                x, y)
-
-    def where(self, criterion, then, else_):
-        def where_inner(inner_crit, inner_then, inner_else):
-            if isinstance(inner_crit, bool):
-                return inner_then if inner_crit else inner_else
-            return cl_array.if_positive(inner_crit != 0, inner_then, inner_else,
-                    queue=self._array_context.queue)
-
-        return rec_multimap_array_container(where_inner, criterion, then, else_)
-
-    def sum(self, a, axis=None, dtype=None):
-
-        if isinstance(axis, int):
-            axis = axis,
-
-        def _rec_sum(ary):
-            if axis not in [None, tuple(range(ary.ndim))]:
-                raise NotImplementedError(f"Sum over '{axis}' axes not supported.")
-
-            return cl_array.sum(ary, dtype=dtype, queue=self._array_context.queue)
-
-        result = rec_map_reduce_array_container(sum, _rec_sum, a)
-
-        if not self._array_context._force_device_scalars:
-            result = result.get()[()]
-        return result
-
-    def min(self, a, axis=None):
-        queue = self._array_context.queue
-
-        if isinstance(axis, int):
-            axis = axis,
-
-        def _rec_min(ary):
-            if axis not in [None, tuple(range(ary.ndim))]:
-                raise NotImplementedError(f"Min. over '{axis}' axes not supported.")
-            return cl_array.min(ary, queue=queue)
-
-        result = rec_map_reduce_array_container(
-                partial(reduce, partial(cl_array.minimum, queue=queue)),
-                _rec_min,
-                a)
-
-        if not self._array_context._force_device_scalars:
-            result = result.get()[()]
-        return result
-
-    def max(self, a, axis=None):
-        queue = self._array_context.queue
-
-        if isinstance(axis, int):
-            axis = axis,
+        return self._new_like(ary, _full_like)
 
-        def _rec_max(ary):
-            if axis not in [None, tuple(range(ary.ndim))]:
-                raise NotImplementedError(f"Max. over '{axis}' axes not supported.")
-            return cl_array.max(ary, queue=queue)
+    def copy(self, ary):
+        def _copy(subary):
+            return subary.copy(queue=self._array_context.queue)
 
-        result = rec_map_reduce_array_container(
-                partial(reduce, partial(cl_array.maximum, queue=queue)),
-                _rec_max,
-                a)
+        return self._new_like(ary, _copy)
 
-        if not self._array_context._force_device_scalars:
-            result = result.get()[()]
-        return result
+    # }}}
 
-    def stack(self, arrays, axis=0):
-        return rec_multimap_array_container(
-                lambda *args: cl_array.stack(arrays=args, axis=axis,
-                    queue=self._array_context.queue),
-                *arrays)
+    # {{{ array manipulation routines
 
     def reshape(self, a, newshape, order="C"):
         return rec_map_array_container(
                 lambda ary: ary.reshape(newshape, order=order),
                 a)
 
-    def concatenate(self, arrays, axis=0):
-        return cl_array.concatenate(
-            arrays, axis,
-            self._array_context.queue,
-            self._array_context.allocator
-        )
-
     def ravel(self, a, order="C"):
         def _rec_ravel(a):
             if order in "FC":
@@ -211,6 +111,23 @@ class PyOpenCLFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
 
         return rec_map_array_container(_rec_ravel, a)
 
+    def concatenate(self, arrays, axis=0):
+        return cl_array.concatenate(
+            arrays, axis,
+            self._array_context.queue,
+            self._array_context.allocator
+        )
+
+    def stack(self, arrays, axis=0):
+        return rec_multimap_array_container(
+                lambda *args: cl_array.stack(arrays=args, axis=axis,
+                    queue=self._array_context.queue),
+                *arrays)
+
+    # }}}
+
+    # {{{ linear algebra
+
     def vdot(self, x, y, dtype=None):
         result = rec_multimap_reduce_array_container(
                 sum,
@@ -221,22 +138,26 @@ class PyOpenCLFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
             result = result.get()[()]
         return result
 
-    def any(self, a):
+    # }}}
+
+    # {{{ logic functions
+
+    def all(self, a):
         queue = self._array_context.queue
         result = rec_map_reduce_array_container(
-                partial(reduce, partial(cl_array.maximum, queue=queue)),
-                lambda subary: subary.any(queue=queue),
+                partial(reduce, partial(cl_array.minimum, queue=queue)),
+                lambda subary: subary.all(queue=queue),
                 a)
 
         if not self._array_context._force_device_scalars:
             result = result.get()[()]
         return result
 
-    def all(self, a):
+    def any(self, a):
         queue = self._array_context.queue
         result = rec_map_reduce_array_container(
-                partial(reduce, partial(cl_array.minimum, queue=queue)),
-                lambda subary: subary.all(queue=queue),
+                partial(reduce, partial(cl_array.maximum, queue=queue)),
+                lambda subary: subary.any(queue=queue),
                 a)
 
         if not self._array_context._force_device_scalars:
@@ -273,6 +194,124 @@ class PyOpenCLFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
 
         return result
 
+    # FIXME: This should be documentation, not a comment.
+    # These are here mainly because some arrays may choose to interpret
+    # equality comparison as a binary predicate of structural identity,
+    # i.e. more like "are you two equal", and not like numpy semantics.
+    # These operations provide access to numpy-style comparisons in that
+    # case.
+
+    def greater(self, x, y):
+        return rec_multimap_array_container(operator.gt, x, y)
+
+    def greater_equal(self, x, y):
+        return rec_multimap_array_container(operator.ge, x, y)
+
+    def less(self, x, y):
+        return rec_multimap_array_container(operator.lt, x, y)
+
+    def less_equal(self, x, y):
+        return rec_multimap_array_container(operator.le, x, y)
+
+    def equal(self, x, y):
+        return rec_multimap_array_container(operator.eq, x, y)
+
+    def not_equal(self, x, y):
+        return rec_multimap_array_container(operator.ne, x, y)
+
+    # }}}
+
+    # {{{ mathematical functions
+
+    def sum(self, a, axis=None, dtype=None):
+        if isinstance(axis, int):
+            axis = axis,
+
+        def _rec_sum(ary):
+            if axis not in [None, tuple(range(ary.ndim))]:
+                raise NotImplementedError(f"Sum over '{axis}' axes not supported.")
+
+            return cl_array.sum(ary, dtype=dtype, queue=self._array_context.queue)
+
+        result = rec_map_reduce_array_container(sum, _rec_sum, a)
+
+        if not self._array_context._force_device_scalars:
+            result = result.get()[()]
+        return result
+
+    def maximum(self, x, y):
+        return rec_multimap_array_container(
+                partial(cl_array.maximum, queue=self._array_context.queue),
+                x, y)
+
+    def amax(self, a, axis=None):
+        queue = self._array_context.queue
+
+        if isinstance(axis, int):
+            axis = axis,
+
+        def _rec_max(ary):
+            if axis not in [None, tuple(range(ary.ndim))]:
+                raise NotImplementedError(f"Max. over '{axis}' axes not supported.")
+            return cl_array.max(ary, queue=queue)
+
+        result = rec_map_reduce_array_container(
+                partial(reduce, partial(cl_array.maximum, queue=queue)),
+                _rec_max,
+                a)
+
+        if not self._array_context._force_device_scalars:
+            result = result.get()[()]
+        return result
+
+    max = amax
+
+    def minimum(self, x, y):
+        return rec_multimap_array_container(
+                partial(cl_array.minimum, queue=self._array_context.queue),
+                x, y)
+
+    def amin(self, a, axis=None):
+        queue = self._array_context.queue
+
+        if isinstance(axis, int):
+            axis = axis,
+
+        def _rec_min(ary):
+            if axis not in [None, tuple(range(ary.ndim))]:
+                raise NotImplementedError(f"Min. over '{axis}' axes not supported.")
+            return cl_array.min(ary, queue=queue)
+
+        result = rec_map_reduce_array_container(
+                partial(reduce, partial(cl_array.minimum, queue=queue)),
+                _rec_min,
+                a)
+
+        if not self._array_context._force_device_scalars:
+            result = result.get()[()]
+        return result
+
+    min = amin
+
+    def absolute(self, a):
+        return self.abs(a)
+
+    # }}}
+
+    # {{{ sorting, searching, and counting
+
+    def where(self, criterion, then, else_):
+        def where_inner(inner_crit, inner_then, inner_else):
+            if isinstance(inner_crit, bool):
+                return inner_then if inner_crit else inner_else
+            return cl_array.if_positive(inner_crit != 0, inner_then, inner_else,
+                    queue=self._array_context.queue)
+
+        return rec_multimap_array_container(where_inner, criterion, then, else_)
+
+    # }}}
+
+
 # }}}
 
 
diff --git a/arraycontext/impl/pytato/fake_numpy.py b/arraycontext/impl/pytato/fake_numpy.py
index e1730b93c6834f09f69a19958429ddd8ec40fa4a..babdf4a03231bec238c155cbaed2f34f8d30ba37 100644
--- a/arraycontext/impl/pytato/fake_numpy.py
+++ b/arraycontext/impl/pytato/fake_numpy.py
@@ -54,96 +54,45 @@ class PytatoFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
         :mod:`pytato` does not define any memory layout for its arrays. See
         :ref:`Pytato docs <pytato:memory-layout>` for more on this.
     """
+
+    _pt_funcs = frozenset({
+        "sin", "cos", "tan", "arcsin", "arccos", "arctan",
+        "sinh", "cosh", "tanh", "exp", "log", "log10",
+        "sqrt", "abs", "isnan"
+        })
+
     def _get_fake_numpy_linalg_namespace(self):
         return PytatoFakeNumpyLinalgNamespace(self._array_context)
 
     def __getattr__(self, name):
-
-        pt_funcs = ["abs", "sin", "cos", "tan", "arcsin", "arccos", "arctan",
-                    "sinh", "cosh", "tanh", "exp", "log", "log10", "isnan",
-                    "sqrt", "exp"]
-        if name in pt_funcs:
+        if name in self._pt_funcs:
             from functools import partial
             return partial(rec_map_array_container, getattr(pt, name))
 
         return super().__getattr__(name)
 
-    def reshape(self, a, newshape, order="C"):
-        return rec_map_array_container(
-                lambda ary: pt.reshape(a, newshape, order=order),
-                a)
+    # NOTE: the order of these follows the order in numpy docs
+    # NOTE: when adding a function here, also add it to `array_context.rst` docs!
 
-    def transpose(self, a, axes=None):
-        return rec_multimap_array_container(pt.transpose, a, axes)
-
-    def concatenate(self, arrays, axis=0):
-        return rec_multimap_array_container(pt.concatenate, arrays, axis)
+    # {{{ array creation routines
 
     def ones_like(self, ary):
-        def _ones_like(subary):
-            return pt.ones(subary.shape, subary.dtype)
-
-        return self._new_like(ary, _ones_like)
-
-    def maximum(self, x, y):
-        return rec_multimap_array_container(pt.maximum, x, y)
-
-    def minimum(self, x, y):
-        return rec_multimap_array_container(pt.minimum, x, y)
-
-    def where(self, criterion, then, else_):
-        return rec_multimap_array_container(pt.where, criterion, then, else_)
-
-    def sum(self, a, axis=None, dtype=None):
-        def _pt_sum(ary):
-            if dtype not in [ary.dtype, None]:
-                raise NotImplementedError
-
-            return pt.sum(ary, axis=axis)
-
-        return rec_map_reduce_array_container(sum, _pt_sum, a)
-
-    def min(self, a, axis=None):
-        return rec_map_reduce_array_container(
-                partial(reduce, pt.minimum), partial(pt.amin, axis=axis), a)
-
-    def max(self, a, axis=None):
-        return rec_map_reduce_array_container(
-                partial(reduce, pt.maximum), partial(pt.amax, axis=axis), a)
-
-    def stack(self, arrays, axis=0):
-        return rec_multimap_array_container(
-                lambda *args: pt.stack(arrays=args, axis=axis),
-                *arrays)
-
-    def broadcast_to(self, array, shape):
-        return rec_map_array_container(partial(pt.broadcast_to, shape=shape), array)
-
-    # {{{ relational operators
+        return self.full_like(ary, 1)
 
-    def equal(self, x, y):
-        return rec_multimap_array_container(pt.equal, x, y)
-
-    def not_equal(self, x, y):
-        return rec_multimap_array_container(pt.not_equal, x, y)
-
-    def greater(self, x, y):
-        return rec_multimap_array_container(pt.greater, x, y)
-
-    def greater_equal(self, x, y):
-        return rec_multimap_array_container(pt.greater_equal, x, y)
+    def full_like(self, ary, fill_value):
+        def _full_like(subary):
+            return pt.full(subary.shape, fill_value, subary.dtype)
 
-    def less(self, x, y):
-        return rec_multimap_array_container(pt.less, x, y)
+        return self._new_like(ary, _full_like)
 
-    def less_equal(self, x, y):
-        return rec_multimap_array_container(pt.less_equal, x, y)
+    # }}}
 
-    def conj(self, x):
-        return rec_multimap_array_container(pt.conj, x)
+    # {{{ array manipulation routines
 
-    def arctan2(self, y, x):
-        return rec_multimap_array_container(pt.arctan2, y, x)
+    def reshape(self, a, newshape, order="C"):
+        return rec_map_array_container(
+                lambda ary: pt.reshape(a, newshape, order=order),
+                a)
 
     def ravel(self, a, order="C"):
         """
@@ -167,16 +116,34 @@ class PytatoFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
 
         return rec_map_array_container(_rec_ravel, a)
 
-    def any(self, a):
-        return rec_map_reduce_array_container(
-                partial(reduce, pt.logical_or),
-                lambda subary: pt.any(subary), a)
+    def transpose(self, a, axes=None):
+        return rec_multimap_array_container(pt.transpose, a, axes)
+
+    def broadcast_to(self, array, shape):
+        return rec_map_array_container(partial(pt.broadcast_to, shape=shape), array)
+
+    def concatenate(self, arrays, axis=0):
+        return rec_multimap_array_container(pt.concatenate, arrays, axis)
+
+    def stack(self, arrays, axis=0):
+        return rec_multimap_array_container(
+                lambda *args: pt.stack(arrays=args, axis=axis),
+                *arrays)
+
+    # }}}
+
+    # {{{ logic functions
 
     def all(self, a):
         return rec_map_reduce_array_container(
                 partial(reduce, pt.logical_and),
                 lambda subary: pt.all(subary), a)
 
+    def any(self, a):
+        return rec_map_reduce_array_container(
+                partial(reduce, pt.logical_or),
+                lambda subary: pt.any(subary), a)
+
     def array_equal(self, a, b):
         actx = self._array_context
 
@@ -202,4 +169,69 @@ class PytatoFakeNumpyNamespace(LoopyBasedFakeNumpyNamespace):
 
         return rec_equal(a, b)
 
+    def greater(self, x, y):
+        return rec_multimap_array_container(pt.greater, x, y)
+
+    def greater_equal(self, x, y):
+        return rec_multimap_array_container(pt.greater_equal, x, y)
+
+    def less(self, x, y):
+        return rec_multimap_array_container(pt.less, x, y)
+
+    def less_equal(self, x, y):
+        return rec_multimap_array_container(pt.less_equal, x, y)
+
+    def equal(self, x, y):
+        return rec_multimap_array_container(pt.equal, x, y)
+
+    def not_equal(self, x, y):
+        return rec_multimap_array_container(pt.not_equal, x, y)
+
+    # }}}
+
+    # {{{ mathematical functions
+
+    def arctan2(self, y, x):
+        return rec_multimap_array_container(pt.arctan2, y, x)
+
+    def sum(self, a, axis=None, dtype=None):
+        def _pt_sum(ary):
+            if dtype not in [ary.dtype, None]:
+                raise NotImplementedError
+
+            return pt.sum(ary, axis=axis)
+
+        return rec_map_reduce_array_container(sum, _pt_sum, a)
+
+    def conj(self, x):
+        return rec_multimap_array_container(pt.conj, x)
+
+    def maximum(self, x, y):
+        return rec_multimap_array_container(pt.maximum, x, y)
+
+    def amax(self, a, axis=None):
+        return rec_map_reduce_array_container(
+                partial(reduce, pt.maximum), partial(pt.amax, axis=axis), a)
+
+    max = amax
+
+    def minimum(self, x, y):
+        return rec_multimap_array_container(pt.minimum, x, y)
+
+    def amin(self, a, axis=None):
+        return rec_map_reduce_array_container(
+                partial(reduce, pt.minimum), partial(pt.amin, axis=axis), a)
+
+    min = amin
+
+    def absolute(self, a):
+        return self.abs(a)
+
+    # }}}
+
+    # {{{ sorting, searching, and counting
+
+    def where(self, criterion, then, else_):
+        return rec_multimap_array_container(pt.where, criterion, then, else_)
+
     # }}}
diff --git a/arraycontext/loopy.py b/arraycontext/loopy.py
index 728904e44a17423b8d9ae38883bf73ec23a34550..1f9031834825e1419ecb5a4923f134442692dd43 100644
--- a/arraycontext/loopy.py
+++ b/arraycontext/loopy.py
@@ -148,7 +148,9 @@ class LoopyBasedFakeNumpyNamespace(BaseFakeNumpyNamespace):
                 or name in self._c_to_numpy_arc_functions):
             return multimapped_over_array_containers(loopy_implemented_elwise_func)
         else:
-            raise AttributeError(name)
+            raise AttributeError(
+                    f"'{type(self._array_context).__name__}.np' object "
+                    f"has no attribute '{name}'")
 
 # }}}
 
diff --git a/doc/array_context.rst b/doc/array_context.rst
index 7a12dfb607790494910eb9232e2220e4e9bf8ffa..680d49bf793ade7aba8622e50f7f4f2734e3d289 100644
--- a/doc/array_context.rst
+++ b/doc/array_context.rst
@@ -18,3 +18,16 @@ Lazy/Deferred evaluation array context based on :mod:`pytato`
 -------------------------------------------------------------
 
 .. automodule:: arraycontext.impl.pytato
+
+.. _numpy-coverage:
+
+:mod:`numpy` coverage
+---------------------
+
+This is a list of functionality implemented by :attr:`arraycontext.ArrayContext.np`.
+
+.. note::
+
+   Only functions and methods that have at least one implementation are listed.
+
+.. include:: numpy_coverage.rst
diff --git a/doc/make_numpy_coverage_table.py b/doc/make_numpy_coverage_table.py
new file mode 100644
index 0000000000000000000000000000000000000000..e100f0b3235884fc7eb0a569d4b101e5e75da483
--- /dev/null
+++ b/doc/make_numpy_coverage_table.py
@@ -0,0 +1,247 @@
+"""
+Workflow:
+
+1. If a new array context is implemented, it should be added to
+   :func:`initialize_contexts`.
+
+2. If a new function is implemented, it should be added to the
+   corresponding ``write_section_name`` function.
+
+3. Once everything is added, regenerate the tables using
+
+.. code::
+
+    python make_numpy_support_table.py numpy_coverage.rst
+"""
+
+import pathlib
+from mako.template import Template
+
+import arraycontext
+
+
+# {{{ templating
+
+HEADER = """
+.. raw:: html
+
+    <style> .red {color:red} </style>
+    <style> .green {color:green} </style>
+
+.. role:: red
+.. role:: green
+"""
+
+TABLE_TEMPLATE = Template("""
+${title}
+${'~' * len(title)}
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+% for ctx in contexts:
+      - :class:`~arraycontext.${type(ctx).__name__}`
+% endfor
+% for name, (directive, in_context) in numpy_functions_for_context.items():
+    * - :${directive}:`numpy.${name}`
+    % for ctx in contexts:
+    <%
+    flag = in_context.get(type(ctx), "yes").capitalize()
+    color = "green" if flag == "Yes" else "red"
+    %>  - :${color}:`${flag}`
+    % endfor
+% endfor
+""")
+
+
+def initialize_contexts():
+    import pyopencl as cl
+    ctx = cl.create_some_context()
+    queue = cl.CommandQueue(ctx)
+
+    return [
+            arraycontext.PyOpenCLArrayContext(queue, force_device_scalars=True),
+            arraycontext.PytatoPyOpenCLArrayContext(queue),
+            ]
+
+
+def build_supported_functions(funcs, contexts):
+    import numpy as np
+
+    numpy_functions_for_context = {}
+    for directive, name in funcs:
+        if not hasattr(np, name):
+            raise ValueError(f"'{name}' not found in numpy namespace")
+
+        in_context = {}
+        for ctx in contexts:
+            try:
+                _ = getattr(ctx.np, name)
+            except AttributeError:
+                in_context[type(ctx)] = "No"
+
+        numpy_functions_for_context[name] = (directive, in_context)
+
+    return numpy_functions_for_context
+
+# }}}
+
+
+# {{{ writing
+
+def write_array_creation_routines(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.array-creation.html
+    funcs = (
+            # (sphinx-directive, name)
+            ("func", "empty_like"),
+            ("func", "ones_like"),
+            ("func", "zeros_like"),
+            ("func", "full_like"),
+            ("func", "copy"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Array creation routines",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+
+def write_array_manipulation_routines(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.array-manipulation.html
+    funcs = (
+            # (sphinx-directive, name)
+            ("func", "reshape"),
+            ("func", "ravel"),
+            ("func", "transpose"),
+            ("func", "broadcast_to"),
+            ("func", "concatenate"),
+            ("func", "stack"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Array manipulation routines",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+
+def write_linear_algebra(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.linalg.html
+    funcs = (
+            # (sphinx-directive, name)
+            ("func", "vdot"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Linear algebra",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+
+def write_logic_functions(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.logic.html
+    funcs = (
+            # (sphinx-directive, name)
+            ("func", "all"),
+            ("func", "any"),
+            ("data", "greater"),
+            ("data", "greater_equal"),
+            ("data", "less"),
+            ("data", "less_equal"),
+            ("data", "equal"),
+            ("data", "not_equal"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Logic Functions",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+
+def write_mathematical_functions(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.math.html
+    funcs = (
+            ("data", "sin"),
+            ("data", "cos"),
+            ("data", "tan"),
+            ("data", "arcsin"),
+            ("data", "arccos"),
+            ("data", "arctan"),
+            ("data", "arctan2"),
+            ("data", "sinh"),
+            ("data", "cosh"),
+            ("data", "tanh"),
+            ("data", "floor"),
+            ("data", "ceil"),
+            ("func", "sum"),
+            ("data", "exp"),
+            ("data", "log"),
+            ("data", "log10"),
+            ("func", "real"),
+            ("func", "imag"),
+            ("data", "conjugate"),
+            ("data", "maximum"),
+            ("func", "amax"),
+            ("data", "minimum"),
+            ("func", "amin"),
+            ("data", "sqrt"),
+            ("data", "absolute"),
+            ("data", "fabs"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Mathematical functions",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+
+def write_searching_sorting_and_counting(outf, contexts):
+    # https://numpy.org/doc/stable/reference/routines.sort.html
+    funcs = (
+            ("func", "where"),
+            )
+
+    r = TABLE_TEMPLATE.render(
+            title="Sorting, searching, and counting",
+            contexts=contexts,
+            numpy_functions_for_context=build_supported_functions(funcs, contexts),
+            )
+    outf.write(r)
+
+# }}}
+
+
+if __name__ == "__main__":
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument("filename", nargs="?", type=pathlib.Path, default=None)
+    args = parser.parse_args()
+
+    import sys
+    if args.filename is not None:
+        outf = open(args.filename, "w")
+    else:
+        outf = sys.stdout
+
+    ctxs = initialize_contexts()
+
+    outf.write(HEADER)
+    write_array_creation_routines(outf, ctxs)
+    write_array_manipulation_routines(outf, ctxs)
+    write_linear_algebra(outf, ctxs)
+    write_logic_functions(outf, ctxs)
+    write_mathematical_functions(outf, ctxs)
+
+    if args.filename is not None:
+        outf.close()
diff --git a/doc/numpy_coverage.rst b/doc/numpy_coverage.rst
new file mode 100644
index 0000000000000000000000000000000000000000..8d58fbe79d7d2d054855458bf9652e5d093e0274
--- /dev/null
+++ b/doc/numpy_coverage.rst
@@ -0,0 +1,196 @@
+
+.. raw:: html
+
+    <style> .red {color:red} </style>
+    <style> .green {color:green} </style>
+
+.. role:: red
+.. role:: green
+
+Array creation routines
+~~~~~~~~~~~~~~~~~~~~~~~
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+      - :class:`~arraycontext.PyOpenCLArrayContext`
+      - :class:`~arraycontext.PytatoPyOpenCLArrayContext`
+    * - :func:`numpy.empty_like`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.ones_like`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.zeros_like`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.full_like`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.copy`
+      - :green:`Yes`
+      - :red:`No`
+
+Array manipulation routines
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+      - :class:`~arraycontext.PyOpenCLArrayContext`
+      - :class:`~arraycontext.PytatoPyOpenCLArrayContext`
+    * - :func:`numpy.reshape`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.ravel`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.transpose`
+      - :red:`No`
+      - :green:`Yes`
+    * - :func:`numpy.broadcast_to`
+      - :red:`No`
+      - :green:`Yes`
+    * - :func:`numpy.concatenate`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.stack`
+      - :green:`Yes`
+      - :green:`Yes`
+
+Linear algebra
+~~~~~~~~~~~~~~
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+      - :class:`~arraycontext.PyOpenCLArrayContext`
+      - :class:`~arraycontext.PytatoPyOpenCLArrayContext`
+    * - :func:`numpy.vdot`
+      - :green:`Yes`
+      - :red:`No`
+
+Logic Functions
+~~~~~~~~~~~~~~~
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+      - :class:`~arraycontext.PyOpenCLArrayContext`
+      - :class:`~arraycontext.PytatoPyOpenCLArrayContext`
+    * - :func:`numpy.all`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.any`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.greater`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.greater_equal`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.less`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.less_equal`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.equal`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.not_equal`
+      - :green:`Yes`
+      - :green:`Yes`
+
+Mathematical functions
+~~~~~~~~~~~~~~~~~~~~~~
+
+.. list-table::
+    :header-rows: 1
+
+    * - Function
+      - :class:`~arraycontext.PyOpenCLArrayContext`
+      - :class:`~arraycontext.PytatoPyOpenCLArrayContext`
+    * - :data:`numpy.sin`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.cos`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.tan`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.arcsin`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.arccos`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.arctan`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.arctan2`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.sinh`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.cosh`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.tanh`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.floor`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.ceil`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.sum`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.exp`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.log`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.log10`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.real`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.imag`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.conjugate`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.maximum`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.amax`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.minimum`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :func:`numpy.amin`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.sqrt`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.absolute`
+      - :green:`Yes`
+      - :green:`Yes`
+    * - :data:`numpy.fabs`
+      - :green:`Yes`
+      - :green:`Yes`