From 441b7945aaf1ac4bee274e7a2adc69ab29c93aeb Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:53:07 -0400
Subject: [PATCH 01/38] Create .gitignore

---
 .gitignore | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 .gitignore

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..845cda6
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,55 @@
+# Prerequisites
+*.d
+
+# Object files
+*.o
+*.ko
+*.obj
+*.elf
+
+# Linker output
+*.ilk
+*.map
+*.exp
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Libraries
+*.lib
+*.a
+*.la
+*.lo
+
+# Shared objects (inc. Windows DLLs)
+*.dll
+*.so
+*.so.*
+*.dylib
+
+# Executables
+*.exe
+*.out
+*.app
+*.i*86
+*.x86_64
+*.hex
+
+# Debug files
+*.dSYM/
+*.su
+*.idb
+*.pdb
+
+# Kernel Module Compile Results
+*.mod*
+*.cmd
+.tmp_versions/
+modules.order
+Module.symvers
+Mkfile.old
+dkms.conf
+
+# debug information files
+*.dwo

From d4df5c7ac65b0dc2c4e085ecf30795bc84a53264 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:55:06 -0400
Subject: [PATCH 02/38] Update .gitignore

---
 .gitignore | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/.gitignore b/.gitignore
index 845cda6..52defc1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -53,3 +53,7 @@ dkms.conf
 
 # debug information files
 *.dwo
