diff --git a/docs/source/api/cuml.metrics.rst b/docs/source/api/cuml.metrics.rst
index d9ab5cd434..e182b272ee 100644
--- a/docs/source/api/cuml.metrics.rst
+++ b/docs/source/api/cuml.metrics.rst
@@ -19,6 +19,7 @@ Classification and Distance Metrics
    log_loss
    roc_auc_score
    precision_recall_curve
+   precision_score
    trustworthiness
 
 Regression Metrics
diff --git a/python/cuml/cuml/metrics/__init__.py b/python/cuml/cuml/metrics/__init__.py
index 0048025f50..643c60a033 100644
--- a/python/cuml/cuml/metrics/__init__.py
+++ b/python/cuml/cuml/metrics/__init__.py
@@ -3,7 +3,11 @@
 # SPDX-License-Identifier: Apache-2.0
 #
 
-from cuml.metrics._classification import accuracy_score, log_loss
+from cuml.metrics._classification import (
+    accuracy_score,
+    log_loss,
+    precision_score,
+)
 from cuml.metrics._ranking import precision_recall_curve, roc_auc_score
 from cuml.metrics.cluster.adjusted_rand_index import adjusted_rand_score
 from cuml.metrics.cluster.completeness_score import (
@@ -50,6 +54,7 @@
     "adjusted_rand_score",
     "roc_auc_score",
     "precision_recall_curve",
+    "precision_score",
     "log_loss",
     "homogeneity_score",
     "completeness_score",
diff --git a/python/cuml/cuml/metrics/_classification.py b/python/cuml/cuml/metrics/_classification.py
index 38e9e1977d..406fd47102 100644
--- a/python/cuml/cuml/metrics/_classification.py
+++ b/python/cuml/cuml/metrics/_classification.py
@@ -2,9 +2,13 @@
 # SPDX-FileCopyrightText: Copyright (c) 2020-2026, NVIDIA CORPORATION.
 # SPDX-License-Identifier: Apache-2.0
 #
+import numbers
+import warnings
+
 import cudf
 import cupy as cp
 import numpy as np
+from sklearn.exceptions import UndefinedMetricWarning
 
 from cuml.internals.validation import (
     check_array,
@@ -46,7 +50,7 @@ def _input_to_cupy_or_cudf_series(x):
                     raise ValueError(
                         f"Expected 1 column but got {out.shape[1]} columns."
                     )
-                out = out.squeeze()  # ensure 1D
+                out = out.reshape(-1)
 
     return out
 
@@ -108,6 +112,325 @@ def accuracy_score(y_true, y_pred, *, sample_weight=None, normalize=True):
         return float(cp.count_nonzero(correct))
 
 
+def precision_score(
+    y_true,
+    y_pred,
+    *,
+    labels=None,
+    pos_label=1,
+    average="binary",
+    sample_weight=None,
+    zero_division="warn",
+):
+    """
+    Compute the precision.
+
+    The precision is the ratio ``tp / (tp + fp)`` where ``tp`` is the number
+    of true positives and ``fp`` the number of false positives. The precision
+    is intuitively the ability of the classifier not to label as positive a
+    sample that is negative.
+
+    The best value is 1 and the worst value is 0.
+
+    Parameters
+    ----------
+    y_true : array-like (device or host) of shape (n_samples,)
+        Ground truth (correct) target values.
+    y_pred : array-like (device or host) of shape (n_samples,)
+        Estimated target values as returned by a classifier.
+    labels : array-like (device or host), default=None
+        The set of labels to include when ``average != 'binary'``, and their
+        order if ``average is None``. Labels present in the data can be
+        excluded, and labels not present in the data will receive the score
+        given by ``zero_division``. Ignored when ``average == 'binary'``.
+    pos_label : int, float, bool or str, default=1
+        The class to report if ``average='binary'`` and the data is binary,
+        otherwise this parameter is ignored.
+    average : {'micro', 'macro', 'weighted', 'binary'} or None, \
+            default='binary'
+        This parameter is required for multiclass targets.
+        ``'micro'``:
+            Calculate metrics globally by counting the total true positives
+            and false positives.
+        ``'macro'``:
+            Calculate metrics for each label, and find their unweighted mean.
+        ``'weighted'``:
+            Calculate metrics for each label, and find their average weighted
+            by support (the number of true instances for each label).
+        ``'binary'``:
+            Only report results for the class specified by ``pos_label``.
+            Only applicable to binary targets.
+        If ``None``, the scores for each label are returned individually.
+    sample_weight : array-like (device or host) of shape (n_samples,), \
+            default=None
+        Sample weights.
+    zero_division : {"warn", 0.0, 1.0}, default="warn"
+        Sets the value to return when there is a zero division. If set to
+        ``"warn"``, this acts like 0, but a warning is also raised.
+
+    Returns
+    -------
+    score : float or numpy.ndarray of float
+        Precision of the positive class in binary classification or the
+        averaged precision of each class for the multiclass task. A NumPy
+        array with one score per label, ordered following ``labels`` (or the
+        sorted union of the observed labels when ``labels is None``), is
+        returned when ``average is None``.
+
+    See Also
+    --------
+    accuracy_score : Accuracy classification score.
+    confusion_matrix : Compute confusion matrix to evaluate the accuracy of a
+        classification.
+
+    Notes
+    -----
+    Numeric labels (integer, whole-number float and bool dtypes) are
+    counted on the GPU. Fractional float targets are rejected, matching
+    scikit-learn's refusal of continuous targets.
+    String, object and categorical labels are supported through a device-side
+    encoding against the sorted union of the observed labels, which matches
+    scikit-learn's label ordering. Null values are not supported. The
+    ``'samples'`` averaging strategy, multilabel indicator input and
+    ``zero_division=np.nan`` accepted by scikit-learn are not supported.
+    NumPy scalar literals such as ``np.int64(0)`` are accepted and
+    substituted as-is, where scikit-learn 1.9 substitutes ``np.nan``.
+    Unlike scikit-learn versions below 1.8, empty targets raise a
+    ``ValueError`` instead of being scored.
+
+    Examples
+    --------
+    .. code-block:: python
+
+        >>> import cupy as cp
+        >>> from cuml.metrics import precision_score
+        >>> y_true = cp.array([0, 1, 2, 0, 1, 2])
+        >>> y_pred = cp.array([0, 2, 1, 0, 0, 1])
+        >>> precision_score(y_true, y_pred, average='macro')
+        0.2222222222222222
+        >>> precision_score(y_true, y_pred, average='micro')
+        0.3333333333333333
+        >>> precision_score(y_true, y_pred, average=None)
+        array([0.66666667, 0.        , 0.        ])
+    """
+
+    average_options = (None, "micro", "macro", "weighted", "binary")
+    if average not in average_options:
+        raise ValueError(f"average has to be one of {average_options}")
+
+    if isinstance(zero_division, str) and zero_division == "warn":
+        zero_division_value = 0.0
+    elif isinstance(zero_division, numbers.Real) and zero_division in (0, 1):
+        zero_division_value = float(zero_division)
+    else:
+        raise ValueError(
+            'zero_division must be one of {"warn", 0, 1}, got '
+            f"{zero_division!r}"
+        )
+
+    y_true = _input_to_cupy_or_cudf_series(y_true)
+    y_pred = _input_to_cupy_or_cudf_series(y_pred)
+
+    check_consistent_length(y_true, y_pred)
+
+    if len(y_true) == 0 or len(y_pred) == 0:
+        raise ValueError(
+            "Found empty input array (e.g., `y_true` or `y_pred`) while a "
+            "minimum of 1 sample is required."
+        )
+
+    for name, y in (("y_true", y_true), ("y_pred", y_pred)):
+        if isinstance(y, cudf.Series) and y.isna().any():
+            raise ValueError(
+                f"precision_score does not support null values in {name}"
+            )
+
+    if (
+        sample_weight := check_sample_weight(sample_weight, dtype=np.float64)
+    ) is not None:
+        check_consistent_length(y_true, sample_weight)
+
+    numeric = all(
+        not isinstance(y, cudf.Series) or y.dtype.kind in "iufb"
+        for y in (y_true, y_pred)
+    )
+
+    if numeric:
+        y_true_t = (
+            y_true.to_cupy() if isinstance(y_true, cudf.Series) else y_true
+        )
+        y_pred_t = (
+            y_pred.to_cupy() if isinstance(y_pred, cudf.Series) else y_pred
+        )
+        present = cp.unique(
+            cp.concatenate([cp.unique(y_true_t), cp.unique(y_pred_t)])
+        )
+        if average == "binary":
+            present_labels = cp.asnumpy(present).tolist()
+        for name, arr in (("y_true", y_true_t), ("y_pred", y_pred_t)):
+            if arr.dtype.kind == "f":
+                if bool(cp.isnan(arr).any()):
+                    raise ValueError(f"Input {name} contains NaN.")
+                if bool(cp.isinf(arr).any()):
+                    raise ValueError(
+                        f"Input {name} contains infinity or a value too "
+                        "large for dtype('float64')."
+                    )
+                if bool((arr != cp.floor(arr)).any()):
+                    raise ValueError(f"'{name}' can only have integer values")
+    else:
+        y_true = (
+            y_true if isinstance(y_true, cudf.Series) else cudf.Series(y_true)
+        )
+        y_pred = (
+            y_pred if isinstance(y_pred, cudf.Series) else cudf.Series(y_pred)
+        )
+        try:
+            present_labels = sorted(
+                set().union(
+                    *[
+                        set(y.unique().dropna().to_pandas().tolist())
+                        for y in (y_true, y_pred)
+                    ]
+                )
+            )
+        except TypeError:
+            raise ValueError(
+                "Mix of label input types (string and number)"
+            ) from None
+
+    if average == "binary":
+        if len(present_labels) > 2:
+            raise ValueError(
+                "Target is multiclass but average='binary'. Please choose "
+                "another average setting, one of [None, 'micro', 'macro', "
+                "'weighted']."
+            )
+        if len(present_labels) >= 2 and pos_label not in present_labels:
+            raise ValueError(
+                f"pos_label={pos_label} is not a valid label. It should be "
+                f"one of {present_labels}"
+            )
+        if pos_label not in present_labels:
+            # an absent positive label can never be predicted
+            if zero_division == "warn":
+                _warn_precision_undefined(1)
+            return zero_division_value
+        out_labels = [pos_label] if not numeric else cp.array([pos_label])
+    else:
+        if pos_label not in (None, 1):
+            warnings.warn(
+                "Note that pos_label (set to "
+                f"{pos_label!r}) is ignored when average != 'binary' "
+                f"(got {average!r}). You may use labels=[pos_label] to "
+                "specify a single positive class.",
+                UserWarning,
+                stacklevel=2,
+            )
+        if labels is not None:
+            out_labels = _labels_as_device_or_host(labels, numeric)
+        else:
+            out_labels = present if numeric else present_labels
+
+    if numeric:
+        table = cp.unique(cp.concatenate([present, out_labels]))
+        pos = cp.searchsorted(table, out_labels)
+        true_idx = cp.searchsorted(table, y_true_t)
+        pred_idx = cp.searchsorted(table, y_pred_t)
+        n_labels_total = table.shape[0]
+    else:
+        table = sorted(set(present_labels) | set(out_labels))
+        cat_dtype = cudf.CategoricalDtype(categories=table)
+        true_idx = (
+            y_true.astype(cat_dtype).cat.codes.to_cupy().astype(np.int64)
+        )
+        pred_idx = (
+            y_pred.astype(cat_dtype).cat.codes.to_cupy().astype(np.int64)
+        )
+        pos = cp.array(
+            [table.index(label) for label in out_labels], dtype=np.int64
+        )
+        n_labels_total = len(table)
+
+    weights = (
+        cp.ones(y_true.shape[0], dtype=cp.float64)
+        if sample_weight is None
+        else sample_weight.astype(cp.float64, copy=False)
+    )
+
+    diag_mask = true_idx == pred_idx
+    tp_sum = cp.bincount(
+        true_idx[diag_mask],
+        weights=weights[diag_mask],
+        minlength=n_labels_total,
+    )[pos]
+    pred_sum = cp.bincount(
+        pred_idx, weights=weights, minlength=n_labels_total
+    )[pos]
+    true_sum = cp.bincount(
+        true_idx, weights=weights, minlength=n_labels_total
+    )[pos]
+
+    empty = pred_sum == 0
+    per_class = cp.where(
+        empty, zero_division_value, tp_sum / cp.where(empty, 1.0, pred_sum)
+    )
+    n_out = per_class.shape[0]
+
+    if average is None:
+        if zero_division == "warn" and empty.any():
+            _warn_precision_undefined(n_out)
+        return cp.asnumpy(per_class)
+
+    if average == "binary":
+        if zero_division == "warn" and pred_sum[0] == 0:
+            _warn_precision_undefined(1)
+        return float(per_class[0])
+
+    if average == "micro":
+        pred_total = pred_sum.sum()
+        if pred_total == 0:
+            if zero_division == "warn":
+                _warn_precision_undefined(1)
+            return zero_division_value
+        return float(tp_sum.sum() / pred_total)
+
+    if zero_division == "warn" and empty.any():
+        _warn_precision_undefined(n_out)
+
+    if average == "macro":
+        return float(per_class.mean())
+
+    if float(true_sum.sum()) == 0.0:
+        # all-zero support: scikit-learn ignores the weights entirely
+        return float(per_class.mean())
+    return float(cp.average(per_class, weights=true_sum))
+
+
+def _labels_as_device_or_host(labels, numeric):
+    if numeric:
+        labels = _input_to_cupy_or_cudf_series(labels)
+        if isinstance(labels, cudf.Series):
+            labels = labels.to_cupy()
+        return cp.reshape(labels, (-1,))
+    if isinstance(labels, cudf.Series):
+        return labels.to_pandas().tolist()
+    if hasattr(labels, "tolist"):
+        return labels.tolist()
+    return list(labels)
+
+
+def _warn_precision_undefined(n_labels):
+    due_to = "due to" if n_labels == 1 else "in labels with"
+    warnings.warn(
+        "Precision is ill-defined and being set to 0.0 "
+        f"{due_to} no predicted samples. Use `zero_division` parameter "
+        "to control this behavior.",
+        UndefinedMetricWarning,
+        stacklevel=3,
+    )
+
+
 def log_loss(
     y_true, y_pred, eps=1e-15, normalize=True, sample_weight=None
 ) -> float:
diff --git a/python/cuml/tests/test_metrics.py b/python/cuml/tests/test_metrics.py
index 59edee6b1c..330948967c 100644
--- a/python/cuml/tests/test_metrics.py
+++ b/python/cuml/tests/test_metrics.py
@@ -23,7 +23,7 @@
 from scipy.stats import entropy as sp_entropy
 from sklearn import preprocessing
 from sklearn.datasets import make_blobs, make_classification
-from sklearn.exceptions import DataConversionWarning
+from sklearn.exceptions import DataConversionWarning, UndefinedMetricWarning
 from sklearn.metrics import confusion_matrix as sk_confusion_matrix
 from sklearn.metrics import hinge_loss as sk_hinge
 from sklearn.metrics import log_loss as sklearn_log_loss
@@ -31,6 +31,7 @@
 from sklearn.metrics import (
     precision_recall_curve as sklearn_precision_recall_curve,
 )
+from sklearn.metrics import precision_score as sk_precision
 from sklearn.metrics import roc_auc_score as sklearn_roc_auc_score
 from sklearn.metrics.cluster import adjusted_rand_score as sk_ars
 from sklearn.metrics.cluster import completeness_score as sk_completeness_score
@@ -57,6 +58,7 @@
     nan_euclidean_distances,
     pairwise_distances,
     precision_recall_curve,
+    precision_score,
     roc_auc_score,
 )
 from cuml.metrics.cluster import adjusted_rand_score as cu_ars
@@ -285,6 +287,346 @@ def test_accuracy_score_scalar_sample_weight():
     ) == cuml.metrics.accuracy_score(y_true, y_pred, normalize=False)
 
 
