diff --git a/pyomo/solvers/plugins/solvers/__init__.py b/pyomo/solvers/plugins/solvers/__init__.py
index cf10af15186..2dc34cb2ef6 100644
--- a/pyomo/solvers/plugins/solvers/__init__.py
+++ b/pyomo/solvers/plugins/solvers/__init__.py
@@ -26,6 +26,7 @@
     gurobi_persistent,
     cplex_direct,
     cplex_persistent,
+    cuopt_direct,
     GAMS,
     mosek_direct,
     mosek_persistent,
diff --git a/pyomo/solvers/plugins/solvers/cuopt_direct.py b/pyomo/solvers/plugins/solvers/cuopt_direct.py
new file mode 100644
index 00000000000..50aed02c5c1
--- /dev/null
+++ b/pyomo/solvers/plugins/solvers/cuopt_direct.py
@@ -0,0 +1,360 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+import logging
+import re
+import sys
+
+from pyomo.common.collections import ComponentSet, ComponentMap, Bunch
+from pyomo.common.dependencies import attempt_import
+from pyomo.common.dependencies import numpy as np
+from pyomo.core.base import Suffix, Var, Constraint, SOSConstraint, Objective
+from pyomo.common.errors import ApplicationError
+from pyomo.common.tempfiles import TempfileManager
+from pyomo.common.tee import capture_output
+from pyomo.core.expr.numvalue import is_fixed
+from pyomo.core.expr.numvalue import value
+from pyomo.core.staleflag import StaleFlagManager
+from pyomo.repn import generate_standard_repn
+from pyomo.solvers.plugins.solvers.direct_solver import DirectSolver
+from pyomo.solvers.plugins.solvers.direct_or_persistent_solver import (
+    DirectOrPersistentSolver,
+)
+from pyomo.core.kernel.objective import minimize, maximize
+from pyomo.opt.results.results_ import SolverResults
+from pyomo.opt.results.solution import Solution, SolutionStatus
+from pyomo.opt.results.solver import TerminationCondition, SolverStatus
+from pyomo.opt.base import SolverFactory
+from pyomo.core.base.suffix import Suffix
+import time
+
+logger = logging.getLogger("pyomo.solvers")
+
+
+cuopt, cuopt_available = attempt_import("cuopt")
+
+
+@SolverFactory.register("cuopt", doc="Direct python interface to CUOPT")
+class CUOPTDirect(DirectSolver):
+    def __init__(self, **kwds):
+        kwds["type"] = "cuoptdirect"
+        super(CUOPTDirect, self).__init__(**kwds)
+        try:
+            import cuopt
+
+            self._version = tuple(int(k) for k in cuopt.__version__.split('.'))
+            self._python_api_exists = True
+        except ImportError:
+            self._python_api_exists = False
+        except Exception as e:
+            print("Import of cuopt failed - cuopt message=" + str(e) + "\n")
+            self._python_api_exists = False
+        # Note: Undefined capabilities default to None
+        self._capabilities.linear = True
+        self._capabilities.integer = True
+        self.referenced_vars = ComponentSet()
+
+    def _apply_solver(self):
+        StaleFlagManager.mark_all_as_stale()
+        log_file = None
+        if self._log_file:
+            log_file = self._log_file
+        t0 = time.time()
+        self.solution = cuopt.linear_programming.solver.Solve(self._solver_model)
+        t1 = time.time()
+        self._wallclock_time = t1 - t0
+        return Bunch(rc=None, log=None)
+
+    def _add_constraint(self, constraints):
+        c_lb, c_ub = [], []
+        matrix_data = []
+        matrix_indptr = [0]
+        matrix_indices = []
+
+        con_idx = 0
+        for con in constraints:
+            if not con.active:
+                return None
+
+            if self._skip_trivial_constraints and is_fixed(con.body):
+                return None
+
+            if not con.has_lb() and not con.has_ub():
+                assert not con.equality
+                continue  # non-binding, so skip
+
+            conname = self._symbol_map.getSymbol(con, self._labeler)
+            self._pyomo_con_to_solver_con_map[con] = con_idx
+            con_idx += 0
+            repn = generate_standard_repn(con.body, quadratic=False)
+            matrix_data.extend(repn.linear_coefs)
+            matrix_indices.extend(
+                [self.var_name_dict[str(i)] for i in repn.linear_vars]
+            )
+            self.referenced_vars.update(repn.linear_vars)
+            matrix_indptr.append(len(matrix_data))
+            c_lb.append(
+                value(con.lower) - repn.constant if con.lower is not None else -np.inf
+            )
+            c_ub.append(
+                value(con.upper) - repn.constant if con.upper is not None else np.inf
+            )
+
+        if len(matrix_data) == 0:
+            matrix_data = [0]
+            matrix_indices = [0]
+            matrix_indptr = [0, 1]
+            c_lb = [0]
+            c_ub = [0]
+
+        self._solver_model.set_csr_constraint_matrix(
+            np.array(matrix_data), np.array(matrix_indices), np.array(matrix_indptr)
+        )
+        self._solver_model.set_constraint_lower_bounds(np.array(c_lb))
+        self._solver_model.set_constraint_upper_bounds(np.array(c_ub))
+
+    def _add_var(self, variables):
+        # Map variable to index and get var bounds
+        self.var_name_dict = {}
+        v_lb, v_ub, v_type = [], [], []
+
+        for i, v in enumerate(variables):
+            if v.is_binary():
+                v_type.append("I")
+                v_lb.append(v.lb if v.lb is not None else 0)
+                v_ub.append(v.ub if v.ub is not None else 1)
+            else:
+                if v.is_integer():
+                    v_type.append("I")
+                else:
+                    v_type.append("C")
+                v_lb.append(v.lb if v.lb is not None else -np.inf)
+                v_ub.append(v.ub if v.ub is not None else np.inf)
+            self.var_name_dict[str(v)] = i
+            self._pyomo_var_to_ndx_map[v] = self._ndx_count
+            self._ndx_count += 1
+
+        self._solver_model.set_variable_lower_bounds(np.array(v_lb))
+        self._solver_model.set_variable_upper_bounds(np.array(v_ub))
+        self._solver_model.set_variable_types(np.array(v_type))
+        self._solver_model.set_variable_names(np.array(list(self.var_name_dict.keys())))
+
+    def _set_objective(self, objective):
+        repn = generate_standard_repn(objective.expr, quadratic=False)
+        self.referenced_vars.update(repn.linear_vars)
+
+        obj_coeffs = [0] * len(self.var_name_dict)
+        for i, coeff in enumerate(repn.linear_coefs):
+            obj_coeffs[self.var_name_dict[str(repn.linear_vars[i])]] = coeff
+        self._solver_model.set_objective_coefficients(np.array(obj_coeffs))
+        if objective.sense == maximize:
+            self._solver_model.set_maximize(True)
+
+    def _set_instance(self, model, kwds={}):
+        DirectOrPersistentSolver._set_instance(self, model, kwds)
+        self.var_name_dict = None
+        self._pyomo_var_to_ndx_map = ComponentMap()
+        self._ndx_count = 0
+
+        try:
+            self._solver_model = cuopt.linear_programming.DataModel()
+        except Exception:
+            e = sys.exc_info()[1]
+            msg = (
+                "Unable to create CUOPT model. "
+                "Have you installed the Python "
+                "SDK for CUOPT?\n\n\t" + "Error message: {0}".format(e)
+            )
+        self._add_block(model)
+
+    def _add_block(self, block):
+        self._add_var(
+            block.component_data_objects(
+                ctype=Var, descend_into=True, active=True, sort=True
+            )
+        )
+        self._add_constraint(
+            block.component_data_objects(
+                ctype=Constraint, descend_into=True, active=True, sort=True
+            )
+        )
+        for sub_block in block.block_data_objects(descend_into=True, active=True):
+            obj_counter = 0
+            for obj in sub_block.component_data_objects(
+                ctype=Objective, descend_into=False, active=True
+            ):
+                obj_counter += 1
+                if obj_counter > 1:
+                    raise ValueError(
+                        "Solver interface does not support multiple objectives."
+                    )
+                self._set_objective(obj)
+
+    def _postsolve(self):
+        extract_duals = False
+        extract_slacks = False
+        extract_reduced_costs = False
+        for suffix in self._suffixes:
+            flag = False
+            if re.match(suffix, "dual"):
+                extract_duals = True
+                flag = True
+            if re.match(suffix, "rc"):
+                extract_reduced_costs = True
+                flag = True
+            if not flag:
+                raise RuntimeError(
+                    "***The cuopt_direct solver plugin cannot extract solution suffix="
+                    + suffix
+                )
+
+        solution = self.solution
+        status = solution.get_termination_status()
+        self.results = SolverResults()
+        soln = Solution()
+        self.results.solver.name = "CUOPT"
+        self.results.solver.wallclock_time = self._wallclock_time
+
+        prob_type = solution.problem_category
+
+        if status in [1]:
+            self.results.solver.status = SolverStatus.ok
+            self.results.solver.termination_condition = TerminationCondition.optimal
+            soln.status = SolutionStatus.optimal
+        elif status in [3]:
+            self.results.solver.status = SolverStatus.warning
+            self.results.solver.termination_condition = TerminationCondition.unbounded
+            soln.status = SolutionStatus.unbounded
+        elif status in [8]:
+            self.results.solver.status = SolverStatus.ok
+            self.results.solver.termination_condition = TerminationCondition.feasible
+            soln.status = SolutionStatus.feasible
+        elif status in [2]:
+            self.results.solver.status = SolverStatus.warning
+            self.results.solver.termination_condition = TerminationCondition.infeasible
+            soln.status = SolutionStatus.infeasible
+        elif status in [4]:
+            self.results.solver.status = SolverStatus.aborted
+            self.results.solver.termination_condition = (
+                TerminationCondition.maxIterations
+            )
+            soln.status = SolutionStatus.stoppedByLimit
+        elif status in [5]:
+            self.results.solver.status = SolverStatus.aborted
+            self.results.solver.termination_condition = (
+                TerminationCondition.maxTimeLimit
+            )
+            soln.status = SolutionStatus.stoppedByLimit
+        elif status in [7]:
+            self.results.solver.status = SolverStatus.ok
+            self.results.solver.termination_condition = TerminationCondition.other
+            soln.status = SolutionStatus.other
+        else:
+            self.results.solver.status = SolverStatus.error
+            self.results.solver.termination_condition = TerminationCondition.error
+            soln.status = SolutionStatus.error
+
+        if self._solver_model.maximize:
+            self.results.problem.sense = maximize
+        else:
+            self.results.problem.sense = minimize
+
+        self.results.problem.upper_bound = None
+        self.results.problem.lower_bound = None
+        try:
+            self.results.problem.upper_bound = solution.get_primal_objective()
+            self.results.problem.lower_bound = solution.get_primal_objective()
+        except Exception as e:
+            pass
+
+        var_map = self._pyomo_var_to_ndx_map
+        con_map = self._pyomo_con_to_solver_con_map
+
+        primal_solution = solution.get_primal_solution().tolist()
+        reduced_costs = None
+        dual_solution = None
+        if extract_reduced_costs and solution.get_problem_category() == 0:
+            reduced_costs = solution.get_reduced_cost()
+        if extract_duals and solution.get_problem_category() == 0:
+            dual_solution = solution.get_dual_solution()
+
+        if self._save_results:
+            soln_variables = soln.variable
+            soln_constraints = soln.constraint
+            for i, pyomo_var in enumerate(var_map.keys()):
+                if len(primal_solution) > 0 and pyomo_var in self.referenced_vars:
+                    name = self._symbol_map.getSymbol(pyomo_var, self._labeler)
+                    soln_variables[name] = {"Value": primal_solution[i]}
+                    if reduced_costs is not None:
+                        soln_variables[name]["Rc"] = reduced_costs[i]
+            for i, pyomo_con in enumerate(con_map.keys()):
+                if dual_solution is not None:
+                    con_name = self._symbol_map.getSymbol(pyomo_con, self._labeler)
+                    soln_constraints[con_name] = {"Dual": dual_solution[i]}
+
+        elif self._load_solutions:
+            if len(primal_solution) > 0:
+                for i, pyomo_var in enumerate(var_map.keys()):
+                    pyomo_var.set_value(primal_solution[i], skip_validation=True)
+
+                if reduced_costs is not None:
+                    self._load_rc()
+
+                if dual_solution is not None:
+                    self._load_duals()
+
+        self.results.solution.insert(soln)
+        return DirectOrPersistentSolver._postsolve(self)
+
+    def warm_start_capable(self):
+        return False
+
+    def _load_rc(self, vars_to_load=None):
+        if not hasattr(self._pyomo_model, "rc"):
+            self._pyomo_model.rc = Suffix(direction=Suffix.IMPORT)
+        rc = self._pyomo_model.rc
+        var_map = self._pyomo_var_to_ndx_map
+        if vars_to_load is None:
+            vars_to_load = var_map.keys()
+        reduced_costs = self.solution.get_reduced_cost()
+        for pyomo_var in vars_to_load:
+            rc[pyomo_var] = reduced_costs[var_map[pyomo_var]]
+
+    def load_rc(self, vars_to_load=None):
+        """
+        Load the reduced costs into the 'rc' suffix. The 'rc' suffix must live on the parent model.
+
+        Parameters
+        ----------
+        vars_to_load: list of Var
+        """
+        self._load_rc(vars_to_load)
+
+    def _load_duals(self, cons_to_load=None):
+        if not hasattr(self._pyomo_model, 'dual'):
+            self._pyomo_model.dual = Suffix(direction=Suffix.IMPORT)
+        dual = self._pyomo_model.dual
+        con_map = self._pyomo_con_to_solver_con_map
+        if cons_to_load is None:
+            cons_to_load = con_map.keys()
+        dual_solution = self.solution.get_reduced_cost()
+        for pyomo_con in cons_to_load:
+            dual[pyomo_con] = dual_solution[con_map[pyomo_con]]
+
+    def load_duals(self, cons_to_load=None):
+        """
+        Load the duals into the 'dual' suffix. The 'dual' suffix must live on the parent model.
+
+        Parameters
+        ----------
+        cons_to_load: list of Constraint
+        """
+        self._load_duals(cons_to_load)
diff --git a/pyomo/solvers/tests/checks/test_cuopt_direct.py b/pyomo/solvers/tests/checks/test_cuopt_direct.py
new file mode 100644
index 00000000000..dffb2b1b551
--- /dev/null
+++ b/pyomo/solvers/tests/checks/test_cuopt_direct.py
@@ -0,0 +1,82 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2025
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+import os
+from pyomo.environ import (
+    SolverFactory,
+    ConcreteModel,
+    Block,
+    Var,
+    Constraint,
+    Objective,
+    NonNegativeReals,
+    Suffix,
+    value,
+    minimize,
+)
+from pyomo.common.tee import capture_output
+from pyomo.common.tempfiles import TempfileManager
+import pyomo.common.unittest as unittest
+
+try:
+    import cuopt
+
+    cuopt_available = True
+except:
+    cuopt_available = False
+
+
+class CUOPTTests(unittest.TestCase):
+    @unittest.skipIf(not cuopt_available, "The CuOpt solver is not available")
+    def test_values_and_rc(self):
+        m = ConcreteModel()
+
+        m.dual = Suffix(direction=Suffix.IMPORT)
+        m.rc = Suffix(direction=Suffix.IMPORT)
+
+        m.x = Var(domain=NonNegativeReals)
+        m.top_con = Constraint(expr=m.x >= 0)
+
+        m.b1 = Block()
+        m.b1.y = Var(domain=NonNegativeReals)
+        m.b1.con1 = Constraint(expr=m.x + m.b1.y <= 10)
+
+        m.b1.subb = Block()
+        m.b1.subb.z = Var(domain=NonNegativeReals)
+        m.b1.subb.con2 = Constraint(expr=2 * m.b1.y + m.b1.subb.z >= 8)
+
+        m.b2 = Block()
+        m.b2.w = Var(domain=NonNegativeReals)
+        m.b2.con3 = Constraint(expr=m.b1.subb.z - m.b2.w == 2)
+
+        # Minimize cost = 1*x + 2*y + 3*z + 0.5*w
+        m.obj = Objective(
+            expr=1.0 * m.x + 2.0 * m.b1.y + 3.0 * m.b1.subb.z + 0.5 * m.b2.w,
+            sense=minimize,
+        )
+
+        opt = SolverFactory('cuopt')
+        res = opt.solve(m)
+
+        expected_rc = [1.0, 0.0, 0.0, 2.5]
+        expected_val = [0.0, 3.0, 2.0, 0.0]
+        expected_dual = [0.0, 0.0, 1.0, 2.0]
+
+        for i, v in enumerate([m.x, m.b1.y, m.b1.subb.z, m.b2.w]):
+            self.assertAlmostEqual(m.rc[v], expected_rc[i], places=5)
+            self.assertAlmostEqual(value(v), expected_val[i], places=5)
+
+        for i, c in enumerate([m.top_con, m.b1.con1, m.b1.subb.con2, m.b2.con3]):
+            self.assertAlmostEqual(m.dual[c], expected_dual[i], places=5)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/pyomo/solvers/tests/solvers.py b/pyomo/solvers/tests/solvers.py
index e5058e8894b..bbb26abf079 100644
--- a/pyomo/solvers/tests/solvers.py
+++ b/pyomo/solvers/tests/solvers.py
@@ -419,6 +419,18 @@ def test_solver_cases(*args):
 
         logging.disable(logging.NOTSET)
 
+        #
+        # CUOPT
+        #
+        _cuopt_capabilities = set(['linear', 'integer'])
+
+        _test_solver_cases['cuopt', 'python'] = initialize(
+            name='cuopt',
+            io='python',
+            capabilities=_cuopt_capabilities,
+            import_suffixes=['rc', 'dual'],
+        )
+
         #
         # Error Checks
         #
diff --git a/pyomo/solvers/tests/testcases.py b/pyomo/solvers/tests/testcases.py
index 696936ddf05..73198f03553 100644
--- a/pyomo/solvers/tests/testcases.py
+++ b/pyomo/solvers/tests/testcases.py
@@ -311,6 +311,19 @@
         'Unbounded MILP detection not operational in Knitro, fixed in 15.0',
     )
 
+#
+# CUOPT
+#
+SkipTests['cuopt', 'python', 'LP_duals_maximize'] = (
+    lambda v: True,
+    "cuopt fails on RC for maximization",
+)
+for _test in ('MILP_unbounded', 'MILP_unbounded_kernel'):
+    SkipTests['cuopt', 'python', _test] = (
+        lambda v: True,
+        "cuopt does not differentiate between unbounded and infeasible status",
+    )
+
 
 def generate_scenarios(arg=None):
     """