+
+# ignored files
+build/*
+test/*

From e2627ed5f4b7b8ee7fb8ee739d0822a4a085186c Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Tue, 16 Sep 2025 23:04:36 -0400
Subject: [PATCH 03/38] Add .gitignore

---
 .gitignore | 60 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 60 insertions(+)
 create mode 100644 .gitignore

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..9c13904
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,60 @@
+
+# Prerequisites
+*.d
+
+# Object files
+*.o
+*.ko
+*.obj
+*.elf
+
+# Linker output
+*.ilk
+*.map
+*.exp
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Libraries
+*.lib
+*.a
+*.la
+*.lo
+
+# Shared objects (inc. Windows DLLs)
+*.dll
+*.so
+*.so.*
+*.dylib
+
+# Executables
+*.exe
+*.out
+*.app
+*.i*86
+*.x86_64
+*.hex
+
+# Debug files
+*.dSYM/
+*.su
+*.idb
+*.pdb
+
+# Kernel Module Compile Results
+*.mod*
+*.cmd
+.tmp_versions/
+modules.order
+Module.symvers
+Mkfile.old
+dkms.conf
+
+# debug information files
+*.dwo
+
+# ignored files
+build/*
+test/*

From a086822067b726f2e1fe681744d6435764c95fdd Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 02:03:02 -0400
Subject: [PATCH 04/38] Init interface

---
 .gitignore          |  1 +
 Makefile            | 73 +++++++++++++++++++++++++++++++++++++++++++--
 cupdlpx/interface.c | 18 +++++++++++
 cupdlpx/interface.h | 20 +++++++++++++
 4 files changed, 110 insertions(+), 2 deletions(-)
 create mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.h

diff --git a/.gitignore b/.gitignore
index 52defc1..654154c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,3 +57,4 @@ dkms.conf
 # ignored files
 build/*
 test/*
+/.vscode
diff --git a/Makefile b/Makefile
index 3fccf67..74c0c5b 100644
--- a/Makefile
+++ b/Makefile
@@ -16,6 +16,7 @@ NVCCFLAGS = -I. -I$(CUDA_HOME)/include -O3 -g -gencode arch=compute_90,code=sm_9
 # LDFLAGS for the linker
 LDFLAGS = -L$(CONDA_PREFIX)/lib -L$(CUDA_HOME)/lib64 -lcudart -lcusparse -lcublas -lz -lm
 
+# Source discovery (exclude the debug main)
 C_SOURCES = $(filter-out $(SRC_DIR)/cupdlpx.c, $(wildcard $(SRC_DIR)/*.c))
 CU_SOURCES = $(wildcard $(SRC_DIR)/*.cu)
 
@@ -25,22 +26,41 @@ OBJECTS = $(C_OBJECTS) $(CU_OBJECTS)
 
 TARGET_LIB = $(BUILD_DIR)/libcupdlpx.a
 
+# Debug executable (optional)
 DEBUG_SRC = $(SRC_DIR)/cupdlpx.c
 DEBUG_EXEC = $(BUILD_DIR)/cupdlpx
 
-.PHONY: all clean build
+# Tests auto-discovery
+TEST_DIR := ./test
+TEST_BUILD_DIR := $(BUILD_DIR)/tests
 
+TEST_CU_SOURCES := $(wildcard $(TEST_DIR)/*.cu)
+TEST_C_SOURCES := $(wildcard $(TEST_DIR)/*.c)
+
+# Each test source becomes an executable at build/tests/<basename>
+TEST_EXEC_CU := $(patsubst $(TEST_DIR)/%.cu,$(TEST_BUILD_DIR)/%,$(TEST_CU_SOURCES))
+TEST_EXEC_C := $(patsubst $(TEST_DIR)/%.c,$(TEST_BUILD_DIR)/%,$(TEST_C_SOURCES))
+
+# Phony targets
+.PHONY: all clean build tests test run-tests run-test clean-tests
+
+# Default: build the static library
 all: $(TARGET_LIB)
 
+# Archive all objects into the static library
 $(TARGET_LIB): $(OBJECTS)
 	@echo "Archiving objects into $(TARGET_LIB)..."
+	@mkdir -p $(BUILD_DIR)
 	@ar rcs $@ $^
 
+# Build the debug executable (links the library with cupdlpx.c main)
 build: $(DEBUG_EXEC)
+
 $(DEBUG_EXEC): $(DEBUG_SRC) $(TARGET_LIB)
 	@echo "Building debug executable..."
 	@$(LINKER) $(NVCCFLAGS) $(DEBUG_SRC) -o $(DEBUG_EXEC) $(TARGET_LIB) $(LDFLAGS)
 
+# Pattern rules for objects
 $(BUILD_DIR)/%.o: $(SRC_DIR)/%.c
 	@mkdir -p $(BUILD_DIR)
 	@echo "Compiling $< -> $@..."
@@ -51,6 +71,55 @@ $(BUILD_DIR)/%.o: $(SRC_DIR)/%.cu
 	@echo "Compiling $< -> $@..."
 	@$(NVCC) $(NVCCFLAGS) -c $< -o $@
 
+# Build all tests discovered under test/
+test: tests
+tests: $(TARGET_LIB) $(TEST_EXEC_CU) $(TEST_EXEC_C)
+	@echo "All tests built under $(TEST_BUILD_DIR)/"
+
+# Run all tests one by one
+run-tests: tests
+	@echo "Running all tests..."
+	@set -e; \
+	for t in $(TEST_EXEC_CU) $(TEST_EXEC_C); do \
+	  echo "=== $$t ==="; \
+	  "$$t" || exit $$?; \
+	  echo; \
+	done
+
+# Run a single test by basename: make run-test name=<basename>
+run-test: tests
+	@if [ -z "$(name)" ]; then \
+	  echo "Usage: make run-test name=<basename-of-test-file>"; exit 2; \
+	fi
+	@if [ -x "$(TEST_BUILD_DIR)/$(name)" ]; then \
+	  echo "=== $(TEST_BUILD_DIR)/$(name) ==="; \
+	  "$(TEST_BUILD_DIR)/$(name)"; \
+	else \
+	  echo "Executable not found: $(TEST_BUILD_DIR)/$(name)"; \
+	  echo "Did you 'make tests' and is there a test source named '$(TEST_DIR)/$(name).c(u)'?"; \
+	  exit 1; \
+	fi
+
+# Build rule for CUDA tests: compile and link directly with nvcc
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.cu $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building CUDA test $< -> $@..."
+	@$(LINKER) $(NVCCFLAGS) -I$(SRC_DIR) $< -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+# Build rule for C tests: compile with gcc, link with nvcc to get CUDA libs
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.c $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building C test $< -> $@..."
+	@$(CC) $(CFLAGS) -I$(SRC_DIR) -c $< -o $(TEST_BUILD_DIR)/$*.o
+	@$(LINKER) $(NVCCFLAGS) $(TEST_BUILD_DIR)/$*.o -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+
+# Cleaning
+clean-tests:
+	@echo "Cleaning test executables..."
+	@rm -rf $(TEST_BUILD_DIR)
+
 clean:
 	@echo "Cleaning up..."
-	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
\ No newline at end of file
+	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
+	@rm -rf $(TEST_BUILD_DIR)
\ No newline at end of file
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
new file mode 100644
index 0000000..6313091
--- /dev/null
+++ b/cupdlpx/interface.c
@@ -0,0 +1,18 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
new file mode 100644
index 0000000..8d408f5
--- /dev/null
+++ b/cupdlpx/interface.h
@@ -0,0 +1,20 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#pragma once
+
+// Temporary placeholder for interface API
+// Functions will be declared here later
\ No newline at end of file

From 148ca7729768eb27aa9c8c399f6c3c1117a8c44f Mon Sep 17 00:00:00 2001
From: Zedong <zdpeng95@gmail.com>
Date: Wed, 17 Sep 2025 15:46:59 -0400
Subject: [PATCH 05/38] Update command path in README for solver usage

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index be630f0..f0cb248 100644
--- a/README.md
+++ b/README.md
@@ -46,7 +46,7 @@ If the solver runs and creates output files, your installation is successful.
 The solver is invoked with the following syntax, specifying an input file and an output directory.
 
 ```bash
-./cupdlpx [OPTIONS] <mps_file> <output_directory>
+./build/cupdlpx [OPTIONS] <mps_file> <output_directory>
 ```
 
 ### Arguments

From fb1e98c26fe27bf09be132f11aa78d870a5b3559 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 16:14:22 -0400
Subject: [PATCH 06/38] New test: LP model

---
 test/test_interface.c | 63 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100644 test/test_interface.c

diff --git a/test/test_interface.c b/test/test_interface.c
new file mode 100644
index 0000000..e7d0543
--- /dev/null
+++ b/test/test_interface.c
@@ -0,0 +1,63 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "cupdlpx/interface.h"
+#include <stdio.h>
+#include <math.h> 
+
+int main() {
+    // Example: min c^T x
+    // s.t. l <= A x <= u, x >= 0
+
+    int m = 3; // number of constraints
+    int n = 2; // number of variables
+
+    // A as a dense matrix
+    double A[3][2] = {
+        {1.0, 2.0},
+        {0.0, 1.0},
+        {3.0, 2.0}
+    };
+
+    // c: objective coefficients
+    double c[2] = {1.0, 1.0};
+
+    // l&u: constraintbounds
+    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+
+    printf("Objective c = [");
+    for (int j = 0; j < n; j++) {
+        printf(" %g", c[j]);
+    }
+    printf(" ]\n");
+
+    printf("Matrix A:\n");
+    for (int i = 0; i < m; i++) {
+        printf("  [");
+        for (int j = 0; j < n; j++) {
+            printf(" %g", A[i][j]);
+        }
+        printf(" ]\n");
+    }
+
+    printf("Constraint bounds (l <= A x <= b):\n");
+    for (int i = 0; i < m; i++) {
+        printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
+    }
+
+    return 0;
+}

From 610587aa9bb939d08630c196585e5bd71daa579d Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 05:59:37 -0400
Subject: [PATCH 07/38] New feat: interface

---
 cupdlpx/cupdlpx.c     |  25 ------
 cupdlpx/interface.c   |  18 -----
 cupdlpx/interface.cu  | 184 ++++++++++++++++++++++++++++++++++++++++++
 cupdlpx/interface.h   |  34 +++++++-
 cupdlpx/solver.cu     |  25 ++++++
 cupdlpx/solver.h      |   2 +
 cupdlpx/struct.h      |  56 +++++++++++++
 test/test_interface.c |  68 +++++++++++++++-
 8 files changed, 365 insertions(+), 47 deletions(-)
 delete mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.cu

diff --git a/cupdlpx/cupdlpx.c b/cupdlpx/cupdlpx.c
index 25572f7..8900404 100644
--- a/cupdlpx/cupdlpx.c
+++ b/cupdlpx/cupdlpx.c
@@ -134,31 +134,6 @@ void save_solver_summary(const pdhg_solver_state_t *solver_state, const char *ou
     free(file_path);
 }
 
-void set_default_parameters(pdhg_parameters_t *params)
-{
-    params->l_inf_ruiz_iterations = 10;
-    params->has_pock_chambolle_alpha = true;
-    params->pock_chambolle_alpha = 1.0;
-    params->bound_objective_rescaling = true;
-    params->verbose = false;
-    params->termination_evaluation_frequency = 200;
-    params->reflection_coefficient = 1.0;
-
-    params->termination_criteria.eps_optimal_relative = 1e-4;
-    params->termination_criteria.eps_feasible_relative = 1e-4;
-    params->termination_criteria.eps_infeasible = 1e-10;
-    params->termination_criteria.time_sec_limit = 3600.0;
-    params->termination_criteria.iteration_limit = INT32_MAX;
-
-    params->restart_params.artificial_restart_threshold = 0.36;
-    params->restart_params.sufficient_reduction_for_restart = 0.2;
-    params->restart_params.necessary_reduction_for_restart = 0.5;
-    params->restart_params.k_p = 0.99;
-    params->restart_params.k_i = 0.01;
-    params->restart_params.k_d = 0.0;
-    params->restart_params.i_smooth = 0.3;
-}
-
 void print_usage(const char *prog_name)
 {
     fprintf(stderr, "Usage: %s [OPTIONS] <mps_file> <output_dir>\n\n", prog_name);
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
deleted file mode 100644
index 6313091..0000000
--- a/cupdlpx/interface.c
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
-Copyright 2025 Haihao Lu
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "interface.h"
-
diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
new file mode 100644
index 0000000..d361169
--- /dev/null
+++ b/cupdlpx/interface.cu
@@ -0,0 +1,184 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+#include "solver.h"
+#include "utils.h"
+#include <cuda_runtime.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+// helper function to allocate and fill or copy an array
+static void fill_or_copy(double** dst, int n, const double* src, double fill_val) {
+    *dst = (double*)safe_malloc((size_t)n * sizeof(double));
+    if (src) memcpy(*dst, src, (size_t)n * sizeof(double));
+    else     for (int i = 0; i < n; ++i) (*dst)[i] = fill_val;
+}
+
+// convert dense → CSR
+static void dense_to_csr(const matrix_desc_t* desc,
+                         int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    int m = desc->m, n = desc->n;
+    double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 1e-12;
+
+    // count nnz
+    int nnz = 0;
+    for (int i = 0; i < m * n; ++i) {
+        if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    int nz = 0;
+    for (int i = 0; i < m; ++i) {
+        (*row_ptr)[i] = nz;
+        for (int j = 0; j < n; ++j) {
+            double v = desc->data.dense.A[i * n + j];
+            if (fabs(v) > tol) {
+                (*col_ind)[nz] = j;
+                (*vals)[nz] = v;
+                ++nz;
+            }
+        }
+    }
+    (*row_ptr)[m] = nz;
+    *nnz_out = nz;
+}
+
+// create an lp_problem_t from a matrix
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+) {
+    lp_problem_t* prob = (lp_problem_t*)safe_malloc(sizeof(lp_problem_t));
+
+    prob->num_variables   = A_desc->n;
+    prob->num_constraints = A_desc->m;
+
+    // handle matrix by format → convert to CSR if needed
+    switch (A_desc->fmt) {
+        case matrix_dense:
+            dense_to_csr(A_desc,
+                         &prob->constraint_matrix_row_pointers,
+                         &prob->constraint_matrix_col_indices,
+                         &prob->constraint_matrix_values,
+                         &prob->constraint_matrix_num_nonzeros);
+            break;
+
+        case matrix_csr:
+            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
+            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
+            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
+            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
+            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
+            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
+            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
+            break;
+
+        // TODO: other formats
+        default:
+            fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
+            free(prob);
+            return NULL;
+    }
+
+    // default fill values
+    prob->objective_constant = objective_constant ? *objective_constant : 0.0;
+    fill_or_copy(&prob->objective_vector, prob->num_variables, objective_c, 0.0);
+    fill_or_copy(&prob->variable_lower_bound, prob->num_variables, var_lb, 0.0);
+    fill_or_copy(&prob->variable_upper_bound, prob->num_variables, var_ub, INFINITY);
+    fill_or_copy(&prob->constraint_lower_bound, prob->num_constraints, con_lb, -INFINITY);
+    fill_or_copy(&prob->constraint_upper_bound, prob->num_constraints, con_ub, INFINITY);
+
+    return prob;
+}
+
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+) {
+    // initialize output
+    lp_solution_t out = {0};
+
+    // argument checks
+    if (!prob) {
+        fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n");
+        return out;
+    }
+
+    // prepare parameters: use defaults if not provided 
+    pdhg_parameters_t local_params;
+    if (params) {
+        local_params = *params;
+    } else {
+        set_default_parameters(&local_params);
+    }
+
+    // call optimizer
+    pdhg_solver_state_t* state = optimize(&local_params, prob);
+    if (!state) {
+        fprintf(stderr, "[interface] optimize returned NULL.\n");
+        return out;
+    }
+
+    // prepare output
+    out.n = prob->num_variables;
+    out.m = prob->num_constraints;
+    out.x = (double*)safe_malloc((size_t)out.n * sizeof(double));
+    out.y = (double*)safe_malloc((size_t)out.m * sizeof(double));
+
+    // copy solution from GPU to CPU
+    cudaError_t e1 = cudaMemcpy(out.x, state->pdhg_primal_solution,
+                                (size_t)out.n * sizeof(double), cudaMemcpyDeviceToHost);
+    cudaError_t e2 = cudaMemcpy(out.y, state->pdhg_dual_solution,
+                                (size_t)out.m * sizeof(double), cudaMemcpyDeviceToHost);
+    if (e1 != cudaSuccess || e2 != cudaSuccess) {
+        fprintf(stderr, "[interface] cudaMemcpy failed: %s / %s\n",
+                cudaGetErrorName(e1), cudaGetErrorName(e2));
+        free(out.x); out.x = NULL;
+        free(out.y); out.y = NULL;
+        pdhg_solver_state_free(state);
+        return out;
+    }
+
+    out.primal_obj = state->primal_objective_value;
+    out.dual_obj = state->dual_objective_value;
+    out.reason = state->termination_reason;
+
+    pdhg_solver_state_free(state);
+    return out;
+}
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol) {
+    if (!sol) return;
+    if (sol->x) free(sol->x);
+    if (sol->y) free(sol->y);
+    sol->n = 0;
+    sol->m = 0;
+    sol->primal_obj = 0.0;
+    sol->dual_obj   = 0.0;
+    sol->reason = TERMINATION_REASON_UNSPECIFIED;
+}
\ No newline at end of file
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index 8d408f5..d5182c9 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -16,5 +16,35 @@ limitations under the License.
 
 #pragma once
 
-// Temporary placeholder for interface API
-// Functions will be declared here later
\ No newline at end of file
+#include "struct.h"
+#include "utils.h"
+#include "io.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// create an lp_problem_t from a matrix descriptor
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+);
+
+// solve the LP problem using PDHG
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+);
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol);
+
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
\ No newline at end of file
diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu
index 6ea34e2..23cbd90 100644
--- a/cupdlpx/solver.cu
+++ b/cupdlpx/solver.cu
@@ -668,4 +668,29 @@ void rescale_info_free(rescale_info_t *info)
     free(info->var_rescale);
 
     free(info);
+}
+
+void set_default_parameters(pdhg_parameters_t *params)
+{
+    params->l_inf_ruiz_iterations = 10;
+    params->has_pock_chambolle_alpha = true;
+    params->pock_chambolle_alpha = 1.0;
+    params->bound_objective_rescaling = true;
+    params->verbose = false;
+    params->termination_evaluation_frequency = 200;
+    params->reflection_coefficient = 1.0;
+
+    params->termination_criteria.eps_optimal_relative = 1e-4;
+    params->termination_criteria.eps_feasible_relative = 1e-4;
+    params->termination_criteria.eps_infeasible = 1e-10;
+    params->termination_criteria.time_sec_limit = 3600.0;
+    params->termination_criteria.iteration_limit = INT32_MAX;
+
+    params->restart_params.artificial_restart_threshold = 0.36;
+    params->restart_params.sufficient_reduction_for_restart = 0.2;
+    params->restart_params.necessary_reduction_for_restart = 0.5;
+    params->restart_params.k_p = 0.99;
+    params->restart_params.k_i = 0.01;
+    params->restart_params.k_d = 0.0;
+    params->restart_params.i_smooth = 0.3;
 }
\ No newline at end of file
diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index 5038e48..a563dc9 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -34,6 +34,8 @@ extern "C"
 
     void pdhg_solver_state_free(pdhg_solver_state_t *state);
 
+    void set_default_parameters(pdhg_parameters_t *params);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index db01a2b..270ee84 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -196,3 +196,59 @@ typedef struct
 	double *ones_primal_d;
 	double *ones_dual_d;
 } pdhg_solver_state_t;
+
+// matrix formats
+typedef enum {
+    matrix_dense = 0,
+    matrix_csr   = 1,
+    matrix_csc   = 2,
+    matrix_coo   = 3
+} matrix_format_t;
+
+// matrix descriptor
+typedef struct {
+    int m; // num_constraints
+    int n; // num_variables
+    matrix_format_t fmt;
+
+    // treat abs(x) < zero_tolerance as zero
+    double zero_tolerance;
+
+    union {
+        struct { // Dense (row-major)
+            const double* A; // m*n
+        } dense;
+
+        struct { // CSR
+            int nnz;
+            const int*    row_ptr; 
+            const int*    col_ind;
+            const double* vals;
+        } csr;
+
+        struct { // CSC
+            int nnz;
+            const int*    col_ptr;
+            const int*    row_ind;
+            const double* vals;
+        } csc;
+
+        struct { // COO
+            int nnz;
+            const int*    row_ind;
+            const int*    col_ind; 
+            const double* vals; 
+        } coo;
+    } data;
+} matrix_desc_t;
+
+typedef struct {
+    int n;  // number of variables
+    int m;  // number of constraints
+    double* x;  // primal solution
+    double* y;  // dual solution 
+    double primal_obj; // primal objective value
+    double dual_obj; // dual objective value    
+    termination_reason_t reason; // termination reason
+} lp_solution_t;
+
diff --git a/test/test_interface.c b/test/test_interface.c
index e7d0543..6fc67fe 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -18,6 +18,23 @@ limitations under the License.
 #include <stdio.h>
 #include <math.h> 
 
+static const char* term_to_str(termination_reason_t r) {
+    switch (r) {
+        case TERMINATION_REASON_OPTIMAL:           return "OPTIMAL";
+        case TERMINATION_REASON_PRIMAL_INFEASIBLE: return "PRIMAL_INFEASIBLE";
+        case TERMINATION_REASON_DUAL_INFEASIBLE:   return "DUAL_INFEASIBLE";
+        case TERMINATION_REASON_TIME_LIMIT:        return "TIME_LIMIT";
+        case TERMINATION_REASON_ITERATION_LIMIT:   return "ITERATION_LIMIT";
+        default:                                   return "UNSPECIFIED";
+    }
+}
+
+static void print_vec(const char* name, const double* v, int n) {
+    printf("%s:", name);
+    for (int i = 0; i < n; ++i) printf(" % .6g", v[i]);
+    printf("\n");
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -32,12 +49,20 @@ int main() {
         {3.0, 2.0}
     };
 
+    // describe A using matrix_desc_t
+    matrix_desc_t A_desc;
+    A_desc.m = m;
+    A_desc.n = n;
+    A_desc.fmt = matrix_dense;
+    A_desc.zero_tolerance = 0.0;
+    A_desc.data.dense.A = &A[0][0];
+
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
 
     // l&u: constraintbounds
-    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
-    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+    double l[3] = {5.0, -INFINITY, -INFINITY};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0}; // upper bounds
 
     printf("Objective c = [");
     for (int j = 0; j < n; j++) {
@@ -59,5 +84,44 @@ int main() {
         printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
     }
 
+     lp_problem_t* prob = make_problem_from_matrix(
+        &A_desc,        // A
+        c,              // c
+        NULL,           // objective_constant
+        NULL,           // var_lb
+        NULL,           // var_ub
+        l,              // con_lb
+        u               // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed.\n");
+        return 1;
+    }
+
+    // solve problem
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+
+    // free problem
+    lp_problem_free(prob);
+
+    // check solution
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return 2;
+    }
+
+    // print solution
+    printf("Solution:\n");
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    // free solution
+    lp_solution_free(&sol);
+
     return 0;
 }

From 2911b6908d705301928ff11b61d036c08de6e87f Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 06:40:28 -0400
Subject: [PATCH 08/38] =?UTF-8?q?New=20feat:=20CSC/COO=E2=86=92CSR=20for?=
 =?UTF-8?q?=20interface?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 cupdlpx/interface.cu  | 161 +++++++++++++++++++++++++++++++++++++++++-
 test/test_interface.c | 126 +++++++++++++++++++++++++--------
 2 files changed, 254 insertions(+), 33 deletions(-)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
index d361169..3c6cd4b 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.cu
@@ -42,10 +42,12 @@ static void dense_to_csr(const matrix_desc_t* desc,
         if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
     }
 
+    // allocate
     *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
     *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
     *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
 
+    // fill
     int nz = 0;
     for (int i = 0; i < m; ++i) {
         (*row_ptr)[i] = nz;
@@ -62,6 +64,134 @@ static void dense_to_csr(const matrix_desc_t* desc,
     *nnz_out = nz;
 }
 
+// convert CSC → CSR
+static int csc_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz = desc->data.csc.nnz;
+    const int *col_ptr = desc->data.csc.col_ptr;
+    const int *row_ind = desc->data.csc.row_ind; 
+    const double *v    = desc->data.csc.vals;
+
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count entries per row
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+
+    // count nnz
+    int eff_nnz = 0;
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] CSC: row index out of range\n"); return -1; }
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            ++((*row_ptr)[ri + 1]);
+            ++eff_nnz;
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // allocate
+    *col_ind = (int*)safe_malloc((size_t)eff_nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)eff_nnz * sizeof(double));
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            int pos = next[ri]++;
+            (*col_ind)[pos] = j;
+            (*vals)[pos]    = val;
+        }
+    }
+
+    free(next);
+    *nnz_out = eff_nnz;
+    return 0;
+}
+
+// convert COO → CSR
+static int coo_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz_in = desc->data.coo.nnz;
+    const int *r = desc->data.coo.row_ind;
+    const int *c = desc->data.coo.col_ind;
+    const double *v = desc->data.coo.vals;
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count nnz
+    int nnz = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) ++nnz;
+    } else {
+        nnz = nnz_in;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    // count entries per row
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) {
+                int ri = r[k];
+                if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+                ++((*row_ptr)[ri + 1]);
+            }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+            ++((*row_ptr)[ri + 1]);
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k) {
+            if (fabs(v[k]) <= tol) continue;
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    }
+
+    free(next);
+    *nnz_out = nnz;
+    return 0;
+}
+
 // create an lp_problem_t from a matrix
 lp_problem_t* make_problem_from_matrix(
     const matrix_desc_t* A_desc,
@@ -77,7 +207,7 @@ lp_problem_t* make_problem_from_matrix(
     prob->num_variables   = A_desc->n;
     prob->num_constraints = A_desc->m;
 
-    // handle matrix by format → convert to CSR if needed
+    // handle matrix by format
     switch (A_desc->fmt) {
         case matrix_dense:
             dense_to_csr(A_desc,
@@ -97,7 +227,34 @@ lp_problem_t* make_problem_from_matrix(
             memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
             break;
 
-        // TODO: other formats
+        case matrix_csc: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] CSC->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+            break;
+        }
+
+        case matrix_coo: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (coo_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] COO->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+           break;
+        }
+            
         default:
             fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
             free(prob);
diff --git a/test/test_interface.c b/test/test_interface.c
index 6fc67fe..730592c 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -35,6 +35,45 @@ static void print_vec(const char* name, const double* v, int n) {
     printf("\n");
 }
 
+static void run_once(const char* tag,
+                     const matrix_desc_t* A_desc,
+                     const double* c, const double* l, const double* u)
+{
+    printf("\n=== %s ===\n", tag);
+
+    lp_problem_t* prob = make_problem_from_matrix(
+        A_desc,  // A
+        c,       // c
+        NULL,    // objective_constant
+        NULL,    // var_lb (defaults to 0)
+        NULL,    // var_ub (defaults to +inf)
+        l,       // con_lb
+        u        // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed for %s.\n", tag);
+        return;
+    }
+
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+    lp_problem_free(prob);
+
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return;
+    }
+
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    lp_solution_free(&sol);
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -50,12 +89,57 @@ int main() {
     };
 
     // describe A using matrix_desc_t
-    matrix_desc_t A_desc;
-    A_desc.m = m;
-    A_desc.n = n;
-    A_desc.fmt = matrix_dense;
-    A_desc.zero_tolerance = 0.0;
-    A_desc.data.dense.A = &A[0][0];
+    matrix_desc_t A_dense;
+    A_dense.m = m;
+    A_dense.n = n;
+    A_dense.fmt = matrix_dense;
+    A_dense.zero_tolerance = 0.0;
+    A_dense.data.dense.A = &A[0][0];
+
+    // A as a CSR matrix
+    static int csr_row_ptr[4] = {0, 2, 3, 5};
+    static int csr_col_ind[5] = {0, 1, 1, 0, 1};
+    static double csr_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csr;
+    A_csr.m = m; A_csr.n = n;
+    A_csr.fmt = matrix_csr;
+    A_csr.zero_tolerance = 0.0;
+    A_csr.data.csr.nnz = 5;
+    A_csr.data.csr.row_ptr = csr_row_ptr;
+    A_csr.data.csr.col_ind = csr_col_ind;
+    A_csr.data.csr.vals = csr_vals;
+
+    // A as a CSC matrix
+    static int csc_col_ptr[3] = {0, 2, 5};
+    static int csc_row_ind[5] = {0, 2, 0, 1, 2};
+    static double csc_vals[5] = {1, 3, 2, 1, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csc;
+    A_csc.m = m; A_csc.n = n;
+    A_csc.fmt = matrix_csc;
+    A_csc.zero_tolerance = 0.0;
+    A_csc.data.csc.nnz = 5;
+    A_csc.data.csc.col_ptr = csc_col_ptr;
+    A_csc.data.csc.row_ind = csc_row_ind;
+    A_csc.data.csc.vals = csc_vals;
+
+    // A as a COO matrix
+    static int coo_row_ind[5] = {0, 0, 1, 2, 2};
+    static int coo_col_ind[5] = {0, 1, 1, 0, 1};
+    static double coo_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_coo;
+    A_coo.m = m; A_coo.n = n;
+    A_coo.fmt = matrix_coo;
+    A_coo.zero_tolerance = 0.0;
+    A_coo.data.coo.nnz = 5;
+    A_coo.data.coo.row_ind = coo_row_ind;
+    A_coo.data.coo.col_ind = coo_col_ind;
+    A_coo.data.coo.vals = coo_vals;
 
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
@@ -85,7 +169,7 @@ int main() {
     }
 
      lp_problem_t* prob = make_problem_from_matrix(
-        &A_desc,        // A
+        &A_dense,       // A
         c,              // c
         NULL,           // objective_constant
         NULL,           // var_lb
@@ -98,30 +182,10 @@ int main() {
         return 1;
     }
 
-    // solve problem
-    lp_solution_t sol = solve_lp_problem(prob, NULL);
-
-    // free problem
-    lp_problem_free(prob);
-
-    // check solution
-    if (!sol.x || !sol.y) {
-        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
-                term_to_str(sol.reason));
-        lp_solution_free(&sol);
-        return 2;
-    }
-
-    // print solution
-    printf("Solution:\n");
-    printf("Termination: %s\n", term_to_str(sol.reason));
-    printf("Primal obj: %.10g\n", sol.primal_obj);
-    printf("Dual   obj: %.10g\n", sol.dual_obj);
-    print_vec("x", sol.x, sol.n);
-    print_vec("y", sol.y, sol.m);
-
-    // free solution
-    lp_solution_free(&sol);
+    run_once("Test 1: Dense Matrix", &A_dense, c, l, u);
+    run_once("Test 2: CSR Matrix", &A_csr, c, l, u);
+    run_once("Test 3: CSC Matrix", &A_csc, c, l, u);
+    run_once("Test 4: COO Matrix", &A_coo, c, l, u);
 
     return 0;
 }

From ffe947f368d4fb37daf422429b731a1875bd8f27 Mon Sep 17 00:00:00 2001
From: Zedong Peng <peng_zedong@126.com>
Date: Thu, 18 Sep 2025 14:15:13 -0400
Subject: [PATCH 09/38] fix final solution rescaling bug and crerate
 cupdlpx_result_t

---
 cupdlpx/cupdlpx.c | 44 +++++++++++-------------
 cupdlpx/solver.cu | 87 +++++++++++++++++++++++++++++++++++++++++++++--
 cupdlpx/solver.h  |  4 +--
 cupdlpx/struct.h  | 30 ++++++++++++++--
 4 files changed, 135 insertions(+), 30 deletions(-)

diff --git a/cupdlpx/cupdlpx.c b/cupdlpx/cupdlpx.c
index 25572f7..023d931 100644
--- a/cupdlpx/cupdlpx.c
+++ b/cupdlpx/cupdlpx.c
@@ -74,7 +74,7 @@ char *extract_instance_name(const char *filename)
     return instance_name;
 }
 
-void save_solution(const double *data_d, int size, const char *output_dir,
+void save_solution(const double *data, int size, const char *output_dir,
                    const char *instance_name, const char *suffix)
 {
     char *file_path = get_output_path(output_dir, instance_name, suffix);
@@ -91,20 +91,16 @@ void save_solution(const double *data_d, int size, const char *output_dir,
         return;
     }
 
-    double *solution_host = (double *)safe_malloc(size * sizeof(double));
-    CUDA_CHECK(cudaMemcpy(solution_host, data_d, size * sizeof(double), cudaMemcpyDeviceToHost));
-
     for (int i = 0; i < size; ++i)
     {
-        fprintf(outfile, "%.10g\n", solution_host[i]);
+        fprintf(outfile, "%.10g\n", data[i]);
     }
 
     fclose(outfile);
     free(file_path);
-    free(solution_host);
 }
 
-void save_solver_summary(const pdhg_solver_state_t *solver_state, const char *output_dir, const char *instance_name)
+void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir, const char *instance_name)
 {
     char *file_path = get_output_path(output_dir, instance_name, "_summary.txt");
     if (file_path == NULL)
@@ -119,17 +115,17 @@ void save_solver_summary(const pdhg_solver_state_t *solver_state, const char *ou
         free(file_path);
         return;
     }
-    fprintf(outfile, "Termination Reason: %s\n", termination_reason_tToString(solver_state->termination_reason));
-    fprintf(outfile, "Runtime (sec): %e\n", solver_state->cumulative_time_sec);
-    fprintf(outfile, "Iterations Count: %d\n", solver_state->total_count);
-    fprintf(outfile, "Primal Objective Value: %e\n", solver_state->primal_objective_value);
-    fprintf(outfile, "Dual Objective Value: %e\n", solver_state->dual_objective_value);
-    fprintf(outfile, "Absolute Primal Residual: %e\n", solver_state->absolute_primal_residual);
-    fprintf(outfile, "Relative Primal Residual: %e\n", solver_state->relative_primal_residual);
-    fprintf(outfile, "Absolute Dual Residual: %e\n", solver_state->absolute_dual_residual);
-    fprintf(outfile, "Relative Dual Residual: %e\n", solver_state->relative_dual_residual);
-    fprintf(outfile, "Absolute Objective Gap: %e\n", solver_state->objective_gap);
-    fprintf(outfile, "Relative Objective Gap: %e\n", solver_state->relative_objective_gap);
+    fprintf(outfile, "Termination Reason: %s\n", termination_reason_tToString(result->termination_reason));
+    fprintf(outfile, "Runtime (sec): %e\n", result->cumulative_time_sec);
+    fprintf(outfile, "Iterations Count: %d\n", result->total_count);
+    fprintf(outfile, "Primal Objective Value: %e\n", result->primal_objective_value);
+    fprintf(outfile, "Dual Objective Value: %e\n", result->dual_objective_value);
+    fprintf(outfile, "Absolute Primal Residual: %e\n", result->absolute_primal_residual);
+    fprintf(outfile, "Relative Primal Residual: %e\n", result->relative_primal_residual);
+    fprintf(outfile, "Absolute Dual Residual: %e\n", result->absolute_dual_residual);
+    fprintf(outfile, "Relative Dual Residual: %e\n", result->relative_dual_residual);
+    fprintf(outfile, "Absolute Objective Gap: %e\n", result->objective_gap);
+    fprintf(outfile, "Relative Objective Gap: %e\n", result->relative_objective_gap);
     fclose(outfile);
     free(file_path);
 }
@@ -251,20 +247,20 @@ int main(int argc, char *argv[])
         return 1;
     }
 
-    pdhg_solver_state_t *solver_state = optimize(&params, problem);
+    cupdlpx_result_t *result = optimize(&params, problem);
 
-    if (solver_state == NULL)
+    if (result == NULL)
     {
         fprintf(stderr, "Solver failed.\n");
     }
     else
     {
-        save_solver_summary(solver_state, output_dir, instance_name);
-        save_solution(solver_state->pdhg_primal_solution, problem->num_variables, output_dir,
+        save_solver_summary(result, output_dir, instance_name);
+        save_solution(result->primal_solution, problem->num_variables, output_dir,
                       instance_name, "_primal_solution.txt");
-        save_solution(solver_state->pdhg_dual_solution, problem->num_constraints, output_dir,
+        save_solution(result->dual_solution, problem->num_constraints, output_dir,
                       instance_name, "_dual_solution.txt");
-        pdhg_solver_state_free(solver_state);
+        cupdlpx_result_free(result);
     }
 
     lp_problem_free(problem);
diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu
index 6ea34e2..9dc0b45 100644
--- a/cupdlpx/solver.cu
+++ b/cupdlpx/solver.cu
@@ -45,6 +45,10 @@ __global__ void halpern_update_kernel(
     const double *initial_primal, double *current_primal, const double *reflected_primal,
     const double *initial_dual, double *current_dual, const double *reflected_dual,
     int n_vars, int n_cons, double weight, double reflection_coeff);
+__global__ void rescale_solution_kernel(
+    double *primal_solution, double *dual_solution, const double *variable_rescaling, const double *constraint_rescaling,
+    const double objective_vector_rescaling, const double constraint_bound_rescaling,
+    int n_vars, int n_cons);
 __global__ void compute_delta_solution_kernel(
     const double *initial_primal, const double *pdhg_primal, double *delta_primal,
     const double *initial_dual, const double *pdhg_dual, double *delta_dual,
@@ -52,15 +56,18 @@ __global__ void compute_delta_solution_kernel(
 static void compute_next_pdhg_primal_solution(pdhg_solver_state_t *state);
 static void compute_next_pdhg_dual_solution(pdhg_solver_state_t *state);
 static void halpern_update(pdhg_solver_state_t *state, double reflection_coefficient);
+static void rescale_solution(pdhg_solver_state_t *state);
+static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state);
 static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *cudaMemsetParams);
 static void initialize_step_size_and_primal_weight(pdhg_solver_state_t *state, const pdhg_parameters_t *params);
 static pdhg_solver_state_t *initialize_solver_state(
     const lp_problem_t *original_problem,
     const rescale_info_t *rescale_info);
 static void compute_fixed_point_error(pdhg_solver_state_t *state);
+void pdhg_solver_state_free(pdhg_solver_state_t *state);
 void rescale_info_free(rescale_info_t *info);
 
-pdhg_solver_state_t *optimize(const pdhg_parameters_t *params, const lp_problem_t *original_problem)
+cupdlpx_result_t *optimize(const pdhg_parameters_t *params, const lp_problem_t *original_problem)
 {
     print_initial_info(params, original_problem);
     rescale_info_t *rescale_info = rescale_problem(params, original_problem);
@@ -115,7 +122,9 @@ pdhg_solver_state_t *optimize(const pdhg_parameters_t *params, const lp_problem_
     }
 
     pdhg_final_log(state, params->verbose, state->termination_reason);
-    return state;
+    cupdlpx_result_t *results = create_result_from_state(state);
+    pdhg_solver_state_free(state);
+    return results;
 }
 
 static pdhg_solver_state_t *initialize_solver_state(
@@ -146,6 +155,8 @@ static pdhg_solver_state_t *initialize_solver_state(
 
     state->termination_reason = TERMINATION_REASON_UNSPECIFIED;
 
+    state->rescaling_time_sec = rescale_info->rescaling_time_sec;
+
 #define ALLOC_AND_COPY(dest, src, bytes)  \
     CUDA_CHECK(cudaMalloc(&dest, bytes)); \
     CUDA_CHECK(cudaMemcpy(dest, src, bytes, cudaMemcpyHostToDevice));
@@ -387,6 +398,23 @@ __global__ void halpern_update_kernel(
     }
 }
 
+__global__ void rescale_solution_kernel(
+    double *primal_solution, double *dual_solution, const double *variable_rescaling, const double *constraint_rescaling,
+    const double objective_vector_rescaling, const double constraint_bound_rescaling,
+    int n_vars, int n_cons)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < n_vars)
+    {
+        primal_solution[i] = primal_solution[i] / variable_rescaling[i] / constraint_bound_rescaling;
+    }
+    else if (i < n_vars + n_cons)
+    {
+        int idx = i - n_vars;
+        dual_solution[idx] = dual_solution[idx] / constraint_rescaling[idx] / objective_vector_rescaling;
+    }
+}
+
 __global__ void compute_delta_solution_kernel(
     const double *initial_primal, const double *pdhg_primal, double *delta_primal,
     const double *initial_dual, const double *pdhg_dual, double *delta_dual,
@@ -468,6 +496,15 @@ static void halpern_update(pdhg_solver_state_t *state, double reflection_coeffic
         state->num_variables, state->num_constraints, weight, reflection_coefficient);
 }
 
+static void rescale_solution(pdhg_solver_state_t *state)
+{
+    rescale_solution_kernel<<<state->num_blocks_primal_dual, THREADS_PER_BLOCK>>>(
+        state->pdhg_primal_solution, state->pdhg_dual_solution,
+        state->variable_rescaling, state->constraint_rescaling,
+        state->objective_vector_rescaling, state->constraint_bound_rescaling,
+        state->num_variables, state->num_constraints);
+}
+
 static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params)
 {
     compute_delta_solution_kernel<<<state->num_blocks_primal_dual, THREADS_PER_BLOCK>>>(
@@ -668,4 +705,50 @@ void rescale_info_free(rescale_info_t *info)
     free(info->var_rescale);
 
     free(info);
+}
+
+void cupdlpx_result_free(cupdlpx_result_t *results)
+{
+    if (results == NULL)
+    {
+        return;
+    }
+
+    free(results->primal_solution);
+    free(results->dual_solution);
+    free(results);
+}
+
+static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state)
+{
+    cupdlpx_result_t *results = (cupdlpx_result_t *)safe_calloc(1, sizeof(cupdlpx_result_t));
+
+    rescale_solution(state);
+
+    results->primal_solution = (double *)safe_malloc(state->num_variables * sizeof(double));
+    results->dual_solution = (double *)safe_malloc(state->num_constraints * sizeof(double));
+
+    CUDA_CHECK(cudaMemcpy(results->primal_solution, state->pdhg_primal_solution, state->num_variables * sizeof(double), cudaMemcpyDeviceToHost));
+    CUDA_CHECK(cudaMemcpy(results->dual_solution, state->pdhg_dual_solution, state->num_constraints * sizeof(double), cudaMemcpyDeviceToHost));
+
+    results->num_variables = state->num_variables;
+    results->num_constraints = state->num_constraints;
+    results->total_count = state->total_count;
+    results->rescaling_time_sec = state->rescaling_time_sec;
+    results->cumulative_time_sec = state->cumulative_time_sec;
+    results->absolute_primal_residual = state->absolute_primal_residual;
+    results->relative_primal_residual = state->relative_primal_residual;
+    results->absolute_dual_residual = state->absolute_dual_residual;
+    results->relative_dual_residual = state->relative_dual_residual;
+    results->primal_objective_value = state->primal_objective_value;
+    results->dual_objective_value = state->dual_objective_value;
+    results->objective_gap = state->objective_gap;
+    results->relative_objective_gap = state->relative_objective_gap;
+    results->max_primal_ray_infeasibility = state->max_primal_ray_infeasibility;
+    results->max_dual_ray_infeasibility = state->max_dual_ray_infeasibility;
+    results->primal_ray_linear_objective = state->primal_ray_linear_objective;
+    results->dual_ray_objective = state->dual_ray_objective;
+    results->termination_reason = state->termination_reason;
+
+    return results;
 }
\ No newline at end of file
diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index 5038e48..a3b1750 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -28,11 +28,11 @@ extern "C"
 {
 #endif
 
-    pdhg_solver_state_t *optimize(
+    cupdlpx_result_t *optimize(
         const pdhg_parameters_t *params,
         const lp_problem_t *original_problem);
 
-    void pdhg_solver_state_free(pdhg_solver_state_t *state);
+    void cupdlpx_result_free(cupdlpx_result_t *results);
 
 #ifdef __cplusplus
 }
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index db01a2b..e5adc63 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -150,8 +150,6 @@ typedef struct
 	double *primal_slack;
 	double *dual_slack;
 	double rescaling_time_sec;
-	double gpu_to_cpu_time_sec;
-	double basic_time_sec;
 	double cumulative_time_sec;
 
 	double *primal_residual;
@@ -196,3 +194,31 @@ typedef struct
 	double *ones_primal_d;
 	double *ones_dual_d;
 } pdhg_solver_state_t;
+
+typedef struct
+{
+	int num_variables;
+	int num_constraints;
+
+	double *primal_solution;
+	double *dual_solution;
+	
+	int total_count;
+	double rescaling_time_sec;
+	double cumulative_time_sec;
+
+	double absolute_primal_residual;
+	double relative_primal_residual;
+	double absolute_dual_residual;
+	double relative_dual_residual;
+	double primal_objective_value;
+	double dual_objective_value;
+	double objective_gap;
+	double relative_objective_gap;
+	double max_primal_ray_infeasibility;
+	double max_dual_ray_infeasibility;
+	double primal_ray_linear_objective;
+	double dual_ray_objective;
+	termination_reason_t termination_reason;
+
+} cupdlpx_result_t;
\ No newline at end of file

From cda216efeecc99f335dd6d7ab1b37efbbd00a8d7 Mon Sep 17 00:00:00 2001
From: Zedong Peng <peng_zedong@126.com>
Date: Thu, 18 Sep 2025 14:26:13 -0400
Subject: [PATCH 10/38] fix io variable num bug

---
 cupdlpx/io.c | 40 +++++++++++++++++++++++++++++++++-------
 1 file changed, 33 insertions(+), 7 deletions(-)

diff --git a/cupdlpx/io.c b/cupdlpx/io.c
index 273c836..4efac82 100644
--- a/cupdlpx/io.c
+++ b/cupdlpx/io.c
@@ -325,6 +325,7 @@ typedef struct
     size_t constraint_capacity;
 
     char *objective_row_name;
+    char *current_col_name;
     double objective_constant;
     bool is_maximize;
     int error_flag;
@@ -649,18 +650,42 @@ static int finalize_rows(MpsParserState *state)
 
 static int parse_columns_section(MpsParserState *state, char **tokens, int n_tokens)
 {
-    if (n_tokens < 3)
+    if (n_tokens < 2)
         return 0;
-    const char *col_name = tokens[0];
 
-    if (!ensure_column_capacity(state))
-        return -1;
+    if (n_tokens >= 2 && strcmp(tokens[1], "'MARKER'") == 0)
+    {
+        return 0;
+    }
+
+    const char *col_name = NULL;
+    int pair_start_index;
+
+    if (n_tokens % 2 != 0)
+    {
+        free(state->current_col_name);
+        state->current_col_name = strdup(tokens[0]);
+        if (!state->current_col_name) return -1;
+        
+        col_name = state->current_col_name;
+        pair_start_index = 1;
+    }
+    else
+    {
+        if (!state->current_col_name) {
+             fprintf(stderr, "ERROR: Column data found before any column name was defined.\n");
+             return -1;
+        }
+        col_name = state->current_col_name;
+        pair_start_index = 0;
+    }
+
+    if (!ensure_column_capacity(state)) return -1;
 
     int col_idx = namemap_put(&state->col_map, col_name);
-    if (col_idx == -1)
-        return -1;
+    if (col_idx == -1) return -1;
 
-    for (int i = 1; i + 1 < n_tokens; i += 2)
+    for (int i = pair_start_index; i + 1 < n_tokens; i += 2)
     {
         const char *row_name = tokens[i];
         double value = atof(tokens[i + 1]);
@@ -866,6 +891,7 @@ static void free_parser_state(MpsParserState *state)
     free(state->constraint_lower_bounds);
     free(state->constraint_upper_bounds);
     free(state->objective_row_name);
+    free(state->current_col_name);
 }
 
 void lp_problem_free(lp_problem_t *L)

From f19f32e33538111f4b93fdf2fa0fafc3f5ccc243 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:53:07 -0400
Subject: [PATCH 11/38] Create .gitignore

---
 .gitignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 9c13904..a84db7a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,3 @@
-
 # Prerequisites
 *.d
 
@@ -58,3 +57,4 @@ dkms.conf
 # ignored files
 build/*
 test/*
+

From 2b0b9a88bd680724f1db6ec22a23befad72a43b8 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:55:06 -0400
Subject: [PATCH 12/38] Update .gitignore

---
 .gitignore | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index a84db7a..52defc1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,4 +57,3 @@ dkms.conf
 # ignored files
 build/*
 test/*
-

From 2915bfe995e2e8687a4d9598ad7eddd2510326cb Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 02:03:02 -0400
Subject: [PATCH 13/38] Init interface

---
 .gitignore          |  1 +
 Makefile            | 73 +++++++++++++++++++++++++++++++++++++++++++--
 cupdlpx/interface.c | 18 +++++++++++
 cupdlpx/interface.h | 20 +++++++++++++
 4 files changed, 110 insertions(+), 2 deletions(-)
 create mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.h

diff --git a/.gitignore b/.gitignore
index 52defc1..654154c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,3 +57,4 @@ dkms.conf
 # ignored files
 build/*
 test/*
+/.vscode
diff --git a/Makefile b/Makefile
index 3fccf67..74c0c5b 100644
--- a/Makefile
+++ b/Makefile
@@ -16,6 +16,7 @@ NVCCFLAGS = -I. -I$(CUDA_HOME)/include -O3 -g -gencode arch=compute_90,code=sm_9
 # LDFLAGS for the linker
 LDFLAGS = -L$(CONDA_PREFIX)/lib -L$(CUDA_HOME)/lib64 -lcudart -lcusparse -lcublas -lz -lm
 
+# Source discovery (exclude the debug main)
 C_SOURCES = $(filter-out $(SRC_DIR)/cupdlpx.c, $(wildcard $(SRC_DIR)/*.c))
 CU_SOURCES = $(wildcard $(SRC_DIR)/*.cu)
 
@@ -25,22 +26,41 @@ OBJECTS = $(C_OBJECTS) $(CU_OBJECTS)
 
 TARGET_LIB = $(BUILD_DIR)/libcupdlpx.a
 
+# Debug executable (optional)
 DEBUG_SRC = $(SRC_DIR)/cupdlpx.c
 DEBUG_EXEC = $(BUILD_DIR)/cupdlpx
 
-.PHONY: all clean build
+# Tests auto-discovery
+TEST_DIR := ./test
+TEST_BUILD_DIR := $(BUILD_DIR)/tests
 
+TEST_CU_SOURCES := $(wildcard $(TEST_DIR)/*.cu)
+TEST_C_SOURCES := $(wildcard $(TEST_DIR)/*.c)
+
+# Each test source becomes an executable at build/tests/<basename>
+TEST_EXEC_CU := $(patsubst $(TEST_DIR)/%.cu,$(TEST_BUILD_DIR)/%,$(TEST_CU_SOURCES))
+TEST_EXEC_C := $(patsubst $(TEST_DIR)/%.c,$(TEST_BUILD_DIR)/%,$(TEST_C_SOURCES))
+
+# Phony targets
+.PHONY: all clean build tests test run-tests run-test clean-tests
+
+# Default: build the static library
 all: $(TARGET_LIB)
 
+# Archive all objects into the static library
 $(TARGET_LIB): $(OBJECTS)
 	@echo "Archiving objects into $(TARGET_LIB)..."
+	@mkdir -p $(BUILD_DIR)
 	@ar rcs $@ $^
 
+# Build the debug executable (links the library with cupdlpx.c main)
 build: $(DEBUG_EXEC)
+
 $(DEBUG_EXEC): $(DEBUG_SRC) $(TARGET_LIB)
 	@echo "Building debug executable..."
 	@$(LINKER) $(NVCCFLAGS) $(DEBUG_SRC) -o $(DEBUG_EXEC) $(TARGET_LIB) $(LDFLAGS)
 
+# Pattern rules for objects
 $(BUILD_DIR)/%.o: $(SRC_DIR)/%.c
 	@mkdir -p $(BUILD_DIR)
 	@echo "Compiling $< -> $@..."
@@ -51,6 +71,55 @@ $(BUILD_DIR)/%.o: $(SRC_DIR)/%.cu
 	@echo "Compiling $< -> $@..."
 	@$(NVCC) $(NVCCFLAGS) -c $< -o $@
 
+# Build all tests discovered under test/
+test: tests
+tests: $(TARGET_LIB) $(TEST_EXEC_CU) $(TEST_EXEC_C)
+	@echo "All tests built under $(TEST_BUILD_DIR)/"
+
+# Run all tests one by one
+run-tests: tests
+	@echo "Running all tests..."
+	@set -e; \
+	for t in $(TEST_EXEC_CU) $(TEST_EXEC_C); do \
+	  echo "=== $$t ==="; \
+	  "$$t" || exit $$?; \
+	  echo; \
+	done
+
+# Run a single test by basename: make run-test name=<basename>
+run-test: tests
+	@if [ -z "$(name)" ]; then \
+	  echo "Usage: make run-test name=<basename-of-test-file>"; exit 2; \
+	fi
+	@if [ -x "$(TEST_BUILD_DIR)/$(name)" ]; then \
+	  echo "=== $(TEST_BUILD_DIR)/$(name) ==="; \
+	  "$(TEST_BUILD_DIR)/$(name)"; \
+	else \
+	  echo "Executable not found: $(TEST_BUILD_DIR)/$(name)"; \
+	  echo "Did you 'make tests' and is there a test source named '$(TEST_DIR)/$(name).c(u)'?"; \
+	  exit 1; \
+	fi
+
+# Build rule for CUDA tests: compile and link directly with nvcc
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.cu $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building CUDA test $< -> $@..."
+	@$(LINKER) $(NVCCFLAGS) -I$(SRC_DIR) $< -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+# Build rule for C tests: compile with gcc, link with nvcc to get CUDA libs
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.c $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building C test $< -> $@..."
+	@$(CC) $(CFLAGS) -I$(SRC_DIR) -c $< -o $(TEST_BUILD_DIR)/$*.o
+	@$(LINKER) $(NVCCFLAGS) $(TEST_BUILD_DIR)/$*.o -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+
+# Cleaning
+clean-tests:
+	@echo "Cleaning test executables..."
+	@rm -rf $(TEST_BUILD_DIR)
+
 clean:
 	@echo "Cleaning up..."
-	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
\ No newline at end of file
+	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
+	@rm -rf $(TEST_BUILD_DIR)
\ No newline at end of file
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
new file mode 100644
index 0000000..6313091
--- /dev/null
+++ b/cupdlpx/interface.c
@@ -0,0 +1,18 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
new file mode 100644
index 0000000..8d408f5
--- /dev/null
+++ b/cupdlpx/interface.h
@@ -0,0 +1,20 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#pragma once
+
+// Temporary placeholder for interface API
+// Functions will be declared here later
\ No newline at end of file

From e6391cfde0eade4ffa0a92e50684153d228754a3 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 16:14:22 -0400
Subject: [PATCH 14/38] New test: LP model

---
 test/test_interface.c | 63 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100644 test/test_interface.c

diff --git a/test/test_interface.c b/test/test_interface.c
new file mode 100644
index 0000000..e7d0543
--- /dev/null
+++ b/test/test_interface.c
@@ -0,0 +1,63 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "cupdlpx/interface.h"
+#include <stdio.h>
+#include <math.h> 
+
+int main() {
+    // Example: min c^T x
+    // s.t. l <= A x <= u, x >= 0
+
+    int m = 3; // number of constraints
+    int n = 2; // number of variables
+
+    // A as a dense matrix
+    double A[3][2] = {
+        {1.0, 2.0},
+        {0.0, 1.0},
+        {3.0, 2.0}
+    };
+
+    // c: objective coefficients
+    double c[2] = {1.0, 1.0};
+
+    // l&u: constraintbounds
+    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+
+    printf("Objective c = [");
+    for (int j = 0; j < n; j++) {
+        printf(" %g", c[j]);
+    }
+    printf(" ]\n");
+
+    printf("Matrix A:\n");
+    for (int i = 0; i < m; i++) {
+        printf("  [");
+        for (int j = 0; j < n; j++) {
+            printf(" %g", A[i][j]);
+        }
+        printf(" ]\n");
+    }
+
+    printf("Constraint bounds (l <= A x <= b):\n");
+    for (int i = 0; i < m; i++) {
+        printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
+    }
+
+    return 0;
+}

From 29f43e192878e1d2e039646b15a5f72fc254d356 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 05:59:37 -0400
Subject: [PATCH 15/38] New feat: interface

---
 cupdlpx/cupdlpx.c     |  25 ------
 cupdlpx/interface.c   |  18 -----
 cupdlpx/interface.cu  | 184 ++++++++++++++++++++++++++++++++++++++++++
 cupdlpx/interface.h   |  34 +++++++-
 cupdlpx/solver.cu     |  26 ++++++
 cupdlpx/solver.h      |   2 +
 cupdlpx/struct.h      |  57 ++++++++++++-
 test/test_interface.c |  68 +++++++++++++++-
 8 files changed, 366 insertions(+), 48 deletions(-)
 delete mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.cu

diff --git a/cupdlpx/cupdlpx.c b/cupdlpx/cupdlpx.c
index 023d931..de26971 100644
--- a/cupdlpx/cupdlpx.c
+++ b/cupdlpx/cupdlpx.c
@@ -130,31 +130,6 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir,
     free(file_path);
 }
 
-void set_default_parameters(pdhg_parameters_t *params)
-{
-    params->l_inf_ruiz_iterations = 10;
-    params->has_pock_chambolle_alpha = true;
-    params->pock_chambolle_alpha = 1.0;
-    params->bound_objective_rescaling = true;
-    params->verbose = false;
-    params->termination_evaluation_frequency = 200;
-    params->reflection_coefficient = 1.0;
-
-    params->termination_criteria.eps_optimal_relative = 1e-4;
-    params->termination_criteria.eps_feasible_relative = 1e-4;
-    params->termination_criteria.eps_infeasible = 1e-10;
-    params->termination_criteria.time_sec_limit = 3600.0;
-    params->termination_criteria.iteration_limit = INT32_MAX;
-
-    params->restart_params.artificial_restart_threshold = 0.36;
-    params->restart_params.sufficient_reduction_for_restart = 0.2;
-    params->restart_params.necessary_reduction_for_restart = 0.5;
-    params->restart_params.k_p = 0.99;
-    params->restart_params.k_i = 0.01;
-    params->restart_params.k_d = 0.0;
-    params->restart_params.i_smooth = 0.3;
-}
-
 void print_usage(const char *prog_name)
 {
     fprintf(stderr, "Usage: %s [OPTIONS] <mps_file> <output_dir>\n\n", prog_name);
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
deleted file mode 100644
index 6313091..0000000
--- a/cupdlpx/interface.c
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
-Copyright 2025 Haihao Lu
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "interface.h"
-
diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
new file mode 100644
index 0000000..d361169
--- /dev/null
+++ b/cupdlpx/interface.cu
@@ -0,0 +1,184 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+#include "solver.h"
+#include "utils.h"
+#include <cuda_runtime.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+// helper function to allocate and fill or copy an array
+static void fill_or_copy(double** dst, int n, const double* src, double fill_val) {
+    *dst = (double*)safe_malloc((size_t)n * sizeof(double));
+    if (src) memcpy(*dst, src, (size_t)n * sizeof(double));
+    else     for (int i = 0; i < n; ++i) (*dst)[i] = fill_val;
+}
+
+// convert dense → CSR
+static void dense_to_csr(const matrix_desc_t* desc,
+                         int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    int m = desc->m, n = desc->n;
+    double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 1e-12;
+
+    // count nnz
+    int nnz = 0;
+    for (int i = 0; i < m * n; ++i) {
+        if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    int nz = 0;
+    for (int i = 0; i < m; ++i) {
+        (*row_ptr)[i] = nz;
+        for (int j = 0; j < n; ++j) {
+            double v = desc->data.dense.A[i * n + j];
+            if (fabs(v) > tol) {
+                (*col_ind)[nz] = j;
+                (*vals)[nz] = v;
+                ++nz;
+            }
+        }
+    }
+    (*row_ptr)[m] = nz;
+    *nnz_out = nz;
+}
+
+// create an lp_problem_t from a matrix
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+) {
+    lp_problem_t* prob = (lp_problem_t*)safe_malloc(sizeof(lp_problem_t));
+
+    prob->num_variables   = A_desc->n;
+    prob->num_constraints = A_desc->m;
+
+    // handle matrix by format → convert to CSR if needed
+    switch (A_desc->fmt) {
+        case matrix_dense:
+            dense_to_csr(A_desc,
+                         &prob->constraint_matrix_row_pointers,
+                         &prob->constraint_matrix_col_indices,
+                         &prob->constraint_matrix_values,
+                         &prob->constraint_matrix_num_nonzeros);
+            break;
+
+        case matrix_csr:
+            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
+            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
+            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
+            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
+            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
+            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
+            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
+            break;
+
+        // TODO: other formats
+        default:
+            fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
+            free(prob);
+            return NULL;
+    }
+
+    // default fill values
+    prob->objective_constant = objective_constant ? *objective_constant : 0.0;
+    fill_or_copy(&prob->objective_vector, prob->num_variables, objective_c, 0.0);
+    fill_or_copy(&prob->variable_lower_bound, prob->num_variables, var_lb, 0.0);
+    fill_or_copy(&prob->variable_upper_bound, prob->num_variables, var_ub, INFINITY);
+    fill_or_copy(&prob->constraint_lower_bound, prob->num_constraints, con_lb, -INFINITY);
+    fill_or_copy(&prob->constraint_upper_bound, prob->num_constraints, con_ub, INFINITY);
+
+    return prob;
+}
+
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+) {
+    // initialize output
+    lp_solution_t out = {0};
+
+    // argument checks
+    if (!prob) {
+        fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n");
+        return out;
+    }
+
+    // prepare parameters: use defaults if not provided 
+    pdhg_parameters_t local_params;
+    if (params) {
+        local_params = *params;
+    } else {
+        set_default_parameters(&local_params);
+    }
+
+    // call optimizer
+    pdhg_solver_state_t* state = optimize(&local_params, prob);
+    if (!state) {
+        fprintf(stderr, "[interface] optimize returned NULL.\n");
+        return out;
+    }
+
+    // prepare output
+    out.n = prob->num_variables;
+    out.m = prob->num_constraints;
+    out.x = (double*)safe_malloc((size_t)out.n * sizeof(double));
+    out.y = (double*)safe_malloc((size_t)out.m * sizeof(double));
+
+    // copy solution from GPU to CPU
+    cudaError_t e1 = cudaMemcpy(out.x, state->pdhg_primal_solution,
+                                (size_t)out.n * sizeof(double), cudaMemcpyDeviceToHost);
+    cudaError_t e2 = cudaMemcpy(out.y, state->pdhg_dual_solution,
+                                (size_t)out.m * sizeof(double), cudaMemcpyDeviceToHost);
+    if (e1 != cudaSuccess || e2 != cudaSuccess) {
+        fprintf(stderr, "[interface] cudaMemcpy failed: %s / %s\n",
+                cudaGetErrorName(e1), cudaGetErrorName(e2));
+        free(out.x); out.x = NULL;
+        free(out.y); out.y = NULL;
+        pdhg_solver_state_free(state);
+        return out;
+    }
+
+    out.primal_obj = state->primal_objective_value;
+    out.dual_obj = state->dual_objective_value;
+    out.reason = state->termination_reason;
+
+    pdhg_solver_state_free(state);
+    return out;
+}
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol) {
+    if (!sol) return;
+    if (sol->x) free(sol->x);
+    if (sol->y) free(sol->y);
+    sol->n = 0;
+    sol->m = 0;
+    sol->primal_obj = 0.0;
+    sol->dual_obj   = 0.0;
+    sol->reason = TERMINATION_REASON_UNSPECIFIED;
+}
\ No newline at end of file
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index 8d408f5..d5182c9 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -16,5 +16,35 @@ limitations under the License.
 
 #pragma once
 
-// Temporary placeholder for interface API
-// Functions will be declared here later
\ No newline at end of file
+#include "struct.h"
+#include "utils.h"
+#include "io.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// create an lp_problem_t from a matrix descriptor
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+);
+
+// solve the LP problem using PDHG
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+);
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol);
+
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
\ No newline at end of file
diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu
index 9dc0b45..7d6e762 100644
--- a/cupdlpx/solver.cu
+++ b/cupdlpx/solver.cu
@@ -751,4 +751,30 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state)
     results->termination_reason = state->termination_reason;
 
     return results;
+}
+
+void set_default_parameters(pdhg_parameters_t *params)
+{
+    params->l_inf_ruiz_iterations = 10;
+    params->has_pock_chambolle_alpha = true;
+    params->pock_chambolle_alpha = 1.0;
+    params->bound_objective_rescaling = true;
+    params->verbose = false;
+    params->termination_evaluation_frequency = 200;
+    params->reflection_coefficient = 1.0;
+
+    params->termination_criteria.eps_optimal_relative = 1e-4;
+    params->termination_criteria.eps_feasible_relative = 1e-4;
+    params->termination_criteria.eps_infeasible = 1e-10;
+    params->termination_criteria.time_sec_limit = 3600.0;
+    params->termination_criteria.iteration_limit = INT32_MAX;
+
+    params->restart_params.artificial_restart_threshold = 0.36;
+    params->restart_params.sufficient_reduction_for_restart = 0.2;
+    params->restart_params.necessary_reduction_for_restart = 0.5;
+    params->restart_params.k_p = 0.99;
+    params->restart_params.k_i = 0.01;
+    params->restart_params.k_d = 0.0;
+    params->restart_params.i_smooth = 0.3;
+
 }
\ No newline at end of file
diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index a3b1750..56da0bd 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -34,6 +34,8 @@ extern "C"
 
     void cupdlpx_result_free(cupdlpx_result_t *results);
 
+    void set_default_parameters(pdhg_parameters_t *params);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index e5adc63..12daeef 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -221,4 +221,59 @@ typedef struct
 	double dual_ray_objective;
 	termination_reason_t termination_reason;
 
-} cupdlpx_result_t;
\ No newline at end of file
+} cupdlpx_result_t;
+
+// matrix formats
+typedef enum {
+    matrix_dense = 0,
+    matrix_csr   = 1,
+    matrix_csc   = 2,
+    matrix_coo   = 3
+} matrix_format_t;
+
+// matrix descriptor
+typedef struct {
+    int m; // num_constraints
+    int n; // num_variables
+    matrix_format_t fmt;
+
+    // treat abs(x) < zero_tolerance as zero
+    double zero_tolerance;
+
+    union {
+        struct { // Dense (row-major)
+            const double* A; // m*n
+        } dense;
+
+        struct { // CSR
+            int nnz;
+            const int*    row_ptr; 
+            const int*    col_ind;
+            const double* vals;
+        } csr;
+
+        struct { // CSC
+            int nnz;
+            const int*    col_ptr;
+            const int*    row_ind;
+            const double* vals;
+        } csc;
+
+        struct { // COO
+            int nnz;
+            const int*    row_ind;
+            const int*    col_ind; 
+            const double* vals; 
+        } coo;
+    } data;
+} matrix_desc_t;
+
+typedef struct {
+    int n;  // number of variables
+    int m;  // number of constraints
+    double* x;  // primal solution
+    double* y;  // dual solution 
+    double primal_obj; // primal objective value
+    double dual_obj; // dual objective value    
+    termination_reason_t reason; // termination reason
+} lp_solution_t;
diff --git a/test/test_interface.c b/test/test_interface.c
index e7d0543..6fc67fe 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -18,6 +18,23 @@ limitations under the License.
 #include <stdio.h>
 #include <math.h> 
 
+static const char* term_to_str(termination_reason_t r) {
+    switch (r) {
+        case TERMINATION_REASON_OPTIMAL:           return "OPTIMAL";
+        case TERMINATION_REASON_PRIMAL_INFEASIBLE: return "PRIMAL_INFEASIBLE";
+        case TERMINATION_REASON_DUAL_INFEASIBLE:   return "DUAL_INFEASIBLE";
+        case TERMINATION_REASON_TIME_LIMIT:        return "TIME_LIMIT";
+        case TERMINATION_REASON_ITERATION_LIMIT:   return "ITERATION_LIMIT";
+        default:                                   return "UNSPECIFIED";
+    }
+}
+
+static void print_vec(const char* name, const double* v, int n) {
+    printf("%s:", name);
+    for (int i = 0; i < n; ++i) printf(" % .6g", v[i]);
+    printf("\n");
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -32,12 +49,20 @@ int main() {
         {3.0, 2.0}
     };
 
+    // describe A using matrix_desc_t
+    matrix_desc_t A_desc;
+    A_desc.m = m;
+    A_desc.n = n;
+    A_desc.fmt = matrix_dense;
+    A_desc.zero_tolerance = 0.0;
+    A_desc.data.dense.A = &A[0][0];
+
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
 
     // l&u: constraintbounds
-    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
-    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+    double l[3] = {5.0, -INFINITY, -INFINITY};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0}; // upper bounds
 
     printf("Objective c = [");
     for (int j = 0; j < n; j++) {
@@ -59,5 +84,44 @@ int main() {
         printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
     }
 
+     lp_problem_t* prob = make_problem_from_matrix(
+        &A_desc,        // A
+        c,              // c
+        NULL,           // objective_constant
+        NULL,           // var_lb
+        NULL,           // var_ub
+        l,              // con_lb
+        u               // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed.\n");
+        return 1;
+    }
+
+    // solve problem
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+
+    // free problem
+    lp_problem_free(prob);
+
+    // check solution
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return 2;
+    }
+
+    // print solution
+    printf("Solution:\n");
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    // free solution
+    lp_solution_free(&sol);
+
     return 0;
 }

From f144d1d4e96d432259c6ca638726b59efce50196 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 06:40:28 -0400
Subject: [PATCH 16/38] =?UTF-8?q?New=20feat:=20CSC/COO=E2=86=92CSR=20for?=
 =?UTF-8?q?=20interface?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 cupdlpx/interface.cu  | 161 +++++++++++++++++++++++++++++++++++++++++-
 test/test_interface.c | 126 +++++++++++++++++++++++++--------
 2 files changed, 254 insertions(+), 33 deletions(-)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
index d361169..3c6cd4b 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.cu
@@ -42,10 +42,12 @@ static void dense_to_csr(const matrix_desc_t* desc,
         if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
     }
 
+    // allocate
     *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
     *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
     *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
 
+    // fill
     int nz = 0;
     for (int i = 0; i < m; ++i) {
         (*row_ptr)[i] = nz;
@@ -62,6 +64,134 @@ static void dense_to_csr(const matrix_desc_t* desc,
     *nnz_out = nz;
 }
 
+// convert CSC → CSR
+static int csc_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz = desc->data.csc.nnz;
+    const int *col_ptr = desc->data.csc.col_ptr;
+    const int *row_ind = desc->data.csc.row_ind; 
+    const double *v    = desc->data.csc.vals;
+
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count entries per row
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+
+    // count nnz
+    int eff_nnz = 0;
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] CSC: row index out of range\n"); return -1; }
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            ++((*row_ptr)[ri + 1]);
+            ++eff_nnz;
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // allocate
+    *col_ind = (int*)safe_malloc((size_t)eff_nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)eff_nnz * sizeof(double));
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            int pos = next[ri]++;
+            (*col_ind)[pos] = j;
+            (*vals)[pos]    = val;
+        }
+    }
+
+    free(next);
+    *nnz_out = eff_nnz;
+    return 0;
+}
+
+// convert COO → CSR
+static int coo_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz_in = desc->data.coo.nnz;
+    const int *r = desc->data.coo.row_ind;
+    const int *c = desc->data.coo.col_ind;
+    const double *v = desc->data.coo.vals;
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count nnz
+    int nnz = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) ++nnz;
+    } else {
+        nnz = nnz_in;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    // count entries per row
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) {
+                int ri = r[k];
+                if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+                ++((*row_ptr)[ri + 1]);
+            }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+            ++((*row_ptr)[ri + 1]);
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k) {
+            if (fabs(v[k]) <= tol) continue;
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    }
+
+    free(next);
+    *nnz_out = nnz;
+    return 0;
+}
+
 // create an lp_problem_t from a matrix
 lp_problem_t* make_problem_from_matrix(
     const matrix_desc_t* A_desc,
@@ -77,7 +207,7 @@ lp_problem_t* make_problem_from_matrix(
     prob->num_variables   = A_desc->n;
     prob->num_constraints = A_desc->m;
 
-    // handle matrix by format → convert to CSR if needed
+    // handle matrix by format
     switch (A_desc->fmt) {
         case matrix_dense:
             dense_to_csr(A_desc,
@@ -97,7 +227,34 @@ lp_problem_t* make_problem_from_matrix(
             memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
             break;
 
-        // TODO: other formats
+        case matrix_csc: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] CSC->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+            break;
+        }
+
+        case matrix_coo: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (coo_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] COO->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+           break;
+        }
+            
         default:
             fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
             free(prob);
diff --git a/test/test_interface.c b/test/test_interface.c
index 6fc67fe..730592c 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -35,6 +35,45 @@ static void print_vec(const char* name, const double* v, int n) {
     printf("\n");
 }
 
+static void run_once(const char* tag,
+                     const matrix_desc_t* A_desc,
+                     const double* c, const double* l, const double* u)
+{
+    printf("\n=== %s ===\n", tag);
+
+    lp_problem_t* prob = make_problem_from_matrix(
+        A_desc,  // A
+        c,       // c
+        NULL,    // objective_constant
+        NULL,    // var_lb (defaults to 0)
+        NULL,    // var_ub (defaults to +inf)
+        l,       // con_lb
+        u        // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed for %s.\n", tag);
+        return;
+    }
+
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+    lp_problem_free(prob);
+
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return;
+    }
+
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    lp_solution_free(&sol);
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -50,12 +89,57 @@ int main() {
     };
 
     // describe A using matrix_desc_t
-    matrix_desc_t A_desc;
-    A_desc.m = m;
-    A_desc.n = n;
-    A_desc.fmt = matrix_dense;
-    A_desc.zero_tolerance = 0.0;
-    A_desc.data.dense.A = &A[0][0];
+    matrix_desc_t A_dense;
+    A_dense.m = m;
+    A_dense.n = n;
+    A_dense.fmt = matrix_dense;
+    A_dense.zero_tolerance = 0.0;
+    A_dense.data.dense.A = &A[0][0];
+
+    // A as a CSR matrix
+    static int csr_row_ptr[4] = {0, 2, 3, 5};
+    static int csr_col_ind[5] = {0, 1, 1, 0, 1};
+    static double csr_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csr;
+    A_csr.m = m; A_csr.n = n;
+    A_csr.fmt = matrix_csr;
+    A_csr.zero_tolerance = 0.0;
+    A_csr.data.csr.nnz = 5;
+    A_csr.data.csr.row_ptr = csr_row_ptr;
+    A_csr.data.csr.col_ind = csr_col_ind;
+    A_csr.data.csr.vals = csr_vals;
+
+    // A as a CSC matrix
+    static int csc_col_ptr[3] = {0, 2, 5};
+    static int csc_row_ind[5] = {0, 2, 0, 1, 2};
+    static double csc_vals[5] = {1, 3, 2, 1, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csc;
+    A_csc.m = m; A_csc.n = n;
+    A_csc.fmt = matrix_csc;
+    A_csc.zero_tolerance = 0.0;
+    A_csc.data.csc.nnz = 5;
+    A_csc.data.csc.col_ptr = csc_col_ptr;
+    A_csc.data.csc.row_ind = csc_row_ind;
+    A_csc.data.csc.vals = csc_vals;
+
+    // A as a COO matrix
+    static int coo_row_ind[5] = {0, 0, 1, 2, 2};
+    static int coo_col_ind[5] = {0, 1, 1, 0, 1};
+    static double coo_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_coo;
+    A_coo.m = m; A_coo.n = n;
+    A_coo.fmt = matrix_coo;
+    A_coo.zero_tolerance = 0.0;
+    A_coo.data.coo.nnz = 5;
+    A_coo.data.coo.row_ind = coo_row_ind;
+    A_coo.data.coo.col_ind = coo_col_ind;
+    A_coo.data.coo.vals = coo_vals;
 
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
@@ -85,7 +169,7 @@ int main() {
     }
 
      lp_problem_t* prob = make_problem_from_matrix(
-        &A_desc,        // A
+        &A_dense,       // A
         c,              // c
         NULL,           // objective_constant
         NULL,           // var_lb
@@ -98,30 +182,10 @@ int main() {
         return 1;
     }
 
-    // solve problem
-    lp_solution_t sol = solve_lp_problem(prob, NULL);
-
-    // free problem
-    lp_problem_free(prob);
-
-    // check solution
-    if (!sol.x || !sol.y) {
-        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
-                term_to_str(sol.reason));
-        lp_solution_free(&sol);
-        return 2;
-    }
-
-    // print solution
-    printf("Solution:\n");
-    printf("Termination: %s\n", term_to_str(sol.reason));
-    printf("Primal obj: %.10g\n", sol.primal_obj);
-    printf("Dual   obj: %.10g\n", sol.dual_obj);
-    print_vec("x", sol.x, sol.n);
-    print_vec("y", sol.y, sol.m);
-
-    // free solution
-    lp_solution_free(&sol);
+    run_once("Test 1: Dense Matrix", &A_dense, c, l, u);
+    run_once("Test 2: CSR Matrix", &A_csr, c, l, u);
+    run_once("Test 3: CSC Matrix", &A_csc, c, l, u);
+    run_once("Test 4: COO Matrix", &A_coo, c, l, u);
 
     return 0;
 }

From cdc53bc303753e13e28848e35fa276e7ff395a19 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 19:55:30 -0400
Subject: [PATCH 17/38] Bug fixed: output interface

Unify `solve_lp_problem` API to return `cupdlpx_result_t*`
---
 cupdlpx/{interface.cu => interface.c} | 74 ++++++---------------------
 cupdlpx/interface.h                   |  5 +-
 cupdlpx/struct.h                      | 12 +----
 test/test_interface.c                 | 29 ++++++-----
 4 files changed, 34 insertions(+), 86 deletions(-)
 rename cupdlpx/{interface.cu => interface.c} (86%)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.c
similarity index 86%
rename from cupdlpx/interface.cu
rename to cupdlpx/interface.c
index 3c6cd4b..edb7441 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.c
@@ -17,7 +17,6 @@ limitations under the License.
 #include "interface.h"
 #include "solver.h"
 #include "utils.h"
-#include <cuda_runtime.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -217,16 +216,6 @@ lp_problem_t* make_problem_from_matrix(
                          &prob->constraint_matrix_num_nonzeros);
             break;
 
-        case matrix_csr:
-            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
-            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
-            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
-            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
-            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
-            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
-            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
-            break;
-
         case matrix_csc: {
             int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
             if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
@@ -254,6 +243,16 @@ lp_problem_t* make_problem_from_matrix(
             prob->constraint_matrix_values       = vals;
            break;
         }
+
+        case matrix_csr:
+            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
+            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
+            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
+            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
+            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
+            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
+            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
+            break;
             
         default:
             fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
@@ -272,17 +271,15 @@ lp_problem_t* make_problem_from_matrix(
     return prob;
 }
 
-lp_solution_t solve_lp_problem(
+cupdlpx_result_t* solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 ) {
-    // initialize output
-    lp_solution_t out = {0};
 
     // argument checks
     if (!prob) {
         fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n");
-        return out;
+        return NULL;
     }
 
     // prepare parameters: use defaults if not provided 
@@ -294,48 +291,11 @@ lp_solution_t solve_lp_problem(
     }
 
     // call optimizer
-    pdhg_solver_state_t* state = optimize(&local_params, prob);
-    if (!state) {
+    cupdlpx_result_t* res = optimize(&local_params, prob);
+    if (!res) {
         fprintf(stderr, "[interface] optimize returned NULL.\n");
-        return out;
-    }
-
-    // prepare output
-    out.n = prob->num_variables;
-    out.m = prob->num_constraints;
-    out.x = (double*)safe_malloc((size_t)out.n * sizeof(double));
-    out.y = (double*)safe_malloc((size_t)out.m * sizeof(double));
-
-    // copy solution from GPU to CPU
-    cudaError_t e1 = cudaMemcpy(out.x, state->pdhg_primal_solution,
-                                (size_t)out.n * sizeof(double), cudaMemcpyDeviceToHost);
-    cudaError_t e2 = cudaMemcpy(out.y, state->pdhg_dual_solution,
-                                (size_t)out.m * sizeof(double), cudaMemcpyDeviceToHost);
-    if (e1 != cudaSuccess || e2 != cudaSuccess) {
-        fprintf(stderr, "[interface] cudaMemcpy failed: %s / %s\n",
-                cudaGetErrorName(e1), cudaGetErrorName(e2));
-        free(out.x); out.x = NULL;
-        free(out.y); out.y = NULL;
-        pdhg_solver_state_free(state);
-        return out;
+        return NULL;
     }
-
-    out.primal_obj = state->primal_objective_value;
-    out.dual_obj = state->dual_objective_value;
-    out.reason = state->termination_reason;
-
-    pdhg_solver_state_free(state);
-    return out;
+    
+    return res;
 }
-
-// free lp solution
-void lp_solution_free(lp_solution_t* sol) {
-    if (!sol) return;
-    if (sol->x) free(sol->x);
-    if (sol->y) free(sol->y);
-    sol->n = 0;
-    sol->m = 0;
-    sol->primal_obj = 0.0;
-    sol->dual_obj   = 0.0;
-    sol->reason = TERMINATION_REASON_UNSPECIFIED;
-}
\ No newline at end of file
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index d5182c9..40fcf89 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -36,14 +36,11 @@ lp_problem_t* make_problem_from_matrix(
 );
 
 // solve the LP problem using PDHG
-lp_solution_t solve_lp_problem(
+cupdlpx_result_t* solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 );
 
-// free lp solution
-void lp_solution_free(lp_solution_t* sol);
-
 
 #ifdef __cplusplus
 } // extern "C"
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index 12daeef..300e02b 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -266,14 +266,4 @@ typedef struct {
             const double* vals; 
         } coo;
     } data;
-} matrix_desc_t;
-
-typedef struct {
-    int n;  // number of variables
-    int m;  // number of constraints
-    double* x;  // primal solution
-    double* y;  // dual solution 
-    double primal_obj; // primal objective value
-    double dual_obj; // dual objective value    
-    termination_reason_t reason; // termination reason
-} lp_solution_t;
+} matrix_desc_t;
\ No newline at end of file
diff --git a/test/test_interface.c b/test/test_interface.c
index 730592c..731e1ac 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -40,7 +40,8 @@ static void run_once(const char* tag,
                      const double* c, const double* l, const double* u)
 {
     printf("\n=== %s ===\n", tag);
-
+    
+    // build problem
     lp_problem_t* prob = make_problem_from_matrix(
         A_desc,  // A
         c,       // c
@@ -55,23 +56,23 @@ static void run_once(const char* tag,
         return;
     }
 
-    lp_solution_t sol = solve_lp_problem(prob, NULL);
+    // solve
+    cupdlpx_result_t* res = solve_lp_problem(prob, NULL);
     lp_problem_free(prob);
-
-    if (!sol.x || !sol.y) {
-        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
-                term_to_str(sol.reason));
-        lp_solution_free(&sol);
+    if (!res) {
+        fprintf(stderr, "[test] solve_lp_problem failed for %s.\n", tag);
         return;
     }
 
-    printf("Termination: %s\n", term_to_str(sol.reason));
-    printf("Primal obj: %.10g\n", sol.primal_obj);
-    printf("Dual   obj: %.10g\n", sol.dual_obj);
-    print_vec("x", sol.x, sol.n);
-    print_vec("y", sol.y, sol.m);
-
-    lp_solution_free(&sol);
+    // print results
+    printf("Termination: %s\n", term_to_str(res->termination_reason));
+    printf("Primal obj: %.10g\n", res->primal_objective_value);
+    printf("Dual   obj: %.10g\n", res->dual_objective_value);
+    print_vec("x", res->primal_solution, res->num_variables);
+    print_vec("y", res->dual_solution, res->num_constraints);
+    
+    // free
+    cupdlpx_result_free(res);
 }
 
 int main() {

From 58841ff013a4c7f3ef8e72fc9b959cc778c12277 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 20:23:51 -0400
Subject: [PATCH 18/38] Document C API and example for LP solving

Added C API section with example for solving LPs.
---
 README.md | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/README.md b/README.md
index f0cb248..fceeb42 100644
--- a/README.md
+++ b/README.md
@@ -71,6 +71,96 @@ The solver generates three text files in the specified <output_directory>. The f
 - `INSTANCE_primal_solution.txt`: The primal solution vector.
 - `INSTANCE_dual_solution.txt`: The dual solution vector.
 
+### C API Example
+
+Besides the command-line interface, cuPDLPx also provides a C API (interface.c) for directly solving LPs in memory. This is useful when integrating cuPDLPx into other C/C++ projects or when generating problems programmatically.
+
+#### Functions and Parameters
+
+The C API involves two main functions:
+
+```c
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,       // constraint matrix A
+    const double* objective_c,         // objective vector c (length n)
+    const double* objective_constant,  // scalar objective offset
+    const double* var_lb,              // variable lower bounds (length n)
+    const double* var_ub,              // variable upper bounds (length n)
+    const double* con_lb,              // constraint lower bounds (length m)
+    const double* con_ub               // constraint upper bounds (length m)
+);
+
+cupdlpx_result_t* solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params    // NULL → use default parameters
+);
+```
+
+`make_problem_from_matrix` parameters:
+- `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
+- `objective_c`: Objective vector. If NULL, defaults to all zeros.
+- `objective_constant`: Scalar constant term added to the objective value. If NULL, defaults to 0.0.
+- `var_lb`: Variable lower bounds. If NULL, defaults to all 0.0.
+- `var_ub`: Variable upper bounds. If NULL, defaults to all +INFINITY.
+- `con_lb`: Constraint lower bounds. If NULL, defaults to all -INFINITY.
+- `con_ub`: Constraint upper bounds. If NULL, defaults to all +INFINITY.
+
+`solve_lp_problem` parameters:
+- `prob`: An LP problem built with make_problem_from_matrix.
+- `params`: Solver parameters. If `NULL`, the solver will use default parameters.
+
+#### Example: Solving a Small LP
+```c
+#include "cupdlpx/interface.h"
+#include <math.h>
+#include <stdio.h>
+
+int main() {
+    int m = 3; // number of constraints
+    int n = 2; // number of variables
+
+    // Dense matrix A
+    double A[3][2] = {
+        {1.0, 2.0},
+        {0.0, 1.0},
+        {3.0, 2.0}
+    };
+
+    // Describe A
+    matrix_desc_t A_desc;
+    A_desc.m = m; A_desc.n = n;
+    A_desc.fmt = matrix_dense;
+    A_desc.zero_tolerance = 0.0;
+    A_desc.data.dense.A = &A[0][0];
+
+    // Objective coefficients
+    double c[2] = {1.0, 1.0};
+
+    // Constraint bounds: l <= A x <= u
+    double l[3] = {5.0, -INFINITY, -INFINITY};
+    double u[3] = {5.0, 2.0, 8.0};
+
+    // Build the problem
+    lp_problem_t* prob = make_problem_from_matrix(
+        &A_desc, c, NULL, NULL, NULL, l, u);
+
+    // Solve (NULL → use default parameters)
+    cupdlpx_result_t* res = solve_lp_problem(prob, NULL);
+
+    printf("Termination reason: %d\n", res->termination_reason);
+    printf("Primal objective: %.6f\n", res->primal_objective_value);
+    printf("Dual objective:   %.6f\n", res->dual_objective_value);
+    for (int j = 0; j < res->num_variables; ++j) {
+        printf("x[%d] = %.6f\n", j, res->primal_solution[j]);
+    }
+
+    lp_problem_free(prob);
+    cupdlpx_result_free(res);
+    
+    return 0;
+}
+```
+
 ## Reference
 If you use cuPDLPx or the ideas in your work, please cite the source below.
 

From 7dfe090b1609d2a0dbeba8eeae2dd5177d3e1217 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 21:18:31 -0400
Subject: [PATCH 19/38] Warning fixed: .c -> .cu for interface

---
 cupdlpx/{interface.c => interface.cu} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename cupdlpx/{interface.c => interface.cu} (100%)

diff --git a/cupdlpx/interface.c b/cupdlpx/interface.cu
similarity index 100%
rename from cupdlpx/interface.c
rename to cupdlpx/interface.cu

From 51254b4b386ecc8b99f42dd62362cb24ae6b21eb Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Fri, 19 Sep 2025 17:56:55 -0400
Subject: [PATCH 20/38] Update: .cu -> .c for interface

---
 cupdlpx/{interface.cu => interface.c} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename cupdlpx/{interface.cu => interface.c} (100%)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.c
similarity index 100%
rename from cupdlpx/interface.cu
rename to cupdlpx/interface.c

From 24d7ed1264c7e14998401c085a03f473122fb9f6 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Fri, 19 Sep 2025 22:50:18 -0400
Subject: [PATCH 21/38] Clean code: remove the repeating declaration

---
 cupdlpx/solver.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index 59af8b3..56da0bd 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -36,8 +36,6 @@ extern "C"
 
     void set_default_parameters(pdhg_parameters_t *params);
 
-    void set_default_parameters(pdhg_parameters_t *params);
-
 #ifdef __cplusplus
 }
 #endif

From b5d8233098dba113e8d55c3c6da4eaed9c8b09cf Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Sat, 20 Sep 2025 02:16:49 -0400
Subject: [PATCH 22/38] Update: `make_problem_from_matrix` ->
 `create_lp_problem`

---
 README.md             | 8 ++++----
 cupdlpx/interface.c   | 2 +-
 cupdlpx/interface.h   | 2 +-
 test/test_interface.c | 8 ++++----
 4 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/README.md b/README.md
index fceeb42..17a22d7 100644
--- a/README.md
+++ b/README.md
@@ -80,7 +80,7 @@ Besides the command-line interface, cuPDLPx also provides a C API (interface.c)
 The C API involves two main functions:
 
 ```c
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,       // constraint matrix A
     const double* objective_c,         // objective vector c (length n)
     const double* objective_constant,  // scalar objective offset
@@ -96,7 +96,7 @@ cupdlpx_result_t* solve_lp_problem(
 );
 ```
 
-`make_problem_from_matrix` parameters:
+`create_lp_problem` parameters:
 - `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
 - `objective_c`: Objective vector. If NULL, defaults to all zeros.
 - `objective_constant`: Scalar constant term added to the objective value. If NULL, defaults to 0.0.
@@ -106,7 +106,7 @@ cupdlpx_result_t* solve_lp_problem(
 - `con_ub`: Constraint upper bounds. If NULL, defaults to all +INFINITY.
 
 `solve_lp_problem` parameters:
-- `prob`: An LP problem built with make_problem_from_matrix.
+- `prob`: An LP problem built with `create_LP_problem`.
 - `params`: Solver parameters. If `NULL`, the solver will use default parameters.
 
 #### Example: Solving a Small LP
@@ -141,7 +141,7 @@ int main() {
     double u[3] = {5.0, 2.0, 8.0};
 
     // Build the problem
-    lp_problem_t* prob = make_problem_from_matrix(
+    lp_problem_t* prob = create_lp_problem(
         &A_desc, c, NULL, NULL, NULL, l, u);
 
     // Solve (NULL → use default parameters)
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
index edb7441..f457c95 100644
--- a/cupdlpx/interface.c
+++ b/cupdlpx/interface.c
@@ -192,7 +192,7 @@ static int coo_to_csr(const matrix_desc_t* desc,
 }
 
 // create an lp_problem_t from a matrix
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,
     const double* objective_c,
     const double* objective_constant,
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index 40fcf89..9263aac 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -25,7 +25,7 @@ extern "C" {
 #endif
 
 // create an lp_problem_t from a matrix descriptor
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,
     const double* objective_c,
     const double* objective_constant,
diff --git a/test/test_interface.c b/test/test_interface.c
index 731e1ac..2ec28bf 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -42,7 +42,7 @@ static void run_once(const char* tag,
     printf("\n=== %s ===\n", tag);
     
     // build problem
-    lp_problem_t* prob = make_problem_from_matrix(
+    lp_problem_t* prob = create_lp_problem(
         A_desc,  // A
         c,       // c
         NULL,    // objective_constant
@@ -52,7 +52,7 @@ static void run_once(const char* tag,
         u        // con_ub
     );
     if (!prob) {
-        fprintf(stderr, "[test] make_problem_from_matrix failed for %s.\n", tag);
+        fprintf(stderr, "[test] create_lp_problem failed for %s.\n", tag);
         return;
     }
 
@@ -169,7 +169,7 @@ int main() {
         printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
     }
 
-     lp_problem_t* prob = make_problem_from_matrix(
+     lp_problem_t* prob = create_lp_problem(
         &A_dense,       // A
         c,              // c
         NULL,           // objective_constant
@@ -179,7 +179,7 @@ int main() {
         u               // con_ub
     );
     if (!prob) {
-        fprintf(stderr, "[test] make_problem_from_matrix failed.\n");
+        fprintf(stderr, "[test] create_lp_problem failed.\n");
         return 1;
     }
 

From 748a07a7370d7ebf7737fd7de38cd85e53d2b752 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Tue, 23 Sep 2025 11:32:21 -0400
Subject: [PATCH 23/38] Update .gitignore

---
 .gitignore | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 515aa66..cd9fb58 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,4 +57,6 @@ dkms.conf
 # ignored files
 build/*
 test/*
-/.vscode
\ No newline at end of file
+/.vscode
+/.venv
+/_b

From 42c5e1c148e9b2c7958c9c97abda28069092213b Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:53:07 -0400
Subject: [PATCH 24/38] Create .gitignore

---
 .gitignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 9c13904..a84db7a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,3 @@
-
 # Prerequisites
 *.d
 
@@ -58,3 +57,4 @@ dkms.conf
 # ignored files
 build/*
 test/*
+

From b75bfb677b2a33611a1ea610612bc41c7feb916b Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 11 Sep 2025 01:55:06 -0400
Subject: [PATCH 25/38] Update .gitignore

---
 .gitignore | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index a84db7a..52defc1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,4 +57,3 @@ dkms.conf
 # ignored files
 build/*
 test/*
-

From 081548ddbafe3a1f62c51883d22efc3cbef6545e Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 02:03:02 -0400
Subject: [PATCH 26/38] Init interface

---
 .gitignore          |  1 +
 Makefile            | 73 +++++++++++++++++++++++++++++++++++++++++++--
 cupdlpx/interface.c | 18 +++++++++++
 cupdlpx/interface.h | 20 +++++++++++++
 4 files changed, 110 insertions(+), 2 deletions(-)
 create mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.h

diff --git a/.gitignore b/.gitignore
index 52defc1..654154c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -57,3 +57,4 @@ dkms.conf
 # ignored files
 build/*
 test/*
+/.vscode
diff --git a/Makefile b/Makefile
index 3fccf67..74c0c5b 100644
--- a/Makefile
+++ b/Makefile
@@ -16,6 +16,7 @@ NVCCFLAGS = -I. -I$(CUDA_HOME)/include -O3 -g -gencode arch=compute_90,code=sm_9
 # LDFLAGS for the linker
 LDFLAGS = -L$(CONDA_PREFIX)/lib -L$(CUDA_HOME)/lib64 -lcudart -lcusparse -lcublas -lz -lm
 
+# Source discovery (exclude the debug main)
 C_SOURCES = $(filter-out $(SRC_DIR)/cupdlpx.c, $(wildcard $(SRC_DIR)/*.c))
 CU_SOURCES = $(wildcard $(SRC_DIR)/*.cu)
 
@@ -25,22 +26,41 @@ OBJECTS = $(C_OBJECTS) $(CU_OBJECTS)
 
 TARGET_LIB = $(BUILD_DIR)/libcupdlpx.a
 
+# Debug executable (optional)
 DEBUG_SRC = $(SRC_DIR)/cupdlpx.c
 DEBUG_EXEC = $(BUILD_DIR)/cupdlpx
 
-.PHONY: all clean build
+# Tests auto-discovery
+TEST_DIR := ./test
+TEST_BUILD_DIR := $(BUILD_DIR)/tests
 
+TEST_CU_SOURCES := $(wildcard $(TEST_DIR)/*.cu)
+TEST_C_SOURCES := $(wildcard $(TEST_DIR)/*.c)
+
+# Each test source becomes an executable at build/tests/<basename>
+TEST_EXEC_CU := $(patsubst $(TEST_DIR)/%.cu,$(TEST_BUILD_DIR)/%,$(TEST_CU_SOURCES))
+TEST_EXEC_C := $(patsubst $(TEST_DIR)/%.c,$(TEST_BUILD_DIR)/%,$(TEST_C_SOURCES))
+
+# Phony targets
+.PHONY: all clean build tests test run-tests run-test clean-tests
+
+# Default: build the static library
 all: $(TARGET_LIB)
 
+# Archive all objects into the static library
 $(TARGET_LIB): $(OBJECTS)
 	@echo "Archiving objects into $(TARGET_LIB)..."
+	@mkdir -p $(BUILD_DIR)
 	@ar rcs $@ $^
 
+# Build the debug executable (links the library with cupdlpx.c main)
 build: $(DEBUG_EXEC)
+
 $(DEBUG_EXEC): $(DEBUG_SRC) $(TARGET_LIB)
 	@echo "Building debug executable..."
 	@$(LINKER) $(NVCCFLAGS) $(DEBUG_SRC) -o $(DEBUG_EXEC) $(TARGET_LIB) $(LDFLAGS)
 
+# Pattern rules for objects
 $(BUILD_DIR)/%.o: $(SRC_DIR)/%.c
 	@mkdir -p $(BUILD_DIR)
 	@echo "Compiling $< -> $@..."
@@ -51,6 +71,55 @@ $(BUILD_DIR)/%.o: $(SRC_DIR)/%.cu
 	@echo "Compiling $< -> $@..."
 	@$(NVCC) $(NVCCFLAGS) -c $< -o $@
 
+# Build all tests discovered under test/
+test: tests
+tests: $(TARGET_LIB) $(TEST_EXEC_CU) $(TEST_EXEC_C)
+	@echo "All tests built under $(TEST_BUILD_DIR)/"
+
+# Run all tests one by one
+run-tests: tests
+	@echo "Running all tests..."
+	@set -e; \
+	for t in $(TEST_EXEC_CU) $(TEST_EXEC_C); do \
+	  echo "=== $$t ==="; \
+	  "$$t" || exit $$?; \
+	  echo; \
+	done
+
+# Run a single test by basename: make run-test name=<basename>
+run-test: tests
+	@if [ -z "$(name)" ]; then \
+	  echo "Usage: make run-test name=<basename-of-test-file>"; exit 2; \
+	fi
+	@if [ -x "$(TEST_BUILD_DIR)/$(name)" ]; then \
+	  echo "=== $(TEST_BUILD_DIR)/$(name) ==="; \
+	  "$(TEST_BUILD_DIR)/$(name)"; \
+	else \
+	  echo "Executable not found: $(TEST_BUILD_DIR)/$(name)"; \
+	  echo "Did you 'make tests' and is there a test source named '$(TEST_DIR)/$(name).c(u)'?"; \
+	  exit 1; \
+	fi
+
+# Build rule for CUDA tests: compile and link directly with nvcc
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.cu $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building CUDA test $< -> $@..."
+	@$(LINKER) $(NVCCFLAGS) -I$(SRC_DIR) $< -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+# Build rule for C tests: compile with gcc, link with nvcc to get CUDA libs
+$(TEST_BUILD_DIR)/%: $(TEST_DIR)/%.c $(TARGET_LIB)
+	@mkdir -p $(TEST_BUILD_DIR)
+	@echo "Building C test $< -> $@..."
+	@$(CC) $(CFLAGS) -I$(SRC_DIR) -c $< -o $(TEST_BUILD_DIR)/$*.o
+	@$(LINKER) $(NVCCFLAGS) $(TEST_BUILD_DIR)/$*.o -o $@ $(TARGET_LIB) $(LDFLAGS)
+
+
+# Cleaning
+clean-tests:
+	@echo "Cleaning test executables..."
+	@rm -rf $(TEST_BUILD_DIR)
+
 clean:
 	@echo "Cleaning up..."
-	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
\ No newline at end of file
+	@rm -rf $(BUILD_DIR) $(TARGET_LIB) $(DEBUG_EXEC)
+	@rm -rf $(TEST_BUILD_DIR)
\ No newline at end of file
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
new file mode 100644
index 0000000..6313091
--- /dev/null
+++ b/cupdlpx/interface.c
@@ -0,0 +1,18 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
new file mode 100644
index 0000000..8d408f5
--- /dev/null
+++ b/cupdlpx/interface.h
@@ -0,0 +1,20 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#pragma once
+
+// Temporary placeholder for interface API
+// Functions will be declared here later
\ No newline at end of file

From ed457264f2963e9610075e96955385ecf66fc015 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 16:14:22 -0400
Subject: [PATCH 27/38] New test: LP model

---
 test/test_interface.c | 63 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 63 insertions(+)
 create mode 100644 test/test_interface.c

diff --git a/test/test_interface.c b/test/test_interface.c
new file mode 100644
index 0000000..e7d0543
--- /dev/null
+++ b/test/test_interface.c
@@ -0,0 +1,63 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "cupdlpx/interface.h"
+#include <stdio.h>
+#include <math.h> 
+
+int main() {
+    // Example: min c^T x
+    // s.t. l <= A x <= u, x >= 0
+
+    int m = 3; // number of constraints
+    int n = 2; // number of variables
+
+    // A as a dense matrix
+    double A[3][2] = {
+        {1.0, 2.0},
+        {0.0, 1.0},
+        {3.0, 2.0}
+    };
+
+    // c: objective coefficients
+    double c[2] = {1.0, 1.0};
+
+    // l&u: constraintbounds
+    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+
+    printf("Objective c = [");
+    for (int j = 0; j < n; j++) {
+        printf(" %g", c[j]);
+    }
+    printf(" ]\n");
+
+    printf("Matrix A:\n");
+    for (int i = 0; i < m; i++) {
+        printf("  [");
+        for (int j = 0; j < n; j++) {
+            printf(" %g", A[i][j]);
+        }
+        printf(" ]\n");
+    }
+
+    printf("Constraint bounds (l <= A x <= b):\n");
+    for (int i = 0; i < m; i++) {
+        printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
+    }
+
+    return 0;
+}

From a57ef9abe3e71ac673f8055245ffd33b3c9765a3 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 05:59:37 -0400
Subject: [PATCH 28/38] New feat: interface

---
 cupdlpx/cupdlpx.c     |  25 ------
 cupdlpx/interface.c   |  18 -----
 cupdlpx/interface.cu  | 184 ++++++++++++++++++++++++++++++++++++++++++
 cupdlpx/interface.h   |  34 +++++++-
 cupdlpx/solver.cu     |  26 ++++++
 cupdlpx/solver.h      |   2 +
 cupdlpx/struct.h      |  57 ++++++++++++-
 test/test_interface.c |  68 +++++++++++++++-
 8 files changed, 366 insertions(+), 48 deletions(-)
 delete mode 100644 cupdlpx/interface.c
 create mode 100644 cupdlpx/interface.cu

diff --git a/cupdlpx/cupdlpx.c b/cupdlpx/cupdlpx.c
index 023d931..de26971 100644
--- a/cupdlpx/cupdlpx.c
+++ b/cupdlpx/cupdlpx.c
@@ -130,31 +130,6 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir,
     free(file_path);
 }
 
-void set_default_parameters(pdhg_parameters_t *params)
-{
-    params->l_inf_ruiz_iterations = 10;
-    params->has_pock_chambolle_alpha = true;
-    params->pock_chambolle_alpha = 1.0;
-    params->bound_objective_rescaling = true;
-    params->verbose = false;
-    params->termination_evaluation_frequency = 200;
-    params->reflection_coefficient = 1.0;
-
-    params->termination_criteria.eps_optimal_relative = 1e-4;
-    params->termination_criteria.eps_feasible_relative = 1e-4;
-    params->termination_criteria.eps_infeasible = 1e-10;
-    params->termination_criteria.time_sec_limit = 3600.0;
-    params->termination_criteria.iteration_limit = INT32_MAX;
-
-    params->restart_params.artificial_restart_threshold = 0.36;
-    params->restart_params.sufficient_reduction_for_restart = 0.2;
-    params->restart_params.necessary_reduction_for_restart = 0.5;
-    params->restart_params.k_p = 0.99;
-    params->restart_params.k_i = 0.01;
-    params->restart_params.k_d = 0.0;
-    params->restart_params.i_smooth = 0.3;
-}
-
 void print_usage(const char *prog_name)
 {
     fprintf(stderr, "Usage: %s [OPTIONS] <mps_file> <output_dir>\n\n", prog_name);
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
deleted file mode 100644
index 6313091..0000000
--- a/cupdlpx/interface.c
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
-Copyright 2025 Haihao Lu
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "interface.h"
-
diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
new file mode 100644
index 0000000..d361169
--- /dev/null
+++ b/cupdlpx/interface.cu
@@ -0,0 +1,184 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+#include "solver.h"
+#include "utils.h"
+#include <cuda_runtime.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+// helper function to allocate and fill or copy an array
+static void fill_or_copy(double** dst, int n, const double* src, double fill_val) {
+    *dst = (double*)safe_malloc((size_t)n * sizeof(double));
+    if (src) memcpy(*dst, src, (size_t)n * sizeof(double));
+    else     for (int i = 0; i < n; ++i) (*dst)[i] = fill_val;
+}
+
+// convert dense → CSR
+static void dense_to_csr(const matrix_desc_t* desc,
+                         int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    int m = desc->m, n = desc->n;
+    double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 1e-12;
+
+    // count nnz
+    int nnz = 0;
+    for (int i = 0; i < m * n; ++i) {
+        if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    int nz = 0;
+    for (int i = 0; i < m; ++i) {
+        (*row_ptr)[i] = nz;
+        for (int j = 0; j < n; ++j) {
+            double v = desc->data.dense.A[i * n + j];
+            if (fabs(v) > tol) {
+                (*col_ind)[nz] = j;
+                (*vals)[nz] = v;
+                ++nz;
+            }
+        }
+    }
+    (*row_ptr)[m] = nz;
+    *nnz_out = nz;
+}
+
+// create an lp_problem_t from a matrix
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+) {
+    lp_problem_t* prob = (lp_problem_t*)safe_malloc(sizeof(lp_problem_t));
+
+    prob->num_variables   = A_desc->n;
+    prob->num_constraints = A_desc->m;
+
+    // handle matrix by format → convert to CSR if needed
+    switch (A_desc->fmt) {
+        case matrix_dense:
+            dense_to_csr(A_desc,
+                         &prob->constraint_matrix_row_pointers,
+                         &prob->constraint_matrix_col_indices,
+                         &prob->constraint_matrix_values,
+                         &prob->constraint_matrix_num_nonzeros);
+            break;
+
+        case matrix_csr:
+            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
+            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
+            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
+            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
+            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
+            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
+            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
+            break;
+
+        // TODO: other formats
+        default:
+            fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
+            free(prob);
+            return NULL;
+    }
+
+    // default fill values
+    prob->objective_constant = objective_constant ? *objective_constant : 0.0;
+    fill_or_copy(&prob->objective_vector, prob->num_variables, objective_c, 0.0);
+    fill_or_copy(&prob->variable_lower_bound, prob->num_variables, var_lb, 0.0);
+    fill_or_copy(&prob->variable_upper_bound, prob->num_variables, var_ub, INFINITY);
+    fill_or_copy(&prob->constraint_lower_bound, prob->num_constraints, con_lb, -INFINITY);
+    fill_or_copy(&prob->constraint_upper_bound, prob->num_constraints, con_ub, INFINITY);
+
+    return prob;
+}
+
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+) {
+    // initialize output
+    lp_solution_t out = {0};
+
+    // argument checks
+    if (!prob) {
+        fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n");
+        return out;
+    }
+
+    // prepare parameters: use defaults if not provided 
+    pdhg_parameters_t local_params;
+    if (params) {
+        local_params = *params;
+    } else {
+        set_default_parameters(&local_params);
+    }
+
+    // call optimizer
+    pdhg_solver_state_t* state = optimize(&local_params, prob);
+    if (!state) {
+        fprintf(stderr, "[interface] optimize returned NULL.\n");
+        return out;
+    }
+
+    // prepare output
+    out.n = prob->num_variables;
+    out.m = prob->num_constraints;
+    out.x = (double*)safe_malloc((size_t)out.n * sizeof(double));
+    out.y = (double*)safe_malloc((size_t)out.m * sizeof(double));
+
+    // copy solution from GPU to CPU
+    cudaError_t e1 = cudaMemcpy(out.x, state->pdhg_primal_solution,
+                                (size_t)out.n * sizeof(double), cudaMemcpyDeviceToHost);
+    cudaError_t e2 = cudaMemcpy(out.y, state->pdhg_dual_solution,
+                                (size_t)out.m * sizeof(double), cudaMemcpyDeviceToHost);
+    if (e1 != cudaSuccess || e2 != cudaSuccess) {
+        fprintf(stderr, "[interface] cudaMemcpy failed: %s / %s\n",
+                cudaGetErrorName(e1), cudaGetErrorName(e2));
+        free(out.x); out.x = NULL;
+        free(out.y); out.y = NULL;
+        pdhg_solver_state_free(state);
+        return out;
+    }
+
+    out.primal_obj = state->primal_objective_value;
+    out.dual_obj = state->dual_objective_value;
+    out.reason = state->termination_reason;
+
+    pdhg_solver_state_free(state);
+    return out;
+}
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol) {
+    if (!sol) return;
+    if (sol->x) free(sol->x);
+    if (sol->y) free(sol->y);
+    sol->n = 0;
+    sol->m = 0;
+    sol->primal_obj = 0.0;
+    sol->dual_obj   = 0.0;
+    sol->reason = TERMINATION_REASON_UNSPECIFIED;
+}
\ No newline at end of file
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index 8d408f5..d5182c9 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -16,5 +16,35 @@ limitations under the License.
 
 #pragma once
 
-// Temporary placeholder for interface API
-// Functions will be declared here later
\ No newline at end of file
+#include "struct.h"
+#include "utils.h"
+#include "io.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+// create an lp_problem_t from a matrix descriptor
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,
+    const double* objective_c,
+    const double* objective_constant,
+    const double* var_lb,
+    const double* var_ub,
+    const double* con_lb,
+    const double* con_ub
+);
+
+// solve the LP problem using PDHG
+lp_solution_t solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params
+);
+
+// free lp solution
+void lp_solution_free(lp_solution_t* sol);
+
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
\ No newline at end of file
diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu
index 9dc0b45..7d6e762 100644
--- a/cupdlpx/solver.cu
+++ b/cupdlpx/solver.cu
@@ -751,4 +751,30 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state)
     results->termination_reason = state->termination_reason;
 
     return results;
+}
+
+void set_default_parameters(pdhg_parameters_t *params)
+{
+    params->l_inf_ruiz_iterations = 10;
+    params->has_pock_chambolle_alpha = true;
+    params->pock_chambolle_alpha = 1.0;
+    params->bound_objective_rescaling = true;
+    params->verbose = false;
+    params->termination_evaluation_frequency = 200;
+    params->reflection_coefficient = 1.0;
+
+    params->termination_criteria.eps_optimal_relative = 1e-4;
+    params->termination_criteria.eps_feasible_relative = 1e-4;
+    params->termination_criteria.eps_infeasible = 1e-10;
+    params->termination_criteria.time_sec_limit = 3600.0;
+    params->termination_criteria.iteration_limit = INT32_MAX;
+
+    params->restart_params.artificial_restart_threshold = 0.36;
+    params->restart_params.sufficient_reduction_for_restart = 0.2;
+    params->restart_params.necessary_reduction_for_restart = 0.5;
+    params->restart_params.k_p = 0.99;
+    params->restart_params.k_i = 0.01;
+    params->restart_params.k_d = 0.0;
+    params->restart_params.i_smooth = 0.3;
+
 }
\ No newline at end of file
diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index a3b1750..56da0bd 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -34,6 +34,8 @@ extern "C"
 
     void cupdlpx_result_free(cupdlpx_result_t *results);
 
+    void set_default_parameters(pdhg_parameters_t *params);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index e5adc63..12daeef 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -221,4 +221,59 @@ typedef struct
 	double dual_ray_objective;
 	termination_reason_t termination_reason;
 
-} cupdlpx_result_t;
\ No newline at end of file
+} cupdlpx_result_t;
+
+// matrix formats
+typedef enum {
+    matrix_dense = 0,
+    matrix_csr   = 1,
+    matrix_csc   = 2,
+    matrix_coo   = 3
+} matrix_format_t;
+
+// matrix descriptor
+typedef struct {
+    int m; // num_constraints
+    int n; // num_variables
+    matrix_format_t fmt;
+
+    // treat abs(x) < zero_tolerance as zero
+    double zero_tolerance;
+
+    union {
+        struct { // Dense (row-major)
+            const double* A; // m*n
+        } dense;
+
+        struct { // CSR
+            int nnz;
+            const int*    row_ptr; 
+            const int*    col_ind;
+            const double* vals;
+        } csr;
+
+        struct { // CSC
+            int nnz;
+            const int*    col_ptr;
+            const int*    row_ind;
+            const double* vals;
+        } csc;
+
+        struct { // COO
+            int nnz;
+            const int*    row_ind;
+            const int*    col_ind; 
+            const double* vals; 
+        } coo;
+    } data;
+} matrix_desc_t;
+
+typedef struct {
+    int n;  // number of variables
+    int m;  // number of constraints
+    double* x;  // primal solution
+    double* y;  // dual solution 
+    double primal_obj; // primal objective value
+    double dual_obj; // dual objective value    
+    termination_reason_t reason; // termination reason
+} lp_solution_t;
diff --git a/test/test_interface.c b/test/test_interface.c
index e7d0543..6fc67fe 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -18,6 +18,23 @@ limitations under the License.
 #include <stdio.h>
 #include <math.h> 
 
+static const char* term_to_str(termination_reason_t r) {
+    switch (r) {
+        case TERMINATION_REASON_OPTIMAL:           return "OPTIMAL";
+        case TERMINATION_REASON_PRIMAL_INFEASIBLE: return "PRIMAL_INFEASIBLE";
+        case TERMINATION_REASON_DUAL_INFEASIBLE:   return "DUAL_INFEASIBLE";
+        case TERMINATION_REASON_TIME_LIMIT:        return "TIME_LIMIT";
+        case TERMINATION_REASON_ITERATION_LIMIT:   return "ITERATION_LIMIT";
+        default:                                   return "UNSPECIFIED";
+    }
+}
+
+static void print_vec(const char* name, const double* v, int n) {
+    printf("%s:", name);
+    for (int i = 0; i < n; ++i) printf(" % .6g", v[i]);
+    printf("\n");
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -32,12 +49,20 @@ int main() {
         {3.0, 2.0}
     };
 
+    // describe A using matrix_desc_t
+    matrix_desc_t A_desc;
+    A_desc.m = m;
+    A_desc.n = n;
+    A_desc.fmt = matrix_dense;
+    A_desc.zero_tolerance = 0.0;
+    A_desc.data.dense.A = &A[0][0];
+
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
 
     // l&u: constraintbounds
-    double l[3] = {5.0, -INFINITY, -1e20};  // lower bounds
-    double u[3] = {5.0, 2.0, 8.0};          // upper bounds
+    double l[3] = {5.0, -INFINITY, -INFINITY};  // lower bounds
+    double u[3] = {5.0, 2.0, 8.0}; // upper bounds
 
     printf("Objective c = [");
     for (int j = 0; j < n; j++) {
@@ -59,5 +84,44 @@ int main() {
         printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
     }
 
+     lp_problem_t* prob = make_problem_from_matrix(
+        &A_desc,        // A
+        c,              // c
+        NULL,           // objective_constant
+        NULL,           // var_lb
+        NULL,           // var_ub
+        l,              // con_lb
+        u               // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed.\n");
+        return 1;
+    }
+
+    // solve problem
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+
+    // free problem
+    lp_problem_free(prob);
+
+    // check solution
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return 2;
+    }
+
+    // print solution
+    printf("Solution:\n");
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    // free solution
+    lp_solution_free(&sol);
+
     return 0;
 }

From a8b21cae93e906e4e17cd4ebf7e84f78f73fdbbd Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 06:40:28 -0400
Subject: [PATCH 29/38] =?UTF-8?q?New=20feat:=20CSC/COO=E2=86=92CSR=20for?=
 =?UTF-8?q?=20interface?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 cupdlpx/interface.cu  | 161 +++++++++++++++++++++++++++++++++++++++++-
 test/test_interface.c | 126 +++++++++++++++++++++++++--------
 2 files changed, 254 insertions(+), 33 deletions(-)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
index d361169..3c6cd4b 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.cu
@@ -42,10 +42,12 @@ static void dense_to_csr(const matrix_desc_t* desc,
         if (fabs(desc->data.dense.A[i]) > tol) ++nnz;
     }
 
+    // allocate
     *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
     *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
     *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
 
+    // fill
     int nz = 0;
     for (int i = 0; i < m; ++i) {
         (*row_ptr)[i] = nz;
@@ -62,6 +64,134 @@ static void dense_to_csr(const matrix_desc_t* desc,
     *nnz_out = nz;
 }
 
+// convert CSC → CSR
+static int csc_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz = desc->data.csc.nnz;
+    const int *col_ptr = desc->data.csc.col_ptr;
+    const int *row_ind = desc->data.csc.row_ind; 
+    const double *v    = desc->data.csc.vals;
+
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count entries per row
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+
+    // count nnz
+    int eff_nnz = 0;
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] CSC: row index out of range\n"); return -1; }
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            ++((*row_ptr)[ri + 1]);
+            ++eff_nnz;
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // allocate
+    *col_ind = (int*)safe_malloc((size_t)eff_nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)eff_nnz * sizeof(double));
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    for (int j = 0; j < n; ++j) {
+        for (int k = col_ptr[j]; k < col_ptr[j + 1]; ++k) {
+            int ri = row_ind[k];
+            double val = v[k];
+            if (tol > 0 && fabs(val) <= tol) continue;
+            int pos = next[ri]++;
+            (*col_ind)[pos] = j;
+            (*vals)[pos]    = val;
+        }
+    }
+
+    free(next);
+    *nnz_out = eff_nnz;
+    return 0;
+}
+
+// convert COO → CSR
+static int coo_to_csr(const matrix_desc_t* desc,
+                      int** row_ptr, int** col_ind, double** vals, int* nnz_out) {
+    const int m = desc->m, n = desc->n;
+    const int nnz_in = desc->data.coo.nnz;
+    const int *r = desc->data.coo.row_ind;
+    const int *c = desc->data.coo.col_ind;
+    const double *v = desc->data.coo.vals;
+    const double tol = (desc->zero_tolerance > 0) ? desc->zero_tolerance : 0.0;
+
+    // count nnz
+    int nnz = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) ++nnz;
+    } else {
+        nnz = nnz_in;
+    }
+
+    *row_ptr = (int*)safe_malloc((size_t)(m + 1) * sizeof(int));
+    *col_ind = (int*)safe_malloc((size_t)nnz * sizeof(int));
+    *vals    = (double*)safe_malloc((size_t)nnz * sizeof(double));
+
+    // count entries per row
+    for (int i = 0; i <= m; ++i) (*row_ptr)[i] = 0;
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k)
+            if (fabs(v[k]) > tol) {
+                int ri = r[k];
+                if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+                ++((*row_ptr)[ri + 1]);
+            }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k];
+            if (ri < 0 || ri >= m) { fprintf(stderr, "[interface] COO: row index out of range\n"); return -1; }
+            ++((*row_ptr)[ri + 1]);
+        }
+    }
+
+    // exclusive scan
+    for (int i = 0; i < m; ++i) (*row_ptr)[i + 1] += (*row_ptr)[i];
+
+    // next position to fill in each row
+    int *next = (int*)safe_malloc((size_t)m * sizeof(int));
+    for (int i = 0; i < m; ++i) next[i] = (*row_ptr)[i];
+
+    // fill column indices and values
+    if (tol > 0) {
+        for (int k = 0; k < nnz_in; ++k) {
+            if (fabs(v[k]) <= tol) continue;
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    } else {
+        for (int k = 0; k < nnz_in; ++k) {
+            int ri = r[k], cj = c[k];
+            if (cj < 0 || cj >= n) { fprintf(stderr, "[interface] COO: col index out of range\n"); free(next); return -1; }
+            int pos = next[ri]++;
+            (*col_ind)[pos] = cj;
+            (*vals)[pos]    = v[k];
+        }
+    }
+
+    free(next);
+    *nnz_out = nnz;
+    return 0;
+}
+
 // create an lp_problem_t from a matrix
 lp_problem_t* make_problem_from_matrix(
     const matrix_desc_t* A_desc,
@@ -77,7 +207,7 @@ lp_problem_t* make_problem_from_matrix(
     prob->num_variables   = A_desc->n;
     prob->num_constraints = A_desc->m;
 
-    // handle matrix by format → convert to CSR if needed
+    // handle matrix by format
     switch (A_desc->fmt) {
         case matrix_dense:
             dense_to_csr(A_desc,
@@ -97,7 +227,34 @@ lp_problem_t* make_problem_from_matrix(
             memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
             break;
 
-        // TODO: other formats
+        case matrix_csc: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] CSC->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+            break;
+        }
+
+        case matrix_coo: {
+            int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
+            if (coo_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
+                fprintf(stderr, "[interface] COO->CSR failed.\n");
+                free(prob);
+                return NULL;
+            }
+            prob->constraint_matrix_num_nonzeros = nnz;
+            prob->constraint_matrix_row_pointers = row_ptr;
+            prob->constraint_matrix_col_indices  = col_ind;
+            prob->constraint_matrix_values       = vals;
+           break;
+        }
+            
         default:
             fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
             free(prob);
diff --git a/test/test_interface.c b/test/test_interface.c
index 6fc67fe..730592c 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -35,6 +35,45 @@ static void print_vec(const char* name, const double* v, int n) {
     printf("\n");
 }
 
+static void run_once(const char* tag,
+                     const matrix_desc_t* A_desc,
+                     const double* c, const double* l, const double* u)
+{
+    printf("\n=== %s ===\n", tag);
+
+    lp_problem_t* prob = make_problem_from_matrix(
+        A_desc,  // A
+        c,       // c
+        NULL,    // objective_constant
+        NULL,    // var_lb (defaults to 0)
+        NULL,    // var_ub (defaults to +inf)
+        l,       // con_lb
+        u        // con_ub
+    );
+    if (!prob) {
+        fprintf(stderr, "[test] make_problem_from_matrix failed for %s.\n", tag);
+        return;
+    }
+
+    lp_solution_t sol = solve_lp_problem(prob, NULL);
+    lp_problem_free(prob);
+
+    if (!sol.x || !sol.y) {
+        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
+                term_to_str(sol.reason));
+        lp_solution_free(&sol);
+        return;
+    }
+
+    printf("Termination: %s\n", term_to_str(sol.reason));
+    printf("Primal obj: %.10g\n", sol.primal_obj);
+    printf("Dual   obj: %.10g\n", sol.dual_obj);
+    print_vec("x", sol.x, sol.n);
+    print_vec("y", sol.y, sol.m);
+
+    lp_solution_free(&sol);
+}
+
 int main() {
     // Example: min c^T x
     // s.t. l <= A x <= u, x >= 0
@@ -50,12 +89,57 @@ int main() {
     };
 
     // describe A using matrix_desc_t
-    matrix_desc_t A_desc;
-    A_desc.m = m;
-    A_desc.n = n;
-    A_desc.fmt = matrix_dense;
-    A_desc.zero_tolerance = 0.0;
-    A_desc.data.dense.A = &A[0][0];
+    matrix_desc_t A_dense;
+    A_dense.m = m;
+    A_dense.n = n;
+    A_dense.fmt = matrix_dense;
+    A_dense.zero_tolerance = 0.0;
+    A_dense.data.dense.A = &A[0][0];
+
+    // A as a CSR matrix
+    static int csr_row_ptr[4] = {0, 2, 3, 5};
+    static int csr_col_ind[5] = {0, 1, 1, 0, 1};
+    static double csr_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csr;
+    A_csr.m = m; A_csr.n = n;
+    A_csr.fmt = matrix_csr;
+    A_csr.zero_tolerance = 0.0;
+    A_csr.data.csr.nnz = 5;
+    A_csr.data.csr.row_ptr = csr_row_ptr;
+    A_csr.data.csr.col_ind = csr_col_ind;
+    A_csr.data.csr.vals = csr_vals;
+
+    // A as a CSC matrix
+    static int csc_col_ptr[3] = {0, 2, 5};
+    static int csc_row_ind[5] = {0, 2, 0, 1, 2};
+    static double csc_vals[5] = {1, 3, 2, 1, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_csc;
+    A_csc.m = m; A_csc.n = n;
+    A_csc.fmt = matrix_csc;
+    A_csc.zero_tolerance = 0.0;
+    A_csc.data.csc.nnz = 5;
+    A_csc.data.csc.col_ptr = csc_col_ptr;
+    A_csc.data.csc.row_ind = csc_row_ind;
+    A_csc.data.csc.vals = csc_vals;
+
+    // A as a COO matrix
+    static int coo_row_ind[5] = {0, 0, 1, 2, 2};
+    static int coo_col_ind[5] = {0, 1, 1, 0, 1};
+    static double coo_vals[5] = {1, 2, 1, 3, 2};
+
+    // describe A using matrix_desc_t
+    matrix_desc_t A_coo;
+    A_coo.m = m; A_coo.n = n;
+    A_coo.fmt = matrix_coo;
+    A_coo.zero_tolerance = 0.0;
+    A_coo.data.coo.nnz = 5;
+    A_coo.data.coo.row_ind = coo_row_ind;
+    A_coo.data.coo.col_ind = coo_col_ind;
+    A_coo.data.coo.vals = coo_vals;
 
     // c: objective coefficients
     double c[2] = {1.0, 1.0};
@@ -85,7 +169,7 @@ int main() {
     }
 
      lp_problem_t* prob = make_problem_from_matrix(
-        &A_desc,        // A
+        &A_dense,       // A
         c,              // c
         NULL,           // objective_constant
         NULL,           // var_lb
@@ -98,30 +182,10 @@ int main() {
         return 1;
     }
 
-    // solve problem
-    lp_solution_t sol = solve_lp_problem(prob, NULL);
-
-    // free problem
-    lp_problem_free(prob);
-
-    // check solution
-    if (!sol.x || !sol.y) {
-        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
-                term_to_str(sol.reason));
-        lp_solution_free(&sol);
-        return 2;
-    }
-
-    // print solution
-    printf("Solution:\n");
-    printf("Termination: %s\n", term_to_str(sol.reason));
-    printf("Primal obj: %.10g\n", sol.primal_obj);
-    printf("Dual   obj: %.10g\n", sol.dual_obj);
-    print_vec("x", sol.x, sol.n);
-    print_vec("y", sol.y, sol.m);
-
-    // free solution
-    lp_solution_free(&sol);
+    run_once("Test 1: Dense Matrix", &A_dense, c, l, u);
+    run_once("Test 2: CSR Matrix", &A_csr, c, l, u);
+    run_once("Test 3: CSC Matrix", &A_csc, c, l, u);
+    run_once("Test 4: COO Matrix", &A_coo, c, l, u);
 
     return 0;
 }

From df7d32b95a977ed13a562966318622c5203c9d02 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Wed, 17 Sep 2025 02:03:02 -0400
Subject: [PATCH 30/38] Init interface

---
 cupdlpx/interface.c | 18 ++++++++++++++++++
 cupdlpx/interface.h |  4 +++-
 2 files changed, 21 insertions(+), 1 deletion(-)
 create mode 100644 cupdlpx/interface.c

diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
new file mode 100644
index 0000000..6313091
--- /dev/null
+++ b/cupdlpx/interface.c
@@ -0,0 +1,18 @@
+/*
+Copyright 2025 Haihao Lu
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include "interface.h"
+
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index d5182c9..a07ebf4 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -47,4 +47,6 @@ void lp_solution_free(lp_solution_t* sol);
 
 #ifdef __cplusplus
 } // extern "C"
-#endif
\ No newline at end of file
+#endif
+// Temporary placeholder for interface API
+// Functions will be declared here later

From 1c3e8e16879b79697e00eb64b1895cd79f6d491c Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 05:59:37 -0400
Subject: [PATCH 31/38] New feat: interface

---
 cupdlpx/interface.c   | 18 ------------------
 cupdlpx/interface.cu  |  6 +++---
 cupdlpx/interface.h   |  6 +-----
 cupdlpx/solver.cu     |  1 -
 cupdlpx/solver.h      |  2 ++
 cupdlpx/struct.h      | 12 +-----------
 test/test_interface.c |  2 +-
 7 files changed, 8 insertions(+), 39 deletions(-)
 delete mode 100644 cupdlpx/interface.c

diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
deleted file mode 100644
index 6313091..0000000
--- a/cupdlpx/interface.c
+++ /dev/null
@@ -1,18 +0,0 @@
-/*
-Copyright 2025 Haihao Lu
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-*/
-
-#include "interface.h"
-
diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.cu
index 3c6cd4b..8972235 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.cu
@@ -272,12 +272,12 @@ lp_problem_t* make_problem_from_matrix(
     return prob;
 }
 
