diff --git a/src/pgraph/PGraph.py b/src/pgraph/PGraph.py
index ae35fd1e..8ad82069 100644
--- a/src/pgraph/PGraph.py
+++ b/src/pgraph/PGraph.py
@@ -219,9 +219,8 @@ def Adjacency(
         :type names: list(N) of str, optional
         :param kwargs: optional arguments passed to the graph constructor, for
             example ``vertex_class`` or ``edge_class``
-        :raises ValueError: ``A`` is not square, the length of ``names`` or
-            the number of rows of ``coords`` does not match ``A``, or ``A``
-            has a non-zero diagonal element and the graph is a ``DGraph``
+        :raises ValueError: ``A`` is not square, or the length of ``names`` or
+            the number of rows of ``coords`` does not match ``A``
         :raises TypeError: ``vertex_class`` or ``edge_class`` is not a
             suitable class, see :meth:`_BaseGraph.__init__`
         :return: graph
@@ -230,12 +229,14 @@ def Adjacency(
         Create a directed or undirected graph where non-zero elements ``A[i,j]``
         correspond to edges from vertex ``i`` to vertex ``j``.
 
+        A non-zero diagonal element ``A[i,i]`` is a self-loop on vertex ``i``,
+        for both kinds of graph.
+
         .. warning:: For undirected graph ``A`` should be symmetric but this
-            is not checked.  Only the strictly upper triangular part is used,
-            so the diagonal is ignored and no self-loops are created. A
-            directed graph rejects a non-zero diagonal with ``ValueError``.
-            Edge counts greater than one are not turned into parallel edges;
-            a non-zero element is a single edge with that element as its cost.
+            is not checked.  Only the upper triangular part, including the
+            diagonal, is used. Elements are costs, not edge counts: a non-zero
+            element is a single edge with that element as its cost, so
+            parallel edges cannot be created this way.
 
         .. runblock:: pycon
 
@@ -267,18 +268,16 @@ def Adjacency(
             g.add_vertex(name=name, coord=coord)
 
         if isinstance(g, UGraph):
-            # undirected graph
+            # undirected graph: upper triangle including the diagonal
             for i in range(A.shape[0]):
-                for j in range(i + 1, A.shape[1]):
+                for j in range(i, A.shape[1]):
                     if A[i, j] > 0:
                         g[i].connect(g[j], cost=A[i, j])
         else:
-            # directed graph
+            # directed graph: whole matrix, a diagonal element is a self-loop
             for i in range(A.shape[0]):
                 for j in range(A.shape[1]):
                     if A[i, j] > 0:
-                        if i == j:
-                            raise ValueError("loops in graph not supported")
                         g[i].connect(g[j], cost=A[i, j])
 
         return g
diff --git a/tests/test_adjacency_loops.py b/tests/test_adjacency_loops.py
new file mode 100644
index 00000000..3587accb
--- /dev/null
+++ b/tests/test_adjacency_loops.py
@@ -0,0 +1,89 @@
+"""
+A non-zero diagonal element of the matrix given to ``Adjacency()`` is a
+self-loop, for undirected and directed graphs alike.  The package permits
+self-loops everywhere else (``connect``, ``add_edge``), so the factory must
+not be the one place that silently drops, or refuses, them.
+"""
+
+import unittest
+
+import numpy as np
+
+from pgraph import DGraph, UGraph
+
+
+class _AdjacencyLoops:
+    Graph: type
+
+    def setUp(self):
+        #           0  1  2
+        self.A = np.array(
+            [
+                [0, 1, 0],  # 0 -- 1
+                [1, 4, 1],  # 1 has a self-loop of cost 4, and -- 2
+                [0, 1, 0],
+            ]
+        )
+
+    def loops(self, g):
+        return [e for e in g.edges() if e.v1 is e.v2]
+
+    def test_diagonal_makes_a_self_loop(self):
+        g = self.Graph.Adjacency(self.A)
+
+        self.assertEqual(g.n, 3)
+        self.assertEqual(g.ne, self.ne_with_loop)
+        loops = self.loops(g)
+        self.assertEqual(len(loops), 1)
+        self.assertIs(loops[0].v1, g[1])
+        self.assertEqual(loops[0].cost, 4)
+
+    def test_no_diagonal_no_loop(self):
+        A = self.A.copy()
+        A[1, 1] = 0
+        g = self.Graph.Adjacency(A)
+        self.assertEqual(self.loops(g), [])
+        self.assertEqual(g.ne, self.ne_without_loop)
+
+    def test_loop_is_a_real_edge(self):
+        g = self.Graph.Adjacency(self.A)
+        self.assertTrue(g.iscyclic())
+        self.assertEqual(g.distance()[1, 1], 4)
+
+    def test_other_validation_unchanged(self):
+        with self.assertRaises(ValueError):
+            self.Graph.Adjacency(np.zeros((2, 3)))
+        with self.assertRaises(ValueError):
+            self.Graph.Adjacency(np.zeros((2, 2)), names=["a"])
+        with self.assertRaises(ValueError):
+            self.Graph.Adjacency(np.zeros((2, 2)), coords=np.zeros((3, 2)))
+
+
+class TestUGraphAdjacencyLoops(_AdjacencyLoops, unittest.TestCase):
+    Graph = UGraph
+    # the symmetric matrix is read once, upper triangle: 0--1, 1--2 and 1--1
+    ne_with_loop = 3
+    ne_without_loop = 2
+
+    def test_lower_triangle_still_ignored(self):
+        # documented: only the upper triangle (and now the diagonal) is read
+        A = np.array([[0, 0], [5, 0]])
+        g = UGraph.Adjacency(A)
+        self.assertEqual(g.ne, 0)
+
+
+class TestDGraphAdjacencyLoops(_AdjacencyLoops, unittest.TestCase):
+    Graph = DGraph
+    # the whole matrix is read: 0->1, 1->0, 1->2, 2->1 and 1->1
+    ne_with_loop = 5
+    ne_without_loop = 4
+
+    def test_directed_uses_whole_matrix(self):
+        A = np.array([[0, 0], [5, 0]])
+        g = DGraph.Adjacency(A)
+        self.assertEqual(g.ne, 1)
+        self.assertEqual(g[1].neighbours(), [g[0]])
+
+
+if __name__ == "__main__":
+    unittest.main()