scipy
diff --git a/‎scipy/linalg/basic.py
+36-23 b/‎scipy/linalg/basic.py
+36-23
diff --git a/‎scipy/linalg/tests/test_basic.py
+22 b/‎scipy/linalg/tests/test_basic.py
+22
@@ -1256,11 +1256,16 @@ def pinv(a, cond=None, rcond=None, return_rank=False, check_finite=True):
     a : (M, N) array_like
         Matrix to be pseudo-inverted.
     cond, rcond : float, optional
-        Cutoff for 'small' singular values in the least-squares solver.
-        Singular values smaller than ``rcond * largest_singular_value``
-        are considered zero.
-        If None, it is set to ``np.finfo(a.dtype).eps``.
-        If `a` is an array of integers, it is set to ``np.finfo('float64').eps``.
+        Cutoff factor for 'small' singular values. In `lstsq`, 
+        singular values less than ``cond*largest_singular_value`` will be
+        considered as zero. If both are omitted, the default value
+        ``max(M, N) * eps`` is passed to `lstsq` where ``eps`` is the
+        corresponding machine precision value of the datatype of ``a``.
+
+        .. versionchanged:: 1.3.0
+            Previously the default cutoff value was just `eps` without the
+            factor ``max(M, N)``.
+
     return_rank : bool, optional
         if True, return the effective rank of the matrix
     check_finite : bool, optional
@@ -1293,9 +1298,13 @@ def pinv(a, cond=None, rcond=None, return_rank=False, check_finite=True):
     """
     a = _asarray_validated(a, check_finite=check_finite)
     b = np.identity(a.shape[0], dtype=a.dtype)
+
     if rcond is not None:
         cond = rcond
 
+    if cond is None:
+        cond = max(a.shape) * np.spacing(a.real.dtype.type(1))
+
     x, resids, rank, s = lstsq(a, b, cond=cond, check_finite=False)
 
     if return_rank:
@@ -1317,12 +1326,15 @@ def pinv2(a, cond=None, rcond=None, return_rank=False, check_finite=True):
     a : (M, N) array_like
         Matrix to be pseudo-inverted.
     cond, rcond : float or None
-        Cutoff for 'small' singular values.
-        Singular values smaller than ``rcond*largest_singular_value``
-        are considered zero.
-        If None and the dtype of `a` is ``np.float32``, it is set to
-        ``np.finfo('float32').eps * 1e3``.
-        Otherwise, it is set to ``np.finfo('float64').eps * 1e6``.
+        Cutoff for 'small' singular values; singular values smaller than this
+        value are considered as zero. If both are omitted, the default value
+        ``max(M,N)*largest_singular_value*eps`` is used where ``eps`` is the
+        machine precision value of the datatype of ``a``.
+
+        .. versionchanged:: 1.3.0
+            Previously the default cutoff value was just ``eps*f`` where ``f``
+            was ``1e3`` for single precision and ``1e6`` for double precision.
+
     return_rank : bool, optional
         If True, return the effective rank of the matrix.
     check_finite : bool, optional
@@ -1360,10 +1372,9 @@ def pinv2(a, cond=None, rcond=None, return_rank=False, check_finite=True):
         cond = rcond
     if cond in [None, -1]:
         t = u.dtype.char.lower()
-        factor = {'f': 1E3, 'd': 1E6}
-        cond = factor[t] * np.finfo(t).eps
+        cond = np.max(s) * max(a.shape) * np.finfo(t).eps
 
-    rank = np.sum(s > cond * np.max(s))
+    rank = np.sum(s > cond)
 
     u = u[:, :rank]
     u /= s[:rank]
@@ -1389,12 +1400,15 @@ def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False,
     a : (N, N) array_like
         Real symmetric or complex hermetian matrix to be pseudo-inverted
     cond, rcond : float or None
-        Cutoff for 'small' singular values.
-        Singular values smaller than ``rcond*largest_singular_value``
-        are considered zero.
-        If None and the dtype of `a` is ``np.float32``, it is set to
-        ``np.finfo('float32').eps * 1e3``.
-        Otherwise, it is set to ``np.finfo('float64').eps * 1e6``.
+        Cutoff for 'small' singular values; singular values smaller than this
+        value are considered as zero. If both are omitted, the default
+        ``max(M,N)*largest_eigenvalue*eps`` is used where ``eps`` is the
+        machine precision value of the datatype of 
6D40
``a``.
+
+        .. versionchanged:: 1.3.0
+            Previously the default cutoff value was just ``eps*f`` where ``f``
+            was ``1e3`` for single precision and ``1e6`` for double precision.
+
     lower : bool, optional
         Whether the pertinent array data is taken from the lower or upper
         triangle of `a`. (Default: lower)
@@ -1436,11 +1450,10 @@ def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False,
         cond = rcond
     if cond in [None, -1]:
         t = u.dtype.char.lower()
-        factor = {'f': 1E3, 'd': 1E6}
-        cond = factor[t] * np.finfo(t).eps
+        cond = np.max(np.abs(s)) * max(a.shape) * np.finfo(t).eps
 
     # For Hermitian matrices, singular values equal abs(eigenvalues)
-    above_cutoff = (abs(s) > cond * np.max(abs(s)))
+    above_cutoff = (abs(s) > cond)
     psigma_diag = 1.0 / s[above_cutoff]
     u = u[:, above_cutoff]
 
 
@@ -1279,6 +1279,28 @@ def test_native_list_argument(self):
         assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3))
 
 
+def test_pinv_pinv2_comparison():  # As reported in gh-8861
+    I_6 = np.eye(6)
+    Ts = np.diag([-1] * 4 + [-2], k=-1) + np.diag([-2] + [-1] * 4, k=1)
+    T = I_6 + Ts
+    A = 25 * (np.kron(I_6, T) + np.kron(Ts, I_6))
+
+    Ap, Ap2 = pinv(A), pinv2(A)
+
+    tol = 1e-11
+    assert_allclose(A @ Ap @ A - A, A @ Ap2 @ A - A, rtol=0., atol=tol)
+    assert_allclose(Ap @ A @ Ap - Ap, Ap2 @ A @ Ap2 - Ap2, rtol=0., atol=tol)
+
+
+@pytest.mark.parametrize('scale', (1e-20, 1., 1e20))
+@pytest.mark.parametrize('pinv_', (pinv, pinvh, pinv2))
+def test_auto_rcond(scale, pinv_):
+    x = np.array([[1, 0], [0, 1e-10]]) * scale
+    expected = np.diag(1. / np.diag(x))
+    x_inv = pinv_(x)
+    assert_allclose(x_inv, expected)
+
+
 class TestVectorNorms(object):
 
     def test_types(self):