-lp_solution_t solve_lp_problem(
+cupdlpx_result_t solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 ) {
     // initialize output
-    lp_solution_t out = {0};
+    cupdlpx_result_t out = {0};
 
     // argument checks
     if (!prob) {
@@ -329,7 +329,7 @@ lp_solution_t solve_lp_problem(
 }
 
 // free lp solution
-void lp_solution_free(lp_solution_t* sol) {
+void lp_solution_free(cupdlpx_result_t* sol) {
     if (!sol) return;
     if (sol->x) free(sol->x);
     if (sol->y) free(sol->y);
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index a07ebf4..e427e05 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -36,15 +36,11 @@ lp_problem_t* make_problem_from_matrix(
 );
 
 // solve the LP problem using PDHG
-lp_solution_t solve_lp_problem(
+cupdlpx_result_t solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 );
 
-// free lp solution
-void lp_solution_free(lp_solution_t* sol);
-
-
 #ifdef __cplusplus
 } // extern "C"
 #endif
diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu
index 7d6e762..a9de2f0 100644
--- a/cupdlpx/solver.cu
+++ b/cupdlpx/solver.cu
@@ -776,5 +776,4 @@ void set_default_parameters(pdhg_parameters_t *params)
     params->restart_params.k_i = 0.01;
     params->restart_params.k_d = 0.0;
     params->restart_params.i_smooth = 0.3;
-
 }
\ No newline at end of file
diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index 56da0bd..59af8b3 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -36,6 +36,8 @@ extern "C"
 
     void set_default_parameters(pdhg_parameters_t *params);
 
+    void set_default_parameters(pdhg_parameters_t *params);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/cupdlpx/struct.h b/cupdlpx/struct.h
index 12daeef..300e02b 100644
--- a/cupdlpx/struct.h
+++ b/cupdlpx/struct.h
@@ -266,14 +266,4 @@ typedef struct {
             const double* vals; 
         } coo;
     } data;
-} matrix_desc_t;
-
-typedef struct {
-    int n;  // number of variables
-    int m;  // number of constraints
-    double* x;  // primal solution
-    double* y;  // dual solution 
-    double primal_obj; // primal objective value
-    double dual_obj; // dual objective value    
-    termination_reason_t reason; // termination reason
-} lp_solution_t;
+} matrix_desc_t;
\ No newline at end of file
diff --git a/test/test_interface.c b/test/test_interface.c
index 730592c..49c5a4c 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -55,7 +55,7 @@ static void run_once(const char* tag,
         return;
     }
 
