diff --git a/.gitignore b/.gitignore index 9c13904..cd9fb58 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,3 @@ - # Prerequisites *.d @@ -58,3 +57,6 @@ dkms.conf # ignored files build/* test/* +/.vscode +/.venv +/_b 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/ +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= +run-test: tests + @if [ -z "$(name)" ]; then \ + echo "Usage: make run-test name="; 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/README.md b/README.md index f0cb248..17a22d7 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,96 @@ The solver generates three text files in the specified . 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* 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 + 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 +); +``` + +`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. +- `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 `create_LP_problem`. +- `params`: Solver parameters. If `NULL`, the solver will use default parameters. + +#### Example: Solving a Small LP +```c +#include "cupdlpx/interface.h" +#include +#include + +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 = create_lp_problem( + &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. 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] \n\n", prog_name); diff --git a/cupdlpx/interface.c b/cupdlpx/interface.c new file mode 100644 index 0000000..c851ebb --- /dev/null +++ b/cupdlpx/interface.c @@ -0,0 +1,300 @@ +/* +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 +#include +#include +#include + +// 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; + } + + // 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; + 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; +} + +// 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* create_lp_problem( + 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 + 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_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; + } + + 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); + 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; +} + +cupdlpx_result_t* solve_lp_problem( + const lp_problem_t* prob, + const pdhg_parameters_t* params +) { + // argument checks + if (!prob) { + fprintf(stderr, "[interface] solve_lp_problem: invalid arguments.\n"); + return NULL; + } + + // 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 + cupdlpx_result_t* res = optimize(&local_params, prob); + if (!res) { + fprintf(stderr, "[interface] optimize returned NULL.\n"); + return NULL; + } + + return res; +} \ No newline at end of file diff --git a/cupdlpx/interface.h b/cupdlpx/interface.h new file mode 100644 index 0000000..1846ae0 --- /dev/null +++ b/cupdlpx/interface.h @@ -0,0 +1,46 @@ +/* +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 + +#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* create_lp_problem( + 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 +cupdlpx_result_t* solve_lp_problem( + const lp_problem_t* prob, + const pdhg_parameters_t* params +); + +#ifdef __cplusplus +} // extern "C" +#endif \ No newline at end of file diff --git a/cupdlpx/solver.cu b/cupdlpx/solver.cu index 9dc0b45..a9de2f0 100644 --- a/cupdlpx/solver.cu +++ b/cupdlpx/solver.cu @@ -751,4 +751,29 @@ 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..300e02b 100644 --- a/cupdlpx/struct.h +++ b/cupdlpx/struct.h @@ -221,4 +221,49 @@ 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; \ No newline at end of file diff --git a/test/test_interface.c b/test/test_interface.c new file mode 100644 index 0000000..2ec28bf --- /dev/null +++ b/test/test_interface.c @@ -0,0 +1,192 @@ +/* +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 +#include + +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"); +} + +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); + + // build problem + lp_problem_t* prob = create_lp_problem( + 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] create_lp_problem failed for %s.\n", tag); + return; + } + + // solve + cupdlpx_result_t* res = solve_lp_problem(prob, NULL); + lp_problem_free(prob); + if (!res) { + fprintf(stderr, "[test] solve_lp_problem failed for %s.\n", tag); + return; + } + + // 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() { + // 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} + }; + + // describe A using matrix_desc_t + 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}; + + // l&u: constraintbounds + 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++) { + 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]); + } + + lp_problem_t* prob = create_lp_problem( + &A_dense, // A + c, // c + NULL, // objective_constant + NULL, // var_lb + NULL, // var_ub + l, // con_lb + u // con_ub + ); + if (!prob) { + fprintf(stderr, "[test] create_lp_problem failed.\n"); + return 1; + } + + 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; +}