+@pytest.mark.parametrize(
+    "true_kind, pred_kind",
+    [
+        ("numpy", "numpy"),
+        ("cupy", "cupy"),
+        ("numpy", "cupy"),
+        ("cudf", "cudf"),
+        ("pandas", "pandas"),
+        ("cudf", "numpy"),
+    ],
+)
+@pytest.mark.parametrize(
+    "average", ["binary", "micro", "macro", "weighted", None]
+)
+@pytest.mark.parametrize("n_classes", [2, 5])
+def test_precision_score(true_kind, pred_kind, average, n_classes):
+    N = 60
+    rng = np.random.RandomState(42)
+    np_true = rng.randint(0, n_classes, N)
+    np_pred = rng.randint(0, n_classes, N)
+
+    def convert(x, kind):
+        if kind == "cupy":
+            return cp.array(x)
+        elif kind == "cudf":
+            return cudf.Series(x)
+        elif kind == "pandas":
+            return pd.Series(x)
+        return x
+
+    if average == "binary" and n_classes > 2:
+        with pytest.raises(ValueError, match="Target is multiclass"):
+            precision_score(
+                convert(np_true, true_kind),
+                convert(np_pred, pred_kind),
+                average=average,
+            )
+        return
+
+    res = precision_score(
+        convert(np_true, true_kind),
+        convert(np_pred, pred_kind),
+        average=average,
+    )
+    sol = sk_precision(np_true, np_pred, average=average)
+    if average is None:
+        assert isinstance(res, np.ndarray)
+        np.testing.assert_allclose(res, sol)
+    else:
+        assert isinstance(res, float)
+        assert_almost_equal(res, sol)
+
+
+@pytest.mark.parametrize(
+    "weight_case", [None, "ones", "random", "random_device", "scalar"]
+)
+def test_precision_score_sample_weight(weight_case):
+    N = 40
+    rng = np.random.RandomState(0)
+    np_true = rng.randint(0, 3, N)
+    np_pred = rng.randint(0, 3, N)
+
+    y_true = cp.asarray(np_true)
+    y_pred = cp.asarray(np_pred)
+
+    if weight_case is None:
+        sample_weight = None
+    elif weight_case == "ones":
+        sample_weight = np.ones(N)
+    elif weight_case == "random":
+        sample_weight = rng.rand(N)
+    elif weight_case == "random_device":
+        sample_weight = cp.asarray(rng.rand(N))
+    else:
+        sample_weight = 2.5
+
+    sk_weight = (
+        cp.asnumpy(sample_weight)
+        if isinstance(sample_weight, cp.ndarray)
+        else sample_weight
+    )
+    if sk_weight is not None and np.isscalar(sk_weight):
+        # scikit-learn rejects scalars; a uniform weight is equivalent
+        sk_weight = np.full(N, sk_weight)
+
+    for average in ["micro", "macro", "weighted"]:
+        assert_almost_equal(
+            precision_score(
+                y_true,
+                y_pred,
+                average=average,
+                sample_weight=sample_weight,
+            ),
+            sk_precision(
+                np_true, np_pred, average=average, sample_weight=sk_weight
+            ),
+        )
+
+    np.testing.assert_allclose(
+        precision_score(
+            y_true, y_pred, average=None, sample_weight=sample_weight
+        ),
+        sk_precision(np_true, np_pred, average=None, sample_weight=sk_weight),
+    )
+
+
+def test_precision_score_labels_order():
+    y_true = np.array([0, 1, 2, 2, 0, 1])
+    y_pred = np.array([2, 1, 1, 0, 0, 0])
+
+    res = precision_score(y_true, y_pred, labels=[2, 0], average=None)
+    sol = sk_precision(y_true, y_pred, labels=[2, 0], average=None)
+    assert res.shape == (2,)
+    np.testing.assert_allclose(res, sol)
+
+    sorted_scores = sk_precision(
+        y_true, y_pred, labels=[0, 1, 2], average=None
+    )
+    np.testing.assert_allclose(res, [sorted_scores[2], sorted_scores[0]])
+
+    # a label absent from the data exercises zero_division
+    res = precision_score(
+        y_true, y_pred, labels=[0, 1, 2, 9], average=None, zero_division=0
+    )
+    sol = sk_precision(
+        y_true, y_pred, labels=[0, 1, 2, 9], average=None, zero_division=0
+    )
+    np.testing.assert_allclose(res, sol)
+
+
+def test_precision_score_pos_label():
+    y_true = np.array([0, 2, 2, 0])
+    y_pred = np.array([0, 2, 0, 0])
+
+    assert_almost_equal(
+        precision_score(y_true, y_pred, pos_label=2),
+        sk_precision(y_true, y_pred, pos_label=2),
+    )
+
+    with pytest.raises(ValueError, match="not a valid label"):
+        precision_score(y_true, y_pred, pos_label=7)
+
+
+def test_precision_score_zero_division_warn():
+    y_true = cp.array([0, 1, 1, 0])
+    y_pred = cp.array([0, 0, 0, 0])
+
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        res = precision_score(y_true, y_pred, average=None)
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        sol = sk_precision(
+            cp.asnumpy(y_true), cp.asnumpy(y_pred), average=None
+        )
+    np.testing.assert_allclose(res, sol)
+
+
+@pytest.mark.parametrize(
+    "zero_division", [0, 1, np.int64(0), np.int64(1), np.float32(1.0)]
+)
+def test_precision_score_zero_division_literal(zero_division):
+    # classes 1 and 2 are never predicted
+    y_true = np.array([0, 1, 2, 0, 1, 2])
+    y_pred = np.array([0, 0, 0, 0, 0, 0])
+    literal = float(zero_division)
+
+    # scikit-learn 1.9 routes numpy scalar literals through the nan branch
+    # of _check_zero_division, so the substitution is only compared against
+    # sklearn for exact int and float inputs
+    compare_sklearn = isinstance(zero_division, (int, float))
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        res = precision_score(
+            y_true, y_pred, average=None, zero_division=zero_division
+        )
+        if compare_sklearn:
+            sol = sk_precision(
+                y_true, y_pred, average=None, zero_division=zero_division
+            )
+    if compare_sklearn:
+        np.testing.assert_allclose(res, sol)
+    else:
+        np.testing.assert_allclose(res, [1 / 3, literal, literal])
+
+    # binary case where pos_label is never predicted
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        res = precision_score(
+            np.array([0, 1]), np.array([0, 0]), zero_division=zero_division
+        )
+        if compare_sklearn:
+            sol = sk_precision(
+                np.array([0, 1]),
+                np.array([0, 0]),
+                zero_division=zero_division,
+            )
+    assert res == literal
+    if compare_sklearn:
+        assert sol == literal
+
+
+def test_precision_score_errors():
+    y_true = np.array([0, 1, 2, 0, 1, 2])
+    y_pred = np.array([0, 2, 1, 0, 0, 1])
+
+    with pytest.raises(ValueError, match="Target is multiclass"):
+        precision_score(y_true, y_pred)
+
+    with pytest.raises(ValueError, match="average has to be one of"):
+        precision_score([0, 1], [1, 0], average="invalid")
+
+    with pytest.raises(ValueError, match="not a valid label"):
+        precision_score([0, 1], [1, 0], pos_label=7)
+
+    with pytest.raises(ValueError, match="zero_division must be one of"):
+        precision_score([0, 1], [1, 0], zero_division=0.5)
+
+    with pytest.raises(ValueError, match="can only have integer values"):
+        precision_score(np.array([0.5, 1.5]), np.array([1.5, 0.5]))
+
+    with pytest.raises(ValueError, match="contains NaN"):
+        precision_score(
+            np.array([0.0, np.nan]), np.array([1.0, 0.0]), average="macro"
+        )
+
+    with pytest.raises(ValueError, match="contains infinity"):
+        precision_score(
+            np.array([0.0, np.inf]), np.array([1.0, 0.0]), average="macro"
+        )
+
+    with pytest.raises(ValueError, match="Mix of label input types"):
+        precision_score(np.array([0, 1]), cudf.Series(["a", "b"]))
+
+    with pytest.raises(ValueError, match="empty input array"):
+        precision_score(
+            cp.array([], dtype=cp.int32), cp.array([], dtype=cp.int32)
+        )
+
+    with pytest.raises(ValueError, match="null values"):
+        precision_score(cudf.Series([1, None]), cudf.Series([1, 0]))
+
+    with pytest.raises(ValueError, match="average has to be one of"):
+        precision_score([0, 1], [1, 0], average="samples")
+
+    with pytest.raises(ValueError, match="zero_division must be one of"):
+        precision_score([0, 1], [1, 0], zero_division=np.nan)
+
+
+def test_precision_score_single_class():
+    # scikit-learn scores single-class targets rather than raising
+    y_true = cp.full(4, 3, dtype=cp.int32)
+    y_pred = cp.full(4, 3, dtype=cp.int32)
+
+    assert precision_score(y_true, y_pred, average="macro") == 1.0
+
+    # pos_label absent from a single-class target scores zero_division
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        res = precision_score(y_true, y_pred)
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        sol = sk_precision(cp.asnumpy(y_true), cp.asnumpy(y_pred))
+    assert res == sol == 0.0
+
+
+def test_precision_score_one_column_input():
+    # column vectors and single-sample (1, 1) inputs are flattened to 1D
+    y_true = cp.array([[0], [1], [1], [0]], dtype=cp.int32)
+    y_pred = cp.array([[0], [1], [0], [0]], dtype=cp.int32)
+
+    res = precision_score(y_true, y_pred)
+    sol = sk_precision(cp.asnumpy(y_true), cp.asnumpy(y_pred))
+    assert res == sol
+
+    assert precision_score(cp.array([[1]]), cp.array([[1]])) == 1.0
+
+
+def test_precision_score_single_class_default_pos_label():
+    # scikit-learn skips pos_label validation for one-class targets and
+    # scores the missing pos_label as no predicted samples
+    y_true = cudf.Series(["cat"])
+    y_pred = cudf.Series(["cat"])
+
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        res = precision_score(y_true, y_pred)
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        sol = sk_precision(["cat"], ["cat"])
+    assert res == sol == 0.0
+
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        res = precision_score(
+            cp.array([0], dtype=cp.int32), cp.array([0], dtype=cp.int32)
+        )
+    with pytest.warns(UndefinedMetricWarning, match="ill-defined"):
+        sol = sk_precision(np.array([0]), np.array([0]))
+    assert res == sol == 0.0
+
+
+@pytest.mark.parametrize("to_category", [False, True])
+def test_precision_score_string_labels(to_category):
+    labels = np.array(["a", "b", "c"], dtype="object")
+    rng = np.random.RandomState(42)
+    np_true = labels.take(rng.randint(0, 3, 30))
+    np_pred = labels.take(rng.randint(0, 3, 30))
+
+    y_true = cudf.Series(np_true)
+    y_pred = cudf.Series(np_pred)
+    if to_category:
+        y_true = y_true.astype("category")
+        y_pred = y_pred.astype("category")
+
+    np.testing.assert_allclose(
+        precision_score(y_true, y_pred, average=None),
+        sk_precision(np_true, np_pred, average=None),
+    )
+    assert_almost_equal(
+        precision_score(y_true, y_pred, average="weighted"),
+        sk_precision(np_true, np_pred, average="weighted"),
+    )
+
+    res = precision_score(y_true, y_pred, labels=["c", "a"], average=None)
+    sol = sk_precision(np_true, np_pred, labels=["c", "a"], average=None)
+    np.testing.assert_allclose(res, sol)
+
+
+@pytest.mark.parametrize("n_samples", [unit_param(50), stress_param(500000)])
+def test_precision_score_random(n_samples):
+    upper = 10 if n_samples > 1000 else 3
+    y_true, y_pred, np_true, np_pred = generate_random_labels(
+        lambda rng: rng.randint(0, upper, n_samples), as_cupy=True
+    )
+
+    assert_almost_equal(
+        precision_score(y_true, y_pred, average="macro"),
+        sk_precision(np_true, np_pred, average="macro"),
+    )
+    np.testing.assert_allclose(
+        precision_score(y_true, y_pred, average=None),
+        sk_precision(np_true, np_pred, average=None),
+    )
+
+
 dataset_names = ["noisy_circles", "noisy_moons", "aniso"] + [
     pytest.param(ds, marks=pytest.mark.xfail) for ds in ["blobs", "varied"]
 ]