-    lp_solution_t sol = solve_lp_problem(prob, NULL);
+    cupdlpx_result_t sol = solve_lp_problem(prob, NULL);
     lp_problem_free(prob);
 
     if (!sol.x || !sol.y) {

From 389a7f4963e979aa97dac2cd19b18a47c41fa572 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 19:55:30 -0400
Subject: [PATCH 32/38] Bug fixed: output interface

Unify `solve_lp_problem` API to return `cupdlpx_result_t*`
---
 cupdlpx/{interface.cu => interface.c} | 75 ++++++---------------------
 cupdlpx/interface.h                   |  2 +-
 test/test_interface.c                 | 29 ++++++-----
 3 files changed, 33 insertions(+), 73 deletions(-)
 rename cupdlpx/{interface.cu => interface.c} (86%)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.c
similarity index 86%
rename from cupdlpx/interface.cu
rename to cupdlpx/interface.c
index 8972235..02f0341 100644
--- a/cupdlpx/interface.cu
+++ b/cupdlpx/interface.c
@@ -17,7 +17,6 @@ limitations under the License.
 #include "interface.h"
 #include "solver.h"
 #include "utils.h"
-#include <cuda_runtime.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -217,16 +216,6 @@ lp_problem_t* make_problem_from_matrix(
                          &prob->constraint_matrix_num_nonzeros);
             break;
 
