diff --git a/linopy/io.py b/linopy/io.py
index 2213cbb5..4fd58da8 100644
--- a/linopy/io.py
+++ b/linopy/io.py
@@ -124,6 +124,42 @@ def print_constraint(cons: Any) -> str:
         return print_variable, print_constraint
 
 
+def vectorized_label_names(
+    labels: np.ndarray,
+    prefix: str,
+    printer: Callable | None = None,
+) -> np.ndarray:
+    """
+    Generate label name arrays using vectorized string ops when possible.
+
+    For simple prefix-based names (e.g. "x0", "x1", ..., "c0", "c1", ...),
+    uses np.char operations which are ~100x faster than np.vectorize for large
+    arrays (500K+ elements).
+
+    Falls back to np.vectorize for custom printer functions that require
+    per-element lookups (e.g. explicit coordinate names).
+
+    Parameters
+    ----------
+    labels : np.ndarray
+        Integer label array.
+    prefix : str
+        Single-character prefix ("x" for variables, "c" for constraints).
+        Only used when printer is None.
+    printer : callable, optional
+        Custom scalar printer function. If provided, falls back to
+        np.vectorize (needed for explicit_coordinate_names mode).
+
+    Returns
+    -------
+    np.ndarray
+        Object array of string names.
+    """
+    if printer is not None:
+        return np.vectorize(printer)(labels).astype(object)
+    return np.char.add(prefix, labels.astype(str)).astype(object)
+
+
 def get_printers(
     m: Model, explicit_coordinate_names: bool = False
 ) -> tuple[Callable, Callable]:
@@ -665,7 +701,9 @@ def to_mosek(
     # for j, n in enumerate(("x" + M.vlabels.astype(str).astype(object))):
     #    task.putvarname(j, n)
 
-    labels = np.vectorize(print_variable)(M.vlabels).astype(object)
+    var_printer = print_variable if explicit_coordinate_names else None
+    con_printer = print_constraint if explicit_coordinate_names else None
+    labels = vectorized_label_names(M.vlabels, "x", var_printer)
     task.generatevarnames(
         np.arange(0, len(labels)), "%0", [len(labels)], None, [0], list(labels)
     )
@@ -704,7 +742,7 @@ def to_mosek(
     ## Constraints
 
     if len(m.constraints) > 0:
-        names = np.vectorize(print_constraint)(M.clabels).astype(object)
+        names = vectorized_label_names(M.clabels, "c", con_printer)
         for i, n in enumerate(names):
             task.putconname(i, n)
         bkc = [
@@ -773,7 +811,9 @@ def to_gurobipy(
 
     M = m.matrices
 
-    names = np.vectorize(print_variable)(M.vlabels).astype(object)
+    var_printer = print_variable if explicit_coordinate_names else None
+    con_printer = print_constraint if explicit_coordinate_names else None
+    names = vectorized_label_names(M.vlabels, "x", var_printer)
     kwargs = {}
     if (
         len(m.binaries.labels)
@@ -792,7 +832,7 @@ def to_gurobipy(
         model.ModelSense = -1
 
     if len(m.constraints):
-        names = np.vectorize(print_constraint)(M.clabels).astype(object)
+        names = vectorized_label_names(M.clabels, "c", con_printer)
         c = model.addMConstr(M.A, x, M.sense, M.b)  # type: ignore
         c.setAttr("ConstrName", list(names))  # type: ignore
 
@@ -881,9 +921,11 @@ def to_highspy(m: Model, explicit_coordinate_names: bool = False) -> Highs:
         h.addRows(num_cons, lower, upper, A.nnz, A.indptr, A.indices, A.data)
 
     lp = h.getLp()
-    lp.col_names_ = np.vectorize(print_variable)(M.vlabels).astype(object)
+    var_printer = print_variable if explicit_coordinate_names else None
+    con_printer = print_constraint if explicit_coordinate_names else None
+    lp.col_names_ = vectorized_label_names(M.vlabels, "x", var_printer)
     if len(M.clabels):
-        lp.row_names_ = np.vectorize(print_constraint)(M.clabels).astype(object)
+        lp.row_names_ = vectorized_label_names(M.clabels, "c", con_printer)
     h.passModel(lp)
 
     # quadrative objective