-        case matrix_csr:
-            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
-            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
-            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
-            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
-            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
-            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
-            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
-            break;
-
         case matrix_csc: {
             int *row_ptr=NULL, *col_ind=NULL; double *vals=NULL; int nnz=0;
             if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0) {
@@ -254,6 +243,16 @@ lp_problem_t* make_problem_from_matrix(
             prob->constraint_matrix_values       = vals;
            break;
         }
+
+        case matrix_csr:
+            prob->constraint_matrix_num_nonzeros = A_desc->data.csr.nnz;
+            prob->constraint_matrix_row_pointers = (int*)safe_malloc((size_t)(A_desc->m + 1) * sizeof(int));
+            prob->constraint_matrix_col_indices = (int*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(int));
+            prob->constraint_matrix_values = (double*)safe_malloc((size_t)A_desc->data.csr.nnz * sizeof(double));
+            memcpy(prob->constraint_matrix_row_pointers, A_desc->data.csr.row_ptr, (size_t)(A_desc->m + 1) * sizeof(int));
+            memcpy(prob->constraint_matrix_col_indices, A_desc->data.csr.col_ind, (size_t)A_desc->data.csr.nnz * sizeof(int));
+            memcpy(prob->constraint_matrix_values, A_desc->data.csr.vals, (size_t)A_desc->data.csr.nnz * sizeof(double));
+            break;
             
         default:
             fprintf(stderr, "[interface] make_problem_from_matrix: unsupported matrix format %d.\n", A_desc->fmt);
@@ -272,17 +271,14 @@ lp_problem_t* make_problem_from_matrix(
     return prob;
 }
 
-cupdlpx_result_t solve_lp_problem(
+cupdlpx_result_t* solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 ) {
-    // initialize output
-    cupdlpx_result_t out = {0};
-
     // argument checks
     if (!prob) {
         fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n");
-        return out;
+        return NULL;
     }
 
     // prepare parameters: use defaults if not provided 
@@ -294,48 +290,11 @@ cupdlpx_result_t solve_lp_problem(
     }
 
     // call optimizer
-    pdhg_solver_state_t* state = optimize(&local_params, prob);
-    if (!state) {
+    cupdlpx_result_t* res = optimize(&local_params, prob);
+    if (!res) {
         fprintf(stderr, "[interface] optimize returned NULL.\n");
-        return out;
+        return NULL;
     }
-
-    // prepare output
-    out.n = prob->num_variables;
-    out.m = prob->num_constraints;
-    out.x = (double*)safe_malloc((size_t)out.n * sizeof(double));
-    out.y = (double*)safe_malloc((size_t)out.m * sizeof(double));
-
-    // copy solution from GPU to CPU
-    cudaError_t e1 = cudaMemcpy(out.x, state->pdhg_primal_solution,
-                                (size_t)out.n * sizeof(double), cudaMemcpyDeviceToHost);
-    cudaError_t e2 = cudaMemcpy(out.y, state->pdhg_dual_solution,
-                                (size_t)out.m * sizeof(double), cudaMemcpyDeviceToHost);
-    if (e1 != cudaSuccess || e2 != cudaSuccess) {
-        fprintf(stderr, "[interface] cudaMemcpy failed: %s / %s\n",
-                cudaGetErrorName(e1), cudaGetErrorName(e2));
-        free(out.x); out.x = NULL;
-        free(out.y); out.y = NULL;
-        pdhg_solver_state_free(state);
-        return out;
-    }
-
-    out.primal_obj = state->primal_objective_value;
-    out.dual_obj = state->dual_objective_value;
-    out.reason = state->termination_reason;
-
-    pdhg_solver_state_free(state);
-    return out;
-}
-
-// free lp solution
-void lp_solution_free(cupdlpx_result_t* sol) {
-    if (!sol) return;
-    if (sol->x) free(sol->x);
-    if (sol->y) free(sol->y);
-    sol->n = 0;
-    sol->m = 0;
-    sol->primal_obj = 0.0;
-    sol->dual_obj   = 0.0;
-    sol->reason = TERMINATION_REASON_UNSPECIFIED;
+    
+    return res;
 }
\ No newline at end of file
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index e427e05..ae4ce04 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -36,7 +36,7 @@ lp_problem_t* make_problem_from_matrix(
 );
 
 // solve the LP problem using PDHG
-cupdlpx_result_t solve_lp_problem(
+cupdlpx_result_t* solve_lp_problem(
     const lp_problem_t* prob,
     const pdhg_parameters_t* params
 );
diff --git a/test/test_interface.c b/test/test_interface.c
index 49c5a4c..731e1ac 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -40,7 +40,8 @@ static void run_once(const char* tag,
                      const double* c, const double* l, const double* u)
 {
     printf("\n=== %s ===\n", tag);
-
+    
+    // build problem
     lp_problem_t* prob = make_problem_from_matrix(
         A_desc,  // A
         c,       // c
@@ -55,23 +56,23 @@ static void run_once(const char* tag,
         return;
     }
 
-    cupdlpx_result_t sol = solve_lp_problem(prob, NULL);
+    // solve
+    cupdlpx_result_t* res = solve_lp_problem(prob, NULL);
     lp_problem_free(prob);
-
-    if (!sol.x || !sol.y) {
-        fprintf(stderr, "[test] solve_lp_problem failed (x/y null). reason=%s\n",
-                term_to_str(sol.reason));
-        lp_solution_free(&sol);
+    if (!res) {
+        fprintf(stderr, "[test] solve_lp_problem failed for %s.\n", tag);
         return;
     }
 
-    printf("Termination: %s\n", term_to_str(sol.reason));
-    printf("Primal obj: %.10g\n", sol.primal_obj);
-    printf("Dual   obj: %.10g\n", sol.dual_obj);
-    print_vec("x", sol.x, sol.n);
-    print_vec("y", sol.y, sol.m);
-
-    lp_solution_free(&sol);
+    // print results
+    printf("Termination: %s\n", term_to_str(res->termination_reason));
+    printf("Primal obj: %.10g\n", res->primal_objective_value);
+    printf("Dual   obj: %.10g\n", res->dual_objective_value);
+    print_vec("x", res->primal_solution, res->num_variables);
+    print_vec("y", res->dual_solution, res->num_constraints);
+    
+    // free
+    cupdlpx_result_free(res);
 }
 
 int main() {

From 7b5a3f35e6f1f04035faae966c740e2d8789e0ad Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 21:18:31 -0400
Subject: [PATCH 33/38] Warning fixed: .c -> .cu for interface

---
 cupdlpx/{interface.c => interface.cu} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename cupdlpx/{interface.c => interface.cu} (100%)

diff --git a/cupdlpx/interface.c b/cupdlpx/interface.cu
similarity index 100%
rename from cupdlpx/interface.c
rename to cupdlpx/interface.cu

From 40614ae6969fe0f5dcc5545edf0aaca936dfa821 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Thu, 18 Sep 2025 20:23:51 -0400
Subject: [PATCH 34/38] Document C API and example for LP solving

Added C API section with example for solving LPs.
---
 README.md | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)

diff --git a/README.md b/README.md
index f0cb248..fceeb42 100644
--- a/README.md
+++ b/README.md
@@ -71,6 +71,96 @@ The solver generates three text files in the specified <output_directory>. The f
 - `INSTANCE_primal_solution.txt`: The primal solution vector.
 - `INSTANCE_dual_solution.txt`: The dual solution vector.
 
+### C API Example
+
+Besides the command-line interface, cuPDLPx also provides a C API (interface.c) for directly solving LPs in memory. This is useful when integrating cuPDLPx into other C/C++ projects or when generating problems programmatically.
+
+#### Functions and Parameters
+
+The C API involves two main functions:
+
+```c
+lp_problem_t* make_problem_from_matrix(
+    const matrix_desc_t* A_desc,       // constraint matrix A
+    const double* objective_c,         // objective vector c (length n)
+    const double* objective_constant,  // scalar objective offset
+    const double* var_lb,              // variable lower bounds (length n)
+    const double* var_ub,              // variable upper bounds (length n)
+    const double* con_lb,              // constraint lower bounds (length m)
+    const double* con_ub               // constraint upper bounds (length m)
+);
+
+cupdlpx_result_t* solve_lp_problem(
+    const lp_problem_t* prob,
+    const pdhg_parameters_t* params    // NULL → use default parameters
+);
+```
+
+`make_problem_from_matrix` parameters:
+- `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
+- `objective_c`: Objective vector. If NULL, defaults to all zeros.
+- `objective_constant`: Scalar constant term added to the objective value. If NULL, defaults to 0.0.
+- `var_lb`: Variable lower bounds. If NULL, defaults to all 0.0.
+- `var_ub`: Variable upper bounds. If NULL, defaults to all +INFINITY.
+- `con_lb`: Constraint lower bounds. If NULL, defaults to all -INFINITY.
+- `con_ub`: Constraint upper bounds. If NULL, defaults to all +INFINITY.
+
+`solve_lp_problem` parameters:
+- `prob`: An LP problem built with make_problem_from_matrix.
+- `params`: Solver parameters. If `NULL`, the solver will use default parameters.
+
+#### Example: Solving a Small LP
+```c
+#include "cupdlpx/interface.h"
+#include <math.h>
+#include <stdio.h>
+
+int main() {
+    int m = 3; // number of constraints
+    int n = 2; // number of variables
+
+    // Dense matrix A
+    double A[3][2] = {
+        {1.0, 2.0},
+        {0.0, 1.0},
+        {3.0, 2.0}
+    };
+
+    // Describe A
+    matrix_desc_t A_desc;
+    A_desc.m = m; A_desc.n = n;
+    A_desc.fmt = matrix_dense;
+    A_desc.zero_tolerance = 0.0;
+    A_desc.data.dense.A = &A[0][0];
+
+    // Objective coefficients
+    double c[2] = {1.0, 1.0};
+
+    // Constraint bounds: l <= A x <= u
+    double l[3] = {5.0, -INFINITY, -INFINITY};
+    double u[3] = {5.0, 2.0, 8.0};
+
+    // Build the problem
+    lp_problem_t* prob = make_problem_from_matrix(
+        &A_desc, c, NULL, NULL, NULL, l, u);
+
+    // Solve (NULL → use default parameters)
+    cupdlpx_result_t* res = solve_lp_problem(prob, NULL);
+
+    printf("Termination reason: %d\n", res->termination_reason);
+    printf("Primal objective: %.6f\n", res->primal_objective_value);
+    printf("Dual objective:   %.6f\n", res->dual_objective_value);
+    for (int j = 0; j < res->num_variables; ++j) {
+        printf("x[%d] = %.6f\n", j, res->primal_solution[j]);
+    }
+
+    lp_problem_free(prob);
+    cupdlpx_result_free(res);
+    
+    return 0;
+}
+```
+
 ## Reference
 If you use cuPDLPx or the ideas in your work, please cite the source below.
 

From c97be29178bc871f3733a204a1afcc183719d05d Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Fri, 19 Sep 2025 17:56:55 -0400
Subject: [PATCH 35/38] Update: .cu -> .c for interface

---
 cupdlpx/{interface.cu => interface.c} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename cupdlpx/{interface.cu => interface.c} (100%)

diff --git a/cupdlpx/interface.cu b/cupdlpx/interface.c
similarity index 100%
rename from cupdlpx/interface.cu
rename to cupdlpx/interface.c

From 0dc5fe7daa154e23bc159dddcfd703508569d22e Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Fri, 19 Sep 2025 22:50:18 -0400
Subject: [PATCH 36/38] Clean code: remove the repeating declaration

---
 cupdlpx/solver.h | 2 --
 1 file changed, 2 deletions(-)

diff --git a/cupdlpx/solver.h b/cupdlpx/solver.h
index 59af8b3..56da0bd 100644
--- a/cupdlpx/solver.h
+++ b/cupdlpx/solver.h
@@ -36,8 +36,6 @@ extern "C"
 
     void set_default_parameters(pdhg_parameters_t *params);
 
-    void set_default_parameters(pdhg_parameters_t *params);
-
 #ifdef __cplusplus
 }
 #endif

From 4f7e19c6d3c75ee668b69ad308abecd544605e5d Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Sat, 20 Sep 2025 02:16:49 -0400
Subject: [PATCH 37/38] Update: `make_problem_from_matrix` ->
 `create_lp_problem`

---
 README.md             | 8 ++++----
 cupdlpx/interface.c   | 2 +-
 cupdlpx/interface.h   | 2 +-
 test/test_interface.c | 8 ++++----
 4 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/README.md b/README.md
index fceeb42..17a22d7 100644
--- a/README.md
+++ b/README.md
@@ -80,7 +80,7 @@ Besides the command-line interface, cuPDLPx also provides a C API (interface.c)
 The C API involves two main functions:
 
 ```c
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,       // constraint matrix A
     const double* objective_c,         // objective vector c (length n)
     const double* objective_constant,  // scalar objective offset
@@ -96,7 +96,7 @@ cupdlpx_result_t* solve_lp_problem(
 );
 ```
 
-`make_problem_from_matrix` parameters:
+`create_lp_problem` parameters:
 - `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
 - `objective_c`: Objective vector. If NULL, defaults to all zeros.
 - `objective_constant`: Scalar constant term added to the objective value. If NULL, defaults to 0.0.
@@ -106,7 +106,7 @@ cupdlpx_result_t* solve_lp_problem(
 - `con_ub`: Constraint upper bounds. If NULL, defaults to all +INFINITY.
 
 `solve_lp_problem` parameters:
-- `prob`: An LP problem built with make_problem_from_matrix.
+- `prob`: An LP problem built with `create_LP_problem`.
 - `params`: Solver parameters. If `NULL`, the solver will use default parameters.
 
 #### Example: Solving a Small LP
@@ -141,7 +141,7 @@ int main() {
     double u[3] = {5.0, 2.0, 8.0};
 
     // Build the problem
-    lp_problem_t* prob = make_problem_from_matrix(
+    lp_problem_t* prob = create_lp_problem(
         &A_desc, c, NULL, NULL, NULL, l, u);
 
     // Solve (NULL → use default parameters)
diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c
index 02f0341..c851ebb 100644
--- a/cupdlpx/interface.c
+++ b/cupdlpx/interface.c
@@ -192,7 +192,7 @@ static int coo_to_csr(const matrix_desc_t* desc,
 }
 
 // create an lp_problem_t from a matrix
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,
     const double* objective_c,
     const double* objective_constant,
diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h
index ae4ce04..467339a 100644
--- a/cupdlpx/interface.h
+++ b/cupdlpx/interface.h
@@ -25,7 +25,7 @@ extern "C" {
 #endif
 
 // create an lp_problem_t from a matrix descriptor
-lp_problem_t* make_problem_from_matrix(
+lp_problem_t* create_lp_problem(
     const matrix_desc_t* A_desc,
     const double* objective_c,
     const double* objective_constant,
diff --git a/test/test_interface.c b/test/test_interface.c
index 731e1ac..2ec28bf 100644
--- a/test/test_interface.c
+++ b/test/test_interface.c
@@ -42,7 +42,7 @@ static void run_once(const char* tag,
     printf("\n=== %s ===\n", tag);
     
     // build problem
-    lp_problem_t* prob = make_problem_from_matrix(
+    lp_problem_t* prob = create_lp_problem(
         A_desc,  // A
         c,       // c
         NULL,    // objective_constant
@@ -52,7 +52,7 @@ static void run_once(const char* tag,
         u        // con_ub
     );
     if (!prob) {
-        fprintf(stderr, "[test] make_problem_from_matrix failed for %s.\n", tag);
+        fprintf(stderr, "[test] create_lp_problem failed for %s.\n", tag);
         return;
     }
 
@@ -169,7 +169,7 @@ int main() {
         printf("  row %d: l = %g, u = %g\n", i, l[i], u[i]);
     }
 
-     lp_problem_t* prob = make_problem_from_matrix(
+     lp_problem_t* prob = create_lp_problem(
         &A_dense,       // A
         c,              // c
         NULL,           // objective_constant
@@ -179,7 +179,7 @@ int main() {
         u               // con_ub
     );
     if (!prob) {
-        fprintf(stderr, "[test] make_problem_from_matrix failed.\n");
+        fprintf(stderr, "[test] create_lp_problem failed.\n");
         return 1;
     }
 

From 3eb933b48555250407c6319c5a13242df0d8d6d6 Mon Sep 17 00:00:00 2001
From: Bo Tang <botang@mie.utoronto.ca>
Date: Tue, 23 Sep 2025 11:32:21 -0400
Subject: [PATCH 38/38] Update .gitignore

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index 654154c..cd9fb58 100644
--- a/.gitignore
+++ b/.gitignore
@@ -58,3 +58,5 @@ dkms.conf
 build/*
 test/*
 /.vscode
+/.venv
+/_b