diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 128f361..60775a5 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: false matrix: - os: [ubuntu-latest, macos-latest, windows-2019] + os: [ubuntu-latest, macos-latest, windows-2022] runs-on: ${{ matrix.os }} @@ -24,18 +24,15 @@ jobs: submodules: 'recursive' - name: Install MATLAB - uses: matlab-actions/setup-matlab@v1 - with: - # Use R2023a due to a licensing bug with R2023b and the GitHub actions shared runners - release: R2023a + uses: matlab-actions/setup-matlab@v2 - name: Build OSQP interface - uses: matlab-actions/run-command@v1 + uses: matlab-actions/run-command@v2 with: - command: make_osqp + command: osqp.build('osqp_mex') - name: Run tests - uses: matlab-actions/run-tests@v1 + uses: matlab-actions/run-tests@v2 with: source-folder: ./ select-by-folder: unittests diff --git a/.gitignore b/.gitignore index 225d671..c3d34c7 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ # ------------------------------------------------------------------- # Out folder out/ +c_sources/build/ # Prerequisites *.d diff --git a/.gitmodules b/.gitmodules index d8ab27a..a1ed43b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,3 @@ [submodule "osqp_sources"] - path = osqp_sources + path = c_sources/osqp_sources url = https://github.com/osqp/osqp diff --git a/@osqp/build.m b/@osqp/build.m new file mode 100644 index 0000000..2d1b484 --- /dev/null +++ b/@osqp/build.m @@ -0,0 +1,151 @@ +% SPDX-License-Identifier: Apache-2.0 + +function build(varargin) +% Matlab MEX makefile for OSQP. +% +% MAKE_OSQP(VARARGIN) is a make file for OSQP solver. It +% builds OSQP and its components from source. +% +% WHAT is the last element of VARARGIN and cell array of strings, +% with the following options: +% +% {}, '' (empty string) or 'all': build all components and link. +% +% 'osqp_mex': builds the OSQP mex interface and the OSQP library +% +% Additional commands: +% +% 'clean': Delete all compiled files +% 'purge': Delete all compiled files and copied code generation files + + if( nargin == 0 ) + what = {'all'}; + verbose = false; + elseif ( nargin == 1 && ismember('-verbose', varargin) ) + what = {'all'}; + verbose = true; + else + what = varargin{nargin}; + if(isempty(strfind(what, 'all')) && ... + isempty(strfind(what, 'osqp_mex')) && ... + isempty(strfind(what, 'clean')) && ... + isempty(strfind(what, 'purge'))) + fprintf('No rule to make target "%s", exiting.\n', what); + end + + verbose = ismember('-verbose', varargin); + end + + %% Determine where the various files are all located + % Various parts of the build system + [osqp_classpath,~,~] = fileparts( mfilename( 'fullpath' ) ); + osqp_mex_src_dir = fullfile( osqp_classpath, '..', 'c_sources' ); + osqp_mex_build_dir = fullfile( osqp_mex_src_dir, 'build' ); + osqp_cg_dest_dir = fullfile( osqp_classpath, '..', 'codegen', 'sources' ); + + % Determine where CMake should look for MATLAB + Matlab_ROOT = strrep( matlabroot, '\', '/' ); + + %% Try to unlock any pre-existing version of osqp_mex + % this prevents compile errors if a user builds, runs osqp + % and then tries to recompile + if(mislocked('osqp_mex')) + munlock('osqp_mex'); + end + + %% Configure, build and install the OSQP mex interface + if( any(strcmpi(what,'osqp_mex')) || any(strcmpi(what,'all')) ) + fprintf('Compiling OSQP solver mex interface...\n'); + + % Create build for the mex file and go inside + if exist( osqp_mex_build_dir, 'dir' ) + rmdir( osqp_mex_build_dir, 's' ); + end + mkdir( osqp_mex_build_dir ); + % cd( osqp_mex_build_dir ); + + % Extend path for CMake mac (via Homebrew) + PATH = getenv('PATH'); + if( (ismac) && (isempty(strfind(PATH, '/usr/local/bin'))) ) + setenv('PATH', [PATH ':/usr/local/bin']); + end + + + %% Configure CMake for the mex interface + fprintf(' Configuring...' ) + temp_dir = fullfile(tempdir, ['osqp_matlab_build_' datestr(now, 'yyyymmdd_HHMMSS_FFF')]); + [status, output] = system( sprintf( 'cmake -B %s -S %s -DCMAKE_INSTALL_PREFIX=\"%s\" -DCMAKE_BUILD_TYPE=RelWithDebInfo -DMatlab_ROOT_DIR=\"%s\" -DOSQP_CODEGEN_INSTALL_DIR=\"%s\"', osqp_mex_build_dir, osqp_mex_src_dir, temp_dir, Matlab_ROOT, osqp_cg_dest_dir), 'LD_LIBRARY_PATH', '' ); + if( status ) + fprintf( '\n' ); + disp( output ); + error( 'Error configuring CMake environment' ); + elseif( verbose ) + fprintf( '\n' ); + disp( output ); + else + fprintf( '\t\t\t\t\t[done]\n' ); + end + + %% Build the mex interface + fprintf( ' Building...') + [status, output] = system( sprintf( 'cmake --build %s --config Release --target install', osqp_mex_build_dir ), 'LD_LIBRARY_PATH', '' ); + if( status ) + fprintf( '\n' ); + disp( output ); + error( 'Error compiling OSQP mex interface' ); + elseif( verbose ) + fprintf( '\n' ); + disp( output ); + else + fprintf( '\t\t\t\t\t\t[done]\n' ); + end + + + %% Install various files + fprintf( ' Installing...' ) + + % Copy mex file to root directory for use + if( ispc ) + [err, errmsg, ~] = copyfile( [osqp_mex_build_dir, filesep, 'Release', filesep, 'osqp_mex.mex*'], [osqp_classpath, filesep, 'private'] ); + else + [err, errmsg, ~] = copyfile( [osqp_mex_build_dir, filesep, 'osqp_mex.mex*'], [osqp_classpath, filesep, 'private'] ); + end + if( ~err ) + fprintf( '\n' ) + disp( errmsg ) + error( ' Error copying mex file' ) + end + + fprintf( '\t\t\t\t\t\t[done]\n' ); + end + + %% Clean and purge + if( any(strcmpi(what,'clean')) || any(strcmpi(what,'purge')) ) + fprintf('Cleaning OSQP mex files and build directory...'); + + % Delete mex file + mexfiles = dir(['*.', mexext]); + for i = 1 : length(mexfiles) + delete(mexfiles(i).name); + end + + % Delete OSQP build directory + if exist(osqp_mex_build_dir, 'dir') + rmdir(osqp_mex_build_dir, 's'); + end + + fprintf('\t\t[done]\n'); + + %% Purge only + if( any(strcmpi(what,'purge')) ) + fprintf('Cleaning OSQP codegen directories...'); + + % Delete codegen files + if exist(osqp_cg_dest_dir, 'dir') + rmdir(osqp_cg_dest_dir, 's'); + end + + fprintf('\t\t\t[done]\n'); + end + end +end diff --git a/@osqp/codegen.m b/@osqp/codegen.m new file mode 100644 index 0000000..ed471d9 --- /dev/null +++ b/@osqp/codegen.m @@ -0,0 +1,70 @@ +% SPDX-License-Identifier: Apache-2.0 + +%% +function codegen(this, out, varargin) + % CODEGEN generate C code for the parametric problem + % + % codegen(target_dir,options) + + % Parse input arguments + p = inputParser; + defaultPrefix = 'prob1_'; % Prefix for filenames and C variables; useful if generating multiple problems + defaultForceRewrite = true; % Force rewrite if output folder exists? + defaultParameters = 'vectors'; % What do we wish to update in the generated code? + % One of 'vectors' (allowing update of q/l/u through prob.update_data_vec) + % or 'matrices' (allowing update of P/A/q/l/u + % through prob.update_data_vec or prob.update_data_mat) + defaultUseFloat = false; % Use single precision in generated code? + defaultPrintingEnable = false; % Enable solver printing? + defaultProfilingEnable = false; % Enable solver profiling? + defaultInterruptEnable = false; % Enable user interrupt (Ctrl-C)? + defaultEnableDerivatives = false; % Enable derivatives? + + addRequired(p, 'out', @isstr); + addOptional(p, 'prefix', defaultPrefix, @isstr); + addParameter(p, 'force_rewrite', defaultForceRewrite, @isboolean); + addParameter(p, 'parameters', defaultParameters, @isstr); + addParameter(p, 'float_type', defaultUseFloat, @isboolean); + addParameter(p, 'printing_enable', defaultPrintingEnable, @isboolean); + addParameter(p, 'profiling_enable', defaultProfilingEnable, @isboolean); + addParameter(p, 'interrupt_enable', defaultInterruptEnable, @isboolean); + addParameter(p, 'derivatives_enable', defaultEnableDerivatives, @isboolean); + + parse(p, out, varargin{:}); + + % Set internal variables + if strcmp(p.Results.parameters, 'vectors') + embedded = 1; + else + embedded = 2; + end + + + % Check whether the specified directory already exists + if exist(out, 'dir') + while(1) + prompt = sprintf('Directory "%s" already exists. Do you want to replace it? y/n [y]: ', out); + str = input(prompt, 's'); + + if any(strcmpi(str, {'','y'})) + rmdir(out, 's'); + break; + elseif strcmpi(str, 'n') + return; + end + end + end + + % Import OSQP path + [osqp_path,~,~] = fileparts(which('osqp.m')); + + % Path to codegen source + cg_dir = fullfile(osqp_path, '..', 'codegen', 'sources'); + copyfile(cg_dir, out); + + % Update codegen defines + update_codegen_defines(this, 'embedded_mode', embedded, 'float_type', p.Results.float_type, 'printing_enable', p.Results.printing_enable, 'profiling_enable', p.Results.profiling_enable, 'interrupt_enable', p.Results.interrupt_enable, 'derivatives_enable', p.Results.derivatives_enable); + % Call codegen + osqp_mex('codegen', this.objectHandle, out, p.Results.prefix); + +end \ No newline at end of file diff --git a/@osqp/osqp.m b/@osqp/osqp.m new file mode 100644 index 0000000..7d64c82 --- /dev/null +++ b/@osqp/osqp.m @@ -0,0 +1,93 @@ +% SPDX-License-Identifier: Apache-2.0 + +classdef osqp < handle + % osqp interface class for OSQP solver + % This class provides a complete interface to the C implementation + % of the OSQP solver. + % + % osqp Properties: + % objectHandle - pointer to the C structure of OSQP solver + % + % osqp Methods: + % + % setup - configure solver with problem data + % solve - solve the QP + % update - modify problem vectors + % warm_start - set warm starting variables x and y + % + % default_settings - create default settings structure + % current_settings - get the current solver settings structure + % update_settings - update the current solver settings structure + % + % get_dimensions - get the number of variables and constraints + % version - return OSQP version + % constant - return a OSQP internal constant + % + % update_codegen_defines - update the current codegen defines + % codegen - generate embeddable C code for the problem + + + properties(SetAccess = private, Hidden = true) + objectHandle % Handle to underlying C instance + end + + methods(Static) + output = build(varargin) + + %% + function out = default_settings() + % DEFAULT_SETTINGS get the default solver settings structure + out = osqp_mex('default_settings'); + + % Convert linsys solver to string + out.linsys_solver = linsys_solver_to_string(out.linsys_solver); + end + + %% + function out = constant(constant_name) + % CONSTANT Return solver constant + % C = CONSTANT(CONSTANT_NAME) return constant called CONSTANT_NAME + out = osqp_mex('constant', constant_name); + end + + %% + function out = version() + % Return OSQP version + out = osqp_mex('version'); + end + end + + methods(Access = private) + currentSettings = validate_settings(this, isInitialization, varargin) + end + + methods + %% Constructor - Create a new solver instance + function this = osqp(varargin) + % Construct OSQP solver class + this.objectHandle = osqp_mex('new', varargin{:}); + end + + %% Destructor - destroy the solver instance + function delete(this) + % Destroy OSQP solver class + osqp_mex('delete', this.objectHandle); + end + + %% + function out = current_settings(this) + % CURRENT_SETTINGS get the current solver settings structure + out = osqp_mex('current_settings', this.objectHandle); + + % Convert linsys solver to string + out.linsys_solver = linsys_solver_to_string(out.linsys_solver); + end + + %% + function [n,m] = get_dimensions(this) + % GET_DIMENSIONS get the number of variables and constraints + + [n,m] = osqp_mex('get_dimensions', this.objectHandle); + end + end +end \ No newline at end of file diff --git a/@osqp/private/linsys_solver_to_string.m b/@osqp/private/linsys_solver_to_string.m new file mode 100644 index 0000000..d1ea38e --- /dev/null +++ b/@osqp/private/linsys_solver_to_string.m @@ -0,0 +1,15 @@ +% SPDX-License-Identifier: Apache-2.0 + +% Convert linear systme solver integer to string +function [linsys_solver_string] = linsys_solver_to_string(linsys_solver) + switch linsys_solver + case osqp.constant('OSQP_UNKNOWN_SOLVER') + linsys_solver_string = 'unknown solver'; + case osqp.constant('OSQP_DIRECT_SOLVER') + linsys_solver_string = 'direct solver'; + case osqp.constant('OSQP_INDIRECT_SOLVER') + linsys_solver_string = 'indirect solver'; + otherwise + error('Unrecognized linear system solver.'); + end +end diff --git a/@osqp/private/string_to_linsys_solver.m b/@osqp/private/string_to_linsys_solver.m new file mode 100644 index 0000000..3c4f241 --- /dev/null +++ b/@osqp/private/string_to_linsys_solver.m @@ -0,0 +1,20 @@ +% SPDX-License-Identifier: Apache-2.0 + +function [linsys_solver] = string_to_linsys_solver(linsys_solver_string) + linsys_solver_string = lower(linsys_solver_string); + switch linsys_solver_string + case 'unknown solver' + linsys_solver = osqp.constant('OSQP_UNKNOWN_SOLVER'); + case 'direct solver' + linsys_solver = osqp.constant('OSQP_DIRECT_SOLVER'); + case 'indirect solver' + linsys_solver = osqp.constant('OSQP_INDIRECT_SOLVER'); + % Default solver: QDLDL + case '' + linsys_solver = osqp.constant('OSQP_DIRECT_SOLVER'); + otherwise + warning('Linear system solver not recognized. Using default solver OSQP_DIRECT_SOLVER.') + linsys_solver = osqp.constant('OSQP_DIRECT_SOLVER'); + end +end + \ No newline at end of file diff --git a/@osqp/setup.m b/@osqp/setup.m new file mode 100644 index 0000000..f70b662 --- /dev/null +++ b/@osqp/setup.m @@ -0,0 +1,103 @@ +% SPDX-License-Identifier: Apache-2.0 + +%% +function varargout = setup(this, varargin) + % SETUP configure solver with problem data + % + % setup(P,q,A,l,u,options) + + nargin = length(varargin); + + %dimension checks on user data. Mex function does not + %perform any checks on inputs, so check everything here + assert(nargin >= 5, 'incorrect number of inputs'); + [P,q,A,l,u] = deal(varargin{1:5}); + + % + % Get problem dimensions + % + + % Get number of variables n + if (~isempty(P)) + n = size(P, 1); + elseif (~isempty(q)) + n = length(q); + elseif (~isempty(A)) + n = size(A, 2); + else + error('The problem does not have any variables'); + end + + % Get number of constraints m + if (isempty(A)) + m = 0; + else + m = size(A, 1); + assert(size(A, 2) == n, 'Incorrect dimension of A'); + end + + % + % Create sparse matrices and full vectors if they are empty + % + + if (isempty(P)) + P = sparse(n, n); + else + P = sparse(P); + end + if (~istriu(P)) + P = triu(P); + end + if (isempty(q)) + q = zeros(n, 1); + else + q = full(q(:)); + end + + % Create proper constraints if they are not passed + if (isempty(A) && (~isempty(l) || ~isempty(u))) || ... + (~isempty(A) && (isempty(l) && isempty(u))) + error('A must be supplied together with at least one bound l or u'); + end + + if (~isempty(A) && isempty(l)) + l = -Inf(m, 1); + end + + if (~isempty(A) && isempty(u)) + u = Inf(m, 1); + end + + if (isempty(A)) + A = sparse(m, n); + l = -Inf(m, 1); + u = Inf(m, 1); + else + l = full(l(:)); + u = full(u(:)); + A = sparse(A); + end + + + % + % Check vector dimensions (not checked from the C solver) + % + assert(length(q) == n, 'Incorrect dimension of q'); + assert(length(l) == m, 'Incorrect dimension of l'); + assert(length(u) == m, 'Incorrect dimension of u'); + + % + % Convert infinity values to OSQP_INFINITY + % + u = min(u, osqp.constant('OSQP_INFTY')); + l = max(l, -osqp.constant('OSQP_INFTY')); + + + %make a settings structure from the remainder of the arguments. + %'true' means that this is a settings initialization, so all + %parameter/values are allowed. No extra inputs will result + %in default settings being passed back + theSettings = validate_settings(this,true,varargin{6:end}); + + [varargout{1:nargout}] = osqp_mex('setup', this.objectHandle, n,m,P,q,A,l,u,theSettings); +end \ No newline at end of file diff --git a/@osqp/solve.m b/@osqp/solve.m new file mode 100644 index 0000000..8165b01 --- /dev/null +++ b/@osqp/solve.m @@ -0,0 +1,13 @@ +% SPDX-License-Identifier: Apache-2.0 + +%% +function varargout = solve(this, varargin) + % SOLVE solve the QP + + nargoutchk(0,1); %either return nothing (but still solve), or a single output structure + [out.x, out.y, out.prim_inf_cert, out.dual_inf_cert, out.info] = osqp_mex('solve', this.objectHandle); + if(nargout) + varargout{1} = out; + end + return; +end \ No newline at end of file diff --git a/@osqp/update.m b/@osqp/update.m new file mode 100644 index 0000000..4ce0903 --- /dev/null +++ b/@osqp/update.m @@ -0,0 +1,71 @@ +% SPDX-License-Identifier: Apache-2.0 + +%% +function update(this,varargin) + % UPDATE modify the linear cost term and/or lower and upper bounds + + %second input 'false' means that this is *not* a settings + %initialization, so some parameter/values will be disallowed + allowedFields = {'q','l','u','Px','Px_idx','Ax','Ax_idx'}; + + if(isempty(varargin)) + return; + elseif(length(varargin) == 1) + if(~isstruct(varargin{1})) + error('Single input should be a structure with new problem data'); + else + newData = varargin{1}; + end + else % param / value style assumed + newData = struct(varargin{:}); + end + + %check for unknown fields + newFields = fieldnames(newData); + badFieldsIdx = find(~ismember(newFields,allowedFields)); + if(~isempty(badFieldsIdx)) + error('Unrecognized input field ''%s'' detected',newFields{badFieldsIdx(1)}); + end + + %get all of the terms. Nonexistent fields will be passed + %as empty mxArrays + try q = double(full(newData.q(:))); catch q = []; end + try l = double(full(newData.l(:))); catch l = []; end + try u = double(full(newData.u(:))); catch u = []; end + try Px = double(full(newData.Px(:))); catch Px = []; end + try Px_idx = double(full(newData.Px_idx(:))); catch Px_idx = []; end + try Ax = double(full(newData.Ax(:))); catch Ax = []; end + try Ax_idx = double(full(newData.Ax_idx(:))); catch Ax_idx = []; end + + [n,m] = get_dimensions(this); + + assert(isempty(q) || length(q) == n, 'input ''q'' is the wrong size'); + assert(isempty(l) || length(l) == m, 'input ''u'' is the wrong size'); + assert(isempty(u) || length(u) == m, 'input ''l'' is the wrong size'); + assert(isempty(Px) || isempty(Px_idx) || length(Px) == length(Px_idx), ... + 'inputs ''Px'' and ''Px_idx'' must be the same size'); + assert(isempty(Ax) || isempty(Ax_idx) || length(Ax) == length(Ax_idx), ... + 'inputs ''Ax'' and ''Ax_idx'' must be the same size'); + + % Adjust index of Px_idx and Ax_idx to match 0-based indexing + % in C + if (~isempty(Px_idx)) + Px_idx = Px_idx - 1; + end + if (~isempty(Ax_idx)) + Ax_idx = Ax_idx - 1; + end + + % Convert infinity values to OSQP_INFTY + if (~isempty(u)) + u = min(u, osqp.constant('OSQP_INFTY')); + end + if (~isempty(l)) + l = max(l, -osqp.constant('OSQP_INFTY')); + end + + %write the new problem data. C-mex does not protect + %against unknown fields, but will handle empty values + osqp_mex('update', this.objectHandle, ... + q, l, u, Px, Px_idx, length(Px), Ax, Ax_idx, length(Ax)); +end \ No newline at end of file diff --git a/@osqp/update_codegen_defines.m b/@osqp/update_codegen_defines.m new file mode 100644 index 0000000..4ec8aa9 --- /dev/null +++ b/@osqp/update_codegen_defines.m @@ -0,0 +1,18 @@ +% SPDX-License-Identifier: Apache-2.0 + +function update_codegen_defines(this, varargin) + % UPDATE_CODEGEN_DEFINES update the current codegen defines + + % Check for structure style input + if(isstruct(varargin{1})) + newSettings = varargin{1}; + assert(length(varargin) == 1, 'too many input arguments'); + else + newSettings = struct(varargin{:}); + end + + % Write the new codegen defiens. The C-function checks for input + % validity. + osqp_mex('update_codegen_defines', this.objectHandle, newSettings); +end + diff --git a/@osqp/update_settings.m b/@osqp/update_settings.m new file mode 100644 index 0000000..50ec8d6 --- /dev/null +++ b/@osqp/update_settings.m @@ -0,0 +1,27 @@ +% SPDX-License-Identifier: Apache-2.0 + +function update_settings(this, varargin) + % UPDATE_SETTINGS update the current solver settings structure + + % Check for structure style input + if(isstruct(varargin{1})) + newSettings = varargin{1}; + assert(length(varargin) == 1, 'too many input arguments'); + else + newSettings = struct(varargin{:}); + end + + % Rho update must be handled special + if( isfield(newSettings, 'rho') ) + osqp_mex('update_rho', this.objectHandle, newSettings.rho); + newSettings = rmfield(newSettings, 'rho'); + end + + % Second input 'false' means that this is *not* a settings + % initialization, so some parameter/values will be disallowed + newSettings = validate_settings(this, false, varargin{:}); + + % Write the solver settings. C-mex does not check input + % data or protect against disallowed parameter modifications + osqp_mex('update_settings', this.objectHandle, newSettings); +end \ No newline at end of file diff --git a/@osqp/validate_settings.m b/@osqp/validate_settings.m new file mode 100644 index 0000000..665673a --- /dev/null +++ b/@osqp/validate_settings.m @@ -0,0 +1,73 @@ +% SPDX-License-Identifier: Apache-2.0 + +function currentSettings = validate_settings(this, isInitialization, varargin) + % Don't allow these fields to be changed + unmodifiableFields = {'scaling', 'linsys_solver'}; + + % Get the current settings + if(isInitialization) + currentSettings = osqp_mex('default_settings', this.objectHandle); + else + currentSettings = osqp_mex('current_settings', this.objectHandle); + end + + % No settings passed -> return defaults + if(isempty(varargin)) + return; + end + + % Check for structure style input + if(isstruct(varargin{1})) + newSettings = varargin{1}; + assert(length(varargin) == 1, 'too many input arguments'); + else + newSettings = struct(varargin{:}); + end + + % Get the osqp settings fields + currentFields = fieldnames(currentSettings); + + % Get the requested fields in the update + newFields = fieldnames(newSettings); + + % Check for unknown parameters + badFieldsIdx = find(~ismember(newFields,currentFields)); + if(~isempty(badFieldsIdx)) + error('Unrecognized solver setting ''%s'' detected',newFields{badFieldsIdx(1)}); + end + + % Convert linsys_solver string to integer + if ismember('linsys_solver',newFields) + if ~ischar(newSettings.linsys_solver) + error('Setting linsys_solver is required to be a string.'); + end + % Convert linsys_solver to number + newSettings.linsys_solver = string_to_linsys_solver(newSettings.linsys_solver); + end + + + % Check for disallowed fields if this in not an initialization call + if(~isInitialization) + badFieldsIdx = find(ismember(newFields,unmodifiableFields)); + for i = badFieldsIdx(:)' + if(~isequal(newSettings.(newFields{i}),currentSettings.(newFields{i}))) + error('Solver setting ''%s'' can only be changed at solver initialization.', newFields{i}); + end + end + end + + + % Check that everything is a nonnegative scalar (this check is already + % performed in C) + % for i = 1:length(newFields) + % val = double(newSettings.(newFields{i})); + % assert(isscalar(val) & isnumeric(val) & val >= 0, ... + % 'Solver setting ''%s'' not specified as nonnegative scalar', newFields{i}); + % end + + % Everything checks out - merge the newSettings into the current ones + for i = 1:length(newFields) + currentSettings.(newFields{i}) = double(newSettings.(newFields{i})); + end +end + \ No newline at end of file diff --git a/@osqp/warm_start.m b/@osqp/warm_start.m new file mode 100644 index 0000000..1fc9da4 --- /dev/null +++ b/@osqp/warm_start.m @@ -0,0 +1,52 @@ +% SPDX-License-Identifier: Apache-2.0 + +function warm_start(this, varargin) + % WARM_START warm start primal and/or dual variables + % + % warm_start('x', x, 'y', y) + % + % or warm_start('x', x) + % or warm_start('y', y) + + + % Get problem dimensions + [n, m] = get_dimensions(this); + + % Get data + allowedFields = {'x','y'}; + + if(isempty(varargin)) + return; + elseif(length(varargin) == 1) + if(~isstruct(varargin{1})) + error('Single input should be a structure with new problem data'); + else + newData = varargin{1}; + end + else % param / value style assumed + newData = struct(varargin{:}); + end + + %check for unknown fields + newFields = fieldnames(newData); + badFieldsIdx = find(~ismember(newFields,allowedFields)); + if(~isempty(badFieldsIdx)) + error('Unrecognized input field ''%s'' detected',newFields{badFieldsIdx(1)}); + end + + %get all of the terms. Nonexistent fields will be passed + %as empty mxArrays + try x = double(full(newData.x(:))); catch x = []; end + try y = double(full(newData.y(:))); catch y = []; end + + % Check dimensions + assert(isempty(x) || length(x) == n, 'input ''x'' is the wrong size'); + assert(isempty(y) || length(y) == m, 'input ''y'' is the wrong size'); + + % Only call when there is a vector to update + if (~isempty(x) || ~isempty(y)) + osqp_mex('warm_start', this.objectHandle, x, y); + else + error('Unrecognized fields'); + end +end \ No newline at end of file diff --git a/c_sources/CMakeLists.txt b/c_sources/CMakeLists.txt new file mode 100644 index 0000000..c552db2 --- /dev/null +++ b/c_sources/CMakeLists.txt @@ -0,0 +1,76 @@ +cmake_minimum_required( VERSION 3.18 ) + +project( osqp-matlab ) + + +# Find the MATLAB installation directory and various settings for it +find_package( Matlab REQUIRED COMPONENTS MEX_COMPILER) + +message( STATUS "Matlab root is " ${Matlab_ROOT_DIR} ) + +include_directories( ${Matlab_INCLUDE_DIRS} ) + +# The mex interface uses C++11 +set( CMAKE_CXX_STANDARD 11 ) +set( CMAKE_CXX_STANDARD_REQUIRED ON ) + +if( CMAKE_COMPILER_IS_GNUCXX ) + # Add debug symbols and optimizations to the build + set( CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -g -O2" ) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O2" ) +endif() + +# Insist on the pre-2018 complex data API so that mxGetPr will work correctly + add_compile_definitions( MATLAB_MEXSRC_RELEASE=R2017b ) +message( STATUS "Using Matlab pre-2018a API for mxGetPr compatibility" ) + +# Some parts of the main libraries need to know we are building for MATLAB +# So pass a definition to let them know +add_compile_definitions( MATLAB ) + +# Configure the settings of the OSQP library +set( OSQP_ALGEBRA_BACKEND "builtin" ) +set( OSQP_USE_FLOAT OFF ) +set( OSQP_CODEGEN ON ) +set( OSQP_ENABLE_DERIVATIVES ON ) +set( OSQP_ENABLE_PRINTING ON ) +set( OSQP_ENABLE_PROFILING ON ) +set( OSQP_ENABLE_INTERRUPT ON ) + +# Configure the build outputs for the OSQP library +set( OSQP_BUILD_STATIC_LIB ON ) +set( OSQP_BUILD_SHARED_LIB OFF ) +set( OSQP_BUILD_UNITTESTS OFF ) +set( OSQP_BUILD_DEMO_EXE OFF ) + +# Add the custom MATLAB memory handling +set( OSQP_CUSTOM_MEMORY "memory_matlab.h" ) +include_directories( ${CMAKE_CURRENT_SOURCE_DIR} ) + +# Actually pull in the OSQP targets +add_subdirectory( osqp_sources ) + +# Add OSQP include directories that are needed +include_directories( ${CMAKE_CURRENT_SOURCE_DIR}/osqp_sources/include/public + ${CMAKE_CURRENT_SOURCE_DIR}/osqp_sources/include/private ) + +# Interrupt support requires the libut library +get_filename_component( Matlab_LIB_DIR ${Matlab_MEX_LIBRARY} DIRECTORY ) + +if( WIN32 ) + set( UT_LIBRARY "${Matlab_LIB_DIR}/libut.lib" ) +else() + set( UT_LIBRARY "-L${Matlab_LIB_DIR} -lut" ) +endif() + +matlab_add_mex( NAME osqp_mex + SRC ${CMAKE_CURRENT_SOURCE_DIR}/osqp_mex.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/interrupt_matlab.c + ${CMAKE_CURRENT_SOURCE_DIR}/memory_matlab.c + ${CMAKE_CURRENT_SOURCE_DIR}/osqp_struct_info.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/osqp_struct_settings.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/osqp_struct_codegen_defines.cpp + LINK_TO osqpstatic + ${UT_LIBRARY} + # Force compilation in the traditional C API (equivalent to the -R2017b flag) + R2017b ) \ No newline at end of file diff --git a/c_sources/arrays_matlab.h b/c_sources/arrays_matlab.h new file mode 100644 index 0000000..3a76fa4 --- /dev/null +++ b/c_sources/arrays_matlab.h @@ -0,0 +1,45 @@ +// SPDX-License-Identifier: Apache-2.0 + +#ifndef ARRAYS_MATLAB_H_ +#define ARRAYS_MATLAB_H_ + +#include + +/** + * Copy the data from one array to another provided array. + */ +template +void copyVector(outArr* out, inArr* in, OSQPInt numel) { + // Don't bother doing anything if there is no input data + if(!in || !out || (numel == 0)) + return; + + // Copy the data + for(OSQPInt i=0; i < numel; i++){ + out[i] = static_cast(in[i]); + } +} + + +/** + * Copy the data from one array to another newly allocated array. + * The caller gains ownership of the returned array. + */ +template +outArr* cloneVector(inArr* in, OSQPInt numel) { + // Don't bother doing anything if there is no input data + if(!in || (numel == 0)) + return NULL; + + // Allocate new array + outArr* out = static_cast(c_malloc(numel * sizeof(outArr))); + + if(!out) + mexErrMsgTxt("Failed to allocate a vector object."); + + // Copy the data + copyVector(out, in, numel); + return out; +} + +#endif \ No newline at end of file diff --git a/c_sources/interrupt_matlab.c b/c_sources/interrupt_matlab.c new file mode 100644 index 0000000..e61da20 --- /dev/null +++ b/c_sources/interrupt_matlab.c @@ -0,0 +1,27 @@ +// SPDX-License-Identifier: Apache-2.0 + +/* + * Implements interrupt handling using ctrl-c for MATLAB mex files. + */ + +#include "interrupt.h" + +#include + +/* No header file available here; define the prototypes ourselves */ +bool utIsInterruptPending(void); +bool utSetInterruptEnabled(bool); + +static int istate; + +void osqp_start_interrupt_listener(void) { + istate = utSetInterruptEnabled(1); +} + +void osqp_end_interrupt_listener(void) { + utSetInterruptEnabled(istate); +} + +int osqp_is_interrupted(void) { + return utIsInterruptPending(); +} diff --git a/c_sources/memory_matlab.c b/c_sources/memory_matlab.c new file mode 100644 index 0000000..0c25bb9 --- /dev/null +++ b/c_sources/memory_matlab.c @@ -0,0 +1,21 @@ +// SPDX-License-Identifier: Apache-2.0 + +#include + +void* c_calloc(size_t num, size_t size) { + void *m = mxCalloc(num, size); + mexMakeMemoryPersistent(m); + return m; +} + +void* c_malloc(size_t size) { + void *m = mxMalloc(size); + mexMakeMemoryPersistent(m); + return m; +} + +void* c_realloc(void *ptr, size_t size) { + void *m = mxRealloc(ptr, size); + mexMakeMemoryPersistent(m); + return m; +} \ No newline at end of file diff --git a/c_sources/memory_matlab.h b/c_sources/memory_matlab.h new file mode 100644 index 0000000..b63a0cf --- /dev/null +++ b/c_sources/memory_matlab.h @@ -0,0 +1,19 @@ +// SPDX-License-Identifier: Apache-2.0 + +/* Memory managment for MATLAB */ +#include "mex.h" + + +#ifdef __cplusplus +extern "C" { +#endif + + void* c_calloc(size_t num, size_t size); + void* c_malloc(size_t size); + void* c_realloc(void *ptr, size_t size); + +#ifdef __cplusplus +} +#endif + +#define c_free mxFree \ No newline at end of file diff --git a/c_sources/osqp_mex.cpp b/c_sources/osqp_mex.cpp new file mode 100644 index 0000000..69ec67c --- /dev/null +++ b/c_sources/osqp_mex.cpp @@ -0,0 +1,644 @@ +// SPDX-License-Identifier: Apache-2.0 + +#include +#include + +#include "mex.h" +#include "matrix.h" +#include "osqp.h" +// Mex-specific functionality +#include "osqp_mex.hpp" +#include "osqp_struct.h" +#include "arrays_matlab.h" +#include "memory_matlab.h" + +//TODO: Check if this definition is required, and maybe replace it with: +// enum linsys_solver_type { QDLDL_SOLVER, MKL_PARDISO_SOLVER }; +#define QDLDL_SOLVER 0 //Based on the previous API +#define CMD_MAX_LEN 64 +#define OUTPUT_DIR_MAX_LEN 256 +#define PREFIX_MAX_LEN 128 + +using std::string; + +// Wrapper class to pass the OSQP solver back and forth with Matlab +class OsqpData +{ +public: + OsqpData() : + solver(NULL), defines(NULL) + {} + OSQPSolver* solver; + OSQPCodegenDefines* defines; +}; + + +// Internal utility function +static void setToNaN(double* arr_out, OSQPInt len){ + OSQPInt i; + for (i = 0; i < len; i++) { + arr_out[i] = mxGetNaN(); + } +} + +// This is a utility function that uses mexErrMsgTxt using std::string instead of char*. +void mexErrMsgTxt(string str) { + mexErrMsgTxt(str.c_str()); +} + +/** + * This is a utility function that assigns prhs[ind] to char. + * + * @param nrhs Number of input arguments + * @param prhs Matlab input arrays + * @param ind prhs index + * @param str Target string + * @param max_len Maximum allowed length of the string. If passed 0 (default), use sizeof(str) instead. +*/ +void setString(int nrhs, const mxArray *prhs[], int ind, char *str, int max_len){ + int str_max_len = sizeof(str); + if (max_len>0) str_max_len = max_len; + if (nrhs < (ind+1) || mxGetString(prhs[ind], str, str_max_len)){ + mexErrMsgTxt("Input #" + std::to_string(ind) + " should be less than " + std::to_string(str_max_len) + " characters long."); + } +} + +/** + * This function validates the inputs to OSQPCodegenDefines. Calls mexErrMsgTxt if a NULL or an invalid input is passed. + * + * @param defines OSQPCodegenDefines object +*/ +void validateCodegenDefines(const OSQPCodegenDefines* defines) { + if (!defines) mexErrMsgTxt("A NULL defines object has been passed."); + + if (defines->embedded_mode){ + if (defines->embedded_mode != 1 && defines->embedded_mode != 2) mexErrMsgTxt("Invalid embedded mode."); + } + + if (defines->float_type){ + if (defines->float_type != 0 && defines->float_type != 1) mexErrMsgTxt("Invalid float type."); + } + + if (defines->printing_enable){ + if (defines->printing_enable != 0 && defines->printing_enable != 1) mexErrMsgTxt("Invalid printing enable."); + } + + if (defines->profiling_enable) { + if (defines->profiling_enable != 0 && defines->profiling_enable != 1) mexErrMsgTxt("Invalid profiling enable."); + } + + if (defines->interrupt_enable) { + if (defines->interrupt_enable != 0 && defines->interrupt_enable != 1) mexErrMsgTxt("Invalid interrupt enable."); + } + + if (defines->derivatives_enable) { + if (defines->derivatives_enable != 0 && defines->derivatives_enable != 1) mexErrMsgTxt("Invalid derivatives enable."); + } +} + + + +/** + * This function updates CodegenDefines object. Calls mexErrMsgTxt if a NULL or an invalid input is passed. + * + * @param target_defines Pointer to the target OSQPCodegenDefines object. + * @param new_defines Pointer to the OSQPCodegenDefines object witht he new values. +*/ +void updateCodegenDefines(OSQPCodegenDefines* target_defines, + const OSQPCodegenDefines* new_defines) { + if (!target_defines) mexErrMsgTxt("Defines is uninitialized. No codegen defines have been configured."); + validateCodegenDefines(new_defines); + + if (new_defines->embedded_mode){ + target_defines->embedded_mode = new_defines->embedded_mode; + } + + if (new_defines->float_type){ + target_defines->float_type = new_defines->float_type; + } + + if (new_defines->printing_enable){ + target_defines->printing_enable = new_defines->printing_enable; + } + + if (new_defines->profiling_enable) { + target_defines->profiling_enable = new_defines->profiling_enable; + } + + if (new_defines->interrupt_enable) { + target_defines->interrupt_enable = new_defines->interrupt_enable; + } + + if (new_defines->derivatives_enable) { + target_defines->derivatives_enable = new_defines->derivatives_enable; + } +} + +// Main mex function +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + // OSQP solver wrapper + OsqpData* osqpData; + + // Exitflag + OSQPInt exitflag = 0; + + // Get the command string + char cmd[CMD_MAX_LEN]; + + setString(nrhs, prhs, 0, cmd, CMD_MAX_LEN); + + /* + * First check to see if a new object was requested + */ + if (!strcmp("new", cmd)) { + // Check parameters + if (nlhs != 1){ + mexErrMsgTxt("New: One output expected."); + } + // Return a handle to a new C++ wrapper instance + osqpData = new OsqpData; + osqpData->defines = new OSQPCodegenDefines; + osqp_set_default_codegen_defines(osqpData->defines); + plhs[0] = convertPtr2Mat(osqpData); + return; + } + + /* + * Next check to see if any of the static methods were called + */ + // Report the version + if (!strcmp("version", cmd)) { + plhs[0] = mxCreateString(osqp_version()); + return; + } + + // Report the default settings + if (!strcmp("default_settings", cmd)) { + // Warn if other commands were ignored + if (nrhs > 2) + mexWarnMsgTxt("Default settings: unexpected number of arguments."); + + // Create a Settings structure in default form and report the results + // Useful for external solver packages (e.g. Yalmip) that want to + // know which solver settings are supported + OSQPSettingsWrapper settings; + plhs[0] = settings.GetMxStruct(); + return; + } + + // Return solver constants + if (!strcmp("constant", cmd)) { + static std::map floatConstants{ + // Numerical constants + {"OSQP_INFTY", OSQP_INFTY} + }; + + static std::map intConstants{ + // Return codes + {"OSQP_SOLVED", OSQP_SOLVED}, + {"OSQP_SOLVED_INACCURATE", OSQP_SOLVED_INACCURATE}, + {"OSQP_UNSOLVED", OSQP_UNSOLVED}, + {"OSQP_PRIMAL_INFEASIBLE", OSQP_PRIMAL_INFEASIBLE}, + {"OSQP_PRIMAL_INFEASIBLE_INACCURATE", OSQP_PRIMAL_INFEASIBLE_INACCURATE}, + {"OSQP_DUAL_INFEASIBLE", OSQP_DUAL_INFEASIBLE}, + {"OSQP_DUAL_INFEASIBLE_INACCURATE", OSQP_DUAL_INFEASIBLE_INACCURATE}, + {"OSQP_MAX_ITER_REACHED", OSQP_MAX_ITER_REACHED}, + {"OSQP_NON_CVX", OSQP_NON_CVX}, + {"OSQP_TIME_LIMIT_REACHED", OSQP_TIME_LIMIT_REACHED}, + + // Linear system solvers + {"QDLDL_SOLVER", QDLDL_SOLVER}, + {"OSQP_UNKNOWN_SOLVER", OSQP_UNKNOWN_SOLVER}, + {"OSQP_DIRECT_SOLVER", OSQP_DIRECT_SOLVER}, + {"OSQP_INDIRECT_SOLVER", OSQP_INDIRECT_SOLVER} + }; + + char constant[64]; + int constantLength = mxGetN(prhs[1]) + 1; + setString(nrhs, prhs, 1, constant, constantLength); + + auto ci = intConstants.find(constant); + + if(ci != intConstants.end()) { + plhs[0] = mxCreateDoubleScalar(ci->second); + return; + } + + auto cf = floatConstants.find(constant); + + if(cf != floatConstants.end()) { + plhs[0] = mxCreateDoubleScalar(cf->second); + return; + } + + // NaN is special because we need the Matlab version + if (!strcmp("OSQP_NAN", constant)){ + plhs[0] = mxCreateDoubleScalar(mxGetNaN()); + return; + } + + mexErrMsgTxt("Constant not recognized."); + + return; + } + + /* + * Finally, check to see if this is a function operating on a solver instance + */ + + // Check for a second input, which should be the class instance handle + if (nrhs < 2) mexErrMsgTxt("Second input should be a class instance handle."); + + // Get the class instance pointer from the second input + osqpData = convertMat2Ptr(prhs[1]); + + // delete the object and its data + if (!strcmp("delete", cmd)) { + osqp_cleanup(osqpData->solver); + destroyObject(prhs[1]); + // Warn if other commands were ignored + if (nlhs != 0 || nrhs != 2) + mexWarnMsgTxt("Delete: Unexpected arguments ignored."); + return; + } + + // report the current settings + if (!strcmp("current_settings", cmd)) { + // Throw an error if this is called before solver is configured + if(!osqpData->solver) { + mexErrMsgTxt("Solver is uninitialized. No settings have been configured."); + } + if(!osqpData->solver->settings) { + mexErrMsgTxt("Solver settings is uninitialized. No settings have been configured."); + } + + // Report the current settings + OSQPSettingsWrapper settings(osqpData->solver->settings); + plhs[0] = settings.GetMxStruct(); + return; + } + + // write_settings + if (!strcmp("update_settings", cmd)) { + // Overwrite the current settings. Mex function is responsible + // for disallowing overwrite of selected settings after initialization, + // and for all error checking + // throw an error if this is called before solver is configured + if(!osqpData->solver){ + mexErrMsgTxt("Solver is uninitialized. No settings have been configured."); + } + + OSQPSettingsWrapper settings(prhs[2]); + osqp_update_settings(osqpData->solver, settings.GetOSQPStruct()); + return; + } + + // Update rho value + if (!strcmp("update_rho", cmd)) { + //throw an error if this is called before solver is configured + if(!osqpData->solver){ + mexErrMsgTxt("Solver is uninitialized. No settings have been configured."); + } + + OSQPFloat rho = (OSQPFloat)mxGetScalar(prhs[2]); + + osqp_update_rho(osqpData->solver, rho); + return; + } + + // setup + if (!strcmp("setup", cmd)) { + //throw an error if this is called more than once + if(osqpData->solver){ + mexErrMsgTxt("Solver is already initialized with problem data."); + } + + // handle the problem data first. Matlab-side + // class wrapper is responsible for ensuring that + // P and A are sparse matrices, everything + // else is a dense vector and all inputs are + // of compatible dimension + + // Get mxArrays + const mxArray* P = prhs[4]; + const mxArray* q = prhs[5]; + const mxArray* A = prhs[6]; + const mxArray* l = prhs[7]; + const mxArray* u = prhs[8]; + + + OSQPInt dataN = (OSQPInt)mxGetScalar(prhs[2]); + OSQPInt dataM = (OSQPInt)mxGetScalar(prhs[3]); + OSQPFloat* dataQ = cloneVector(mxGetPr(q), dataN); + OSQPFloat* dataL = cloneVector(mxGetPr(l), dataM); + OSQPFloat* dataU = cloneVector(mxGetPr(u), dataM); + + // Matrix P: nnz = P->p[n] + OSQPInt * Pp = cloneVector(mxGetJc(P), dataN + 1); + OSQPInt * Pi = cloneVector(mxGetIr(P), Pp[dataN]); + OSQPFloat * Px = cloneVector(mxGetPr(P), Pp[dataN]); + OSQPCscMatrix* dataP = OSQPCscMatrix_new(dataN, dataN, Pp[dataN], Px, Pi, Pp); + + // Matrix A: nnz = A->p[n] + OSQPInt* Ap = cloneVector(mxGetJc(A), dataN + 1); + OSQPInt* Ai = cloneVector(mxGetIr(A), Ap[dataN]); + OSQPFloat * Ax = cloneVector(mxGetPr(A), Ap[dataN]); + OSQPCscMatrix* dataA = OSQPCscMatrix_new(dataM, dataN, Ap[dataN], Ax, Ai, Ap); + + // Create Settings + OSQPSettingsWrapper settings; + + if(!mxIsEmpty(prhs[9])){ + // Populate settings structure from mxArray input, otherwise the default settings are used + settings.ParseMxStruct(prhs[9]); + } + + // Setup workspace + //exitflag = osqp_setup(&(osqpData->work), data, settings); + exitflag = osqp_setup(&(osqpData->solver), dataP, dataQ, dataA, dataL, dataU, dataM, dataN, settings.GetOSQPStruct()); + //cleanup temporary structures + // Data + if (Px) c_free(Px); + if (Pi) c_free(Pi); + if (Pp) c_free(Pp); + if (Ax) c_free(Ax); + if (Ai) c_free(Ai); + if (Ap) c_free(Ap); + if (dataQ) c_free(dataQ); + if (dataL) c_free(dataL); + if (dataU) c_free(dataU); + if (dataP) OSQPCscMatrix_free(dataP); + if (dataA) OSQPCscMatrix_free(dataA); + + // Report error (if any) + if(exitflag){ + mexErrMsgTxt("Invalid problem setup"); + } + + return; + + } + + // get #constraints and variables + if (!strcmp("get_dimensions", cmd)) { + + //throw an error if this is called before solver is configured + if(!osqpData->solver){ + mexErrMsgTxt("Solver has not been initialized."); + } + OSQPInt n, m; + osqp_get_dimensions(osqpData->solver, &m, &n); + plhs[0] = mxCreateDoubleScalar(n); + plhs[1] = mxCreateDoubleScalar(m); + + return; + } + + // update linear cost and bounds + if (!strcmp("update", cmd)) { + + //throw an error if this is called before solver is configured + if(!osqpData->solver){ + mexErrMsgTxt("Solver has not been initialized."); + } + + // Fill q, l, u + const mxArray *q = prhs[2]; + const mxArray *l = prhs[3]; + const mxArray *u = prhs[4]; + const mxArray *Px = prhs[5]; + const mxArray *Px_idx = prhs[6]; + const mxArray *Ax = prhs[8]; + const mxArray *Ax_idx = prhs[9]; + + int Px_n, Ax_n; + Px_n = mxGetScalar(prhs[7]); + Ax_n = mxGetScalar(prhs[10]); + + // Copy vectors to ensure they are cast as OSQPFloat + OSQPFloat *q_vec = NULL; + OSQPFloat *l_vec = NULL; + OSQPFloat *u_vec = NULL; + OSQPFloat *Px_vec = NULL; + OSQPFloat *Ax_vec = NULL; + OSQPInt *Px_idx_vec = NULL; + OSQPInt *Ax_idx_vec = NULL; + + OSQPInt n, m; + osqp_get_dimensions(osqpData->solver, &m, &n); + if(!mxIsEmpty(q)){ + q_vec = cloneVector(mxGetPr(q), n); + } + if(!mxIsEmpty(l)){ + l_vec = cloneVector(mxGetPr(l), m); + } + if(!mxIsEmpty(u)){ + u_vec = cloneVector(mxGetPr(u), m); + } + if(!mxIsEmpty(Px)){ + Px_vec = cloneVector(mxGetPr(Px), Px_n); + } + if(!mxIsEmpty(Ax)){ + Ax_vec = cloneVector(mxGetPr(Ax), Ax_n); + } + if(!mxIsEmpty(Px_idx)){ + Px_idx_vec = cloneVector(mxGetPr(Px_idx), Px_n); + } + if(!mxIsEmpty(Ax_idx)){ + Ax_idx_vec = cloneVector(mxGetPr(Ax_idx), Ax_n); + } + + if (!exitflag && (!mxIsEmpty(q) || !mxIsEmpty(l) || !mxIsEmpty(u))) { + exitflag = osqp_update_data_vec(osqpData->solver, q_vec, l_vec, u_vec); + if (exitflag) exitflag=1; + } + + if (!exitflag && (!mxIsEmpty(Px) || !mxIsEmpty(Ax))) { + exitflag = osqp_update_data_mat(osqpData->solver, Px_vec, Px_idx_vec, Px_n, Ax_vec, Ax_idx_vec, Ax_n); + if (exitflag) exitflag=2; + } + + + // Free vectors + if(!mxIsEmpty(q)) c_free(q_vec); + if(!mxIsEmpty(l)) c_free(l_vec); + if(!mxIsEmpty(u)) c_free(u_vec); + if(!mxIsEmpty(Px)) c_free(Px_vec); + if(!mxIsEmpty(Ax)) c_free(Ax_vec); + if(!mxIsEmpty(Px_idx)) c_free(Px_idx_vec); + if(!mxIsEmpty(Ax_idx)) c_free(Ax_idx_vec); + + // Report errors (if any) + switch (exitflag) { + case 1: + mexErrMsgTxt("Data update error!"); + case 2: + mexErrMsgTxt("Matrix update error!"); + } + + return; + } + + if (!strcmp("warm_start", cmd)) { + + //throw an error if this is called before solver is configured + if(!osqpData->solver){ + mexErrMsgTxt("Solver has not been initialized."); + } + + // Fill x and y + const mxArray *x = prhs[2]; + const mxArray *y = prhs[3]; + + // Copy vectors to ensure they are cast as OSQPFloat + OSQPFloat *x_vec = NULL; + OSQPFloat *y_vec = NULL; + + OSQPInt n, m; + osqp_get_dimensions(osqpData->solver, &m, &n); + + if(!mxIsEmpty(x)){ + x_vec = cloneVector(mxGetPr(x),n); + } + if(!mxIsEmpty(y)){ + y_vec = cloneVector(mxGetPr(y),m); + } + + // Warm start x and y + osqp_warm_start(osqpData->solver, x_vec, y_vec); + + // Free vectors + if(!x_vec) c_free(x_vec); + if(!y_vec) c_free(y_vec); + + return; + } + + // SOLVE + if (!strcmp("solve", cmd)) { + if (nlhs != 5 || nrhs != 2){ + mexErrMsgTxt("Solve : wrong number of inputs / outputs"); + } + if(!osqpData->solver){ + mexErrMsgTxt("No problem data has been given."); + } + // solve the problem + osqp_solve(osqpData->solver); + + OSQPInt n, m; + osqp_get_dimensions(osqpData->solver, &m, &n); + // Allocate space for solution + // primal variables + plhs[0] = mxCreateDoubleMatrix(n,1,mxREAL); + // dual variables + plhs[1] = mxCreateDoubleMatrix(m,1,mxREAL); + // primal infeasibility certificate + plhs[2] = mxCreateDoubleMatrix(m,1,mxREAL); + // dual infeasibility certificate + plhs[3] = mxCreateDoubleMatrix(n,1,mxREAL); + + //copy results to mxArray outputs + //assume that five outputs will always + //be returned to matlab-side class wrapper + if ((osqpData->solver->info->status_val != OSQP_PRIMAL_INFEASIBLE) && + (osqpData->solver->info->status_val != OSQP_DUAL_INFEASIBLE)){ + + //primal and dual solutions + copyVector(mxGetPr(plhs[0]), osqpData->solver->solution->x, n); + copyVector(mxGetPr(plhs[1]), osqpData->solver->solution->y, m); + + //infeasibility certificates -> NaN values + setToNaN(mxGetPr(plhs[2]), m); + setToNaN(mxGetPr(plhs[3]), n); + + } else if (osqpData->solver->info->status_val == OSQP_PRIMAL_INFEASIBLE || + osqpData->solver->info->status_val == OSQP_PRIMAL_INFEASIBLE_INACCURATE){ //primal infeasible + + //primal and dual solutions -> NaN values + setToNaN(mxGetPr(plhs[0]), n); + setToNaN(mxGetPr(plhs[1]), m); + + //primal infeasibility certificates + copyVector(mxGetPr(plhs[2]), osqpData->solver->solution->prim_inf_cert, m); + + //dual infeasibility certificates -> NaN values + setToNaN(mxGetPr(plhs[3]), n); + + // Set objective value to infinity + osqpData->solver->info->obj_val = mxGetInf(); + + } else { //dual infeasible + + //primal and dual solutions -> NaN values + setToNaN(mxGetPr(plhs[0]), n); + setToNaN(mxGetPr(plhs[1]), m); + + //primal infeasibility certificates -> NaN values + setToNaN(mxGetPr(plhs[2]), m); + + //dual infeasibility certificates + copyVector(mxGetPr(plhs[3]), osqpData->solver->solution->dual_inf_cert, n); + + // Set objective value to -infinity + osqpData->solver->info->obj_val = -mxGetInf(); + } + + if (osqpData->solver->info->status_val == OSQP_NON_CVX) { + osqpData->solver->info->obj_val = mxGetNaN(); + } + + // Populate the info structure + OSQPInfoWrapper info(osqpData->solver->info); + plhs[4] = info.GetMxStruct(); + + return; + } + + if (!strcmp("default_codegen_defines", cmd)) { + // Warn if other commands were ignored + if (nrhs > 2) + mexWarnMsgTxt("Default codegen settings: unexpected number of arguments."); + + // Create a Settings structure in default form and report the results + OSQPCodegenDefinesWrapper codegen_settings; + plhs[0] = codegen_settings.GetMxStruct(); + return; + } + + if (!strcmp("update_codegen_defines", cmd)) { + // Overwrite the current settings. Mex function is responsible + // for disallowing overwrite of selected settings after initialization, + // and for all error checking + // throw an error if this is called before solver is configured + if(!osqpData->defines){ + mexErrMsgTxt("Defines is uninitialized. No codegen defines have been configured."); + } + + OSQPCodegenDefinesWrapper defines(prhs[2]); + updateCodegenDefines(osqpData->defines, defines.GetOSQPStruct()); + return; + } + + if (!strcmp("codegen", cmd)) { + + //Check that the correct number of arguments is passed + if (nrhs != 4) mexErrMsgTxt("Codegen: unexpected number of arguments"); + char output_dir[OUTPUT_DIR_MAX_LEN]; + char prefix[PREFIX_MAX_LEN]; + setString(nrhs, prhs, 2, output_dir, OUTPUT_DIR_MAX_LEN); + setString(nrhs, prhs, 3, prefix, PREFIX_MAX_LEN); + + //Check that the solver was initialized + if(!osqpData->solver) mexErrMsgTxt("Solver has not been initialized."); + if(!osqpData->defines) mexErrMsgTxt("Codegen defines has not been initialized."); + + exitflag = osqp_codegen(osqpData->solver,output_dir, prefix, osqpData->defines); + + if (exitflag) mexErrMsgTxt("Codegen failed with exitflag = " + std::to_string(exitflag)); + return; + } + + // Got here, so command not recognized + mexErrMsgTxt("Command not recognized."); +} diff --git a/osqp_mex.hpp b/c_sources/osqp_mex.hpp similarity index 96% rename from osqp_mex.hpp rename to c_sources/osqp_mex.hpp index ceb3b07..e0d1c0f 100755 --- a/osqp_mex.hpp +++ b/c_sources/osqp_mex.hpp @@ -1,53 +1,53 @@ -#ifndef __OSQP_MEX_HPP__ -#define __OSQP_MEX_HPP__ -#include "mex.h" -#include -#include -#include -#include - -#define OSQP_MEX_SIGNATURE 0x271C1A7A -template class osqp_mex_handle -{ -public: - osqp_mex_handle(base *ptr) : ptr_m(ptr), name_m(typeid(base).name()) { signature_m = OSQP_MEX_SIGNATURE; } - ~osqp_mex_handle() { signature_m = 0; delete ptr_m; } - bool isValid() { return ((signature_m == OSQP_MEX_SIGNATURE) && !strcmp(name_m.c_str(), typeid(base).name())); } - base *ptr() { return ptr_m; } - -private: - uint32_t signature_m; - std::string name_m; - base *ptr_m; -}; - -template inline mxArray *convertPtr2Mat(base *ptr) -{ - mexLock(); - mxArray *out = mxCreateNumericMatrix(1, 1, mxUINT64_CLASS, mxREAL); - *((uint64_t *)mxGetData(out)) = reinterpret_cast(new osqp_mex_handle(ptr)); - return out; -} - -template inline osqp_mex_handle *convertMat2HandlePtr(const mxArray *in) -{ - if (mxGetNumberOfElements(in) != 1 || mxGetClassID(in) != mxUINT64_CLASS || mxIsComplex(in)) - mexErrMsgTxt("Input must be a real uint64 scalar."); - osqp_mex_handle *ptr = reinterpret_cast *>(*((uint64_t *)mxGetData(in))); - if (!ptr->isValid()) - mexErrMsgTxt("Handle not valid."); - return ptr; -} - -template inline base *convertMat2Ptr(const mxArray *in) -{ - return convertMat2HandlePtr(in)->ptr(); -} - -template inline void destroyObject(const mxArray *in) -{ - delete convertMat2HandlePtr(in); - mexUnlock(); -} - -#endif // __OSQP_MEX_HPP__ +#ifndef __OSQP_MEX_HPP__ +#define __OSQP_MEX_HPP__ +#include "mex.h" +#include +#include +#include +#include + +#define OSQP_MEX_SIGNATURE 0x271C1A7A +template class osqp_mex_handle +{ +public: + osqp_mex_handle(base *ptr) : ptr_m(ptr), name_m(typeid(base).name()) { signature_m = OSQP_MEX_SIGNATURE; } + ~osqp_mex_handle() { signature_m = 0; delete ptr_m; } + bool isValid() { return ((signature_m == OSQP_MEX_SIGNATURE) && !strcmp(name_m.c_str(), typeid(base).name())); } + base *ptr() { return ptr_m; } + +private: + uint32_t signature_m; + std::string name_m; + base *ptr_m; +}; + +template inline mxArray *convertPtr2Mat(base *ptr) +{ + mexLock(); + mxArray *out = mxCreateNumericMatrix(1, 1, mxUINT64_CLASS, mxREAL); + *((uint64_t *)mxGetData(out)) = reinterpret_cast(new osqp_mex_handle(ptr)); + return out; +} + +template inline osqp_mex_handle *convertMat2HandlePtr(const mxArray *in) +{ + if (mxGetNumberOfElements(in) != 1 || mxGetClassID(in) != mxUINT64_CLASS || mxIsComplex(in)) + mexErrMsgTxt("Input must be a real uint64 scalar."); + osqp_mex_handle *ptr = reinterpret_cast *>(*((uint64_t *)mxGetData(in))); + if (!ptr->isValid()) + mexErrMsgTxt("Handle not valid."); + return ptr; +} + +template inline base *convertMat2Ptr(const mxArray *in) +{ + return convertMat2HandlePtr(in)->ptr(); +} + +template inline void destroyObject(const mxArray *in) +{ + delete convertMat2HandlePtr(in); + mexUnlock(); +} + +#endif // __OSQP_MEX_HPP__ diff --git a/c_sources/osqp_sources b/c_sources/osqp_sources new file mode 160000 index 0000000..b6bfa62 --- /dev/null +++ b/c_sources/osqp_sources @@ -0,0 +1 @@ +Subproject commit b6bfa62657828790ae088c714915db104ebd6ab3 diff --git a/c_sources/osqp_struct.h b/c_sources/osqp_struct.h new file mode 100644 index 0000000..4ee2a61 --- /dev/null +++ b/c_sources/osqp_struct.h @@ -0,0 +1,213 @@ +// SPDX-License-Identifier: Apache-2.0 + +#ifndef OSQP_STRUCT_H_ +#define OSQP_STRUCT_H_ + +#include +#include +#include +#include + +#include +#include + +#include "memory_matlab.h" +#include + +/** + * Base class used to store the field types for a struct. + */ +class OSQPStructFieldBase { +public: + OSQPStructFieldBase() {} + + /** + * Set the field in the given Matlab struct to the value of this field + */ + virtual void ToMxStruct(mxArray* aStruct) = 0; + + /** + * Set the field in the internal struct with the data from aStruct + */ + virtual void ToOSQPStruct(const mxArray* aStruct) = 0; +}; + +/** + * Class to hold a numeric struct field (e.g. float/double/int/enum, etc.). + */ +template +class OSQPStructField : public OSQPStructFieldBase { +public: + OSQPStructField(T* aStructPtr, std::string aName) : + m_structPtr(aStructPtr), + m_name(aName) { + } + + void ToMxStruct(mxArray* aStruct) override { + mxAddField(aStruct, m_name.data()); + mxSetField(aStruct, 0, m_name.data(), mxCreateDoubleScalar(*m_structPtr)); + } + + void ToOSQPStruct(const mxArray* aStruct) override { + *(m_structPtr) = static_cast(mxGetScalar(mxGetField(aStruct, 0, m_name.data()))); + } + +private: + T* m_structPtr; + std::string m_name; +}; + +/** + * Class to hold a character array (actual array, not char* array) field in a struct. + */ +class OSQPStructFieldCharArray : public OSQPStructFieldBase { +public: + OSQPStructFieldCharArray(char* aStructPtr, size_t aLength, std::string aName) : + m_structPtr(aStructPtr), + m_name(aName), + m_length(aLength) { + } + + void ToMxStruct(mxArray* aStruct) override { + mxAddField(aStruct, m_name.data()); + mxSetField(aStruct, 0, m_name.data(), mxCreateString(m_structPtr)); + } + + void ToOSQPStruct(const mxArray* aStruct) override { + mxArray* tmp = mxGetField(aStruct, 0, m_name.data()); + mxGetString(tmp, m_structPtr, m_length); + } + +private: + char* m_structPtr; + std::string m_name; + size_t m_length; +}; + +/** + * Wrap a struct from OSQP to automatically transfer the data between OSQP and Matlab. + */ +template +class OSQPStructWrapper { +public: + /** + * Initialize the wrapper using the default values. + */ + OSQPStructWrapper() { + // Allocate the default struct and register field handlers + registerFields(); + } + + /** + * Initialize the wrapper using the values from the OSQP struct. + */ + OSQPStructWrapper(const T* aStruct) { + // Allocate the default struct and register field handlers + registerFields(); + ParseOSQPStruct(aStruct); + } + + /** + * Initialize the wrapper using the values from the Matlab struct + */ + OSQPStructWrapper(const mxArray* aStruct) { + // Allocate the default struct and register field handlers + registerFields(); + ParseMxStruct(aStruct); + } + + ~OSQPStructWrapper() { + for(auto& s : m_structFields) { + delete s; + } + + c_free(m_struct); + } + + /** + * Return a Matlab struct populated with the values of the current struct + * contained in this wrapper. + * + * @return a Matlab struct with a copy of the struct (caller owns this copy and must free it) + */ + mxArray* GetMxStruct() { + // No fields are added right now, they are added in the for loop when they are set + mxArray* matStruct = mxCreateStructMatrix(1, 1, 0, NULL); + + // Copy the current struct into the struct to return + for(const auto& s : m_structFields) { + s->ToMxStruct(matStruct); + } + + return matStruct; + } + + /** + * Read a Matlab struct and populate the wrapper with its values. + */ + void ParseMxStruct(const mxArray* aStruct) { + for(const auto& s : m_structFields) { + s->ToOSQPStruct(aStruct); + } + } + + /** + * Get a copy of the struct contained inside this wrapper. + * + * @return a copy of the struct (caller owns this copy and must free it) + */ + T* GetOSQPStructCopy() { + // Allocate the default struct + T* ret = static_cast(c_calloc(1, sizeof(T))); + + // Copy the current values for their return + std::memcpy(ret, m_struct, sizeof(T)); + return ret; + } + + /** + * Get the pointer to the internal struct object. + */ + T* GetOSQPStruct() { + return m_struct; + } + + /* + * Read an existing OSQP struct object into this wrapper. + * The struct elements are copied, so no ownership of the aStruct pointer is transferred. + */ + void ParseOSQPStruct(const T* aStruct) { + std::memcpy(m_struct, aStruct, sizeof(T)); + } + +private: + /** + * Register all the fields for the wrapper. + * This function should be specialized for each struct type to map the fields appropriately. + */ + void registerFields(); + + // All struct fields + std::vector m_structFields; + + // Base OSQP struct object. Owned by this wrapper. + T* m_struct; +}; + +/** + * Wrapper around the OSQPSettings struct + */ +typedef OSQPStructWrapper OSQPSettingsWrapper; + +/** + * Wrapper around the OSQPInfo struct + */ +typedef OSQPStructWrapper OSQPInfoWrapper; + +/** + * Wrapper around the OSQPCodegenDefines struct + */ +typedef OSQPStructWrapper OSQPCodegenDefinesWrapper; + + +#endif \ No newline at end of file diff --git a/c_sources/osqp_struct_codegen_defines.cpp b/c_sources/osqp_struct_codegen_defines.cpp new file mode 100644 index 0000000..7f35a60 --- /dev/null +++ b/c_sources/osqp_struct_codegen_defines.cpp @@ -0,0 +1,30 @@ +// SPDX-License-Identifier: Apache-2.0 + +#include +#include "osqp_struct.h" + +/* + * Specialization for the codegen_defines struct + */ +template<> +void OSQPStructWrapper::registerFields() { + m_struct = static_cast(c_calloc(1, sizeof(OSQPCodegenDefines))); + if(!m_struct) + mexErrMsgTxt("Failed to allocate a OSQPCodegenDefines object."); + + osqp_set_default_codegen_defines(m_struct); + + /* + * Register the mapping between struct field name and the settings memory location + */ + m_structFields.push_back(new OSQPStructField(&m_struct->embedded_mode, "embedded_mode")); + m_structFields.push_back(new OSQPStructField(&m_struct->float_type, "float_type")); + m_structFields.push_back(new OSQPStructField(&m_struct->printing_enable, "printing_enable")); + m_structFields.push_back(new OSQPStructField(&m_struct->profiling_enable, "profiling_enable")); + m_structFields.push_back(new OSQPStructField(&m_struct->interrupt_enable, "interrupt_enable")); + m_structFields.push_back(new OSQPStructField(&m_struct->derivatives_enable,"derivatives_enable")); +} + + +// Instantiate the OSQPCodegenDefines wrapper class +template class OSQPStructWrapper; \ No newline at end of file diff --git a/c_sources/osqp_struct_info.cpp b/c_sources/osqp_struct_info.cpp new file mode 100644 index 0000000..5db5994 --- /dev/null +++ b/c_sources/osqp_struct_info.cpp @@ -0,0 +1,45 @@ +// SPDX-License-Identifier: Apache-2.0 + +#include +#include "osqp_struct.h" + + +/* + * Specialization of the struct wrapper for the OSQPInfo struct. + */ +template<> +void OSQPStructWrapper::registerFields() { + m_struct = static_cast(c_calloc(1, sizeof(OSQPInfo))); + + if(!m_struct) + mexErrMsgTxt("Failed to allocate a OSQPInfo object."); + + /* + * Register the mapping between struct field name and the info struct memory location + */ + // Solver status + m_structFields.push_back(new OSQPStructFieldCharArray(m_struct->status, 32, "status")); + m_structFields.push_back(new OSQPStructField(&m_struct->status_val, "status_val")); + m_structFields.push_back(new OSQPStructField(&m_struct->status_polish, "status_polish")); + + // Solution quality + m_structFields.push_back(new OSQPStructField(&m_struct->obj_val, "obj_val")); + m_structFields.push_back(new OSQPStructField(&m_struct->prim_res, "prim_res")); + m_structFields.push_back(new OSQPStructField(&m_struct->dual_res, "dual_res")); + + // Algorithm information + m_structFields.push_back(new OSQPStructField(&m_struct->iter, "iter")); + m_structFields.push_back(new OSQPStructField(&m_struct->rho_updates, "rho_updates")); + m_structFields.push_back(new OSQPStructField(&m_struct->rho_estimate, "rho_estimate")); + + // Timing information + m_structFields.push_back(new OSQPStructField(&m_struct->setup_time, "setup_time")); + m_structFields.push_back(new OSQPStructField(&m_struct->solve_time, "solve_time")); + m_structFields.push_back(new OSQPStructField(&m_struct->update_time, "update_time")); + m_structFields.push_back(new OSQPStructField(&m_struct->polish_time, "polish_time")); + m_structFields.push_back(new OSQPStructField(&m_struct->run_time, "run_time")); +} + + +// Instantiate the OSQPInfo wrapper class +template class OSQPStructWrapper; \ No newline at end of file diff --git a/c_sources/osqp_struct_settings.cpp b/c_sources/osqp_struct_settings.cpp new file mode 100644 index 0000000..92a8c4b --- /dev/null +++ b/c_sources/osqp_struct_settings.cpp @@ -0,0 +1,63 @@ +// SPDX-License-Identifier: Apache-2.0 + +#include +#include "osqp_struct.h" + +/* + * Specialization for the settings struct + */ +template<> +void OSQPStructWrapper::registerFields() { + m_struct = static_cast(c_calloc(1, sizeof(OSQPSettings))); + + if(!m_struct) + mexErrMsgTxt("Failed to allocate a OSQPSettings object."); + + osqp_set_default_settings(m_struct); + + /* + * Register the mapping between struct field name and the settings memory location + */ + m_structFields.push_back(new OSQPStructField(&m_struct->device, "device")); + m_structFields.push_back(new OSQPStructField(&m_struct->linsys_solver, "linsys_solver")); + m_structFields.push_back(new OSQPStructField(&m_struct->verbose, "verbose")); + m_structFields.push_back(new OSQPStructField(&m_struct->warm_starting, "warm_starting")); + m_structFields.push_back(new OSQPStructField(&m_struct->scaling, "scaling")); + m_structFields.push_back(new OSQPStructField(&m_struct->polishing, "polishing")); + + // ADMM parameters + m_structFields.push_back(new OSQPStructField(&m_struct->rho, "rho")); + m_structFields.push_back(new OSQPStructField(&m_struct->rho_is_vec, "rho_is_vec")); + m_structFields.push_back(new OSQPStructField(&m_struct->sigma, "sigma")); + m_structFields.push_back(new OSQPStructField(&m_struct->alpha, "alpha")); + + // CG settings + m_structFields.push_back(new OSQPStructField(&m_struct->cg_max_iter, "cg_max_iter")); + m_structFields.push_back(new OSQPStructField(&m_struct->cg_tol_reduction, "cg_tol_reduction")); + m_structFields.push_back(new OSQPStructField(&m_struct->cg_tol_fraction, "cg_tol_fraction")); + m_structFields.push_back(new OSQPStructField(&m_struct->cg_precond, "cg_precond")); + + // adaptive rho logic + m_structFields.push_back(new OSQPStructField(&m_struct->adaptive_rho, "adaptive_rho")); + m_structFields.push_back(new OSQPStructField(&m_struct->adaptive_rho_interval, "adaptive_rho_interval")); + m_structFields.push_back(new OSQPStructField(&m_struct->adaptive_rho_fraction, "adaptive_rho_fraction")); + m_structFields.push_back(new OSQPStructField(&m_struct->adaptive_rho_tolerance, "adaptive_rho_tolerance")); + + // termination parameters + m_structFields.push_back(new OSQPStructField(&m_struct->max_iter, "max_iter")); + m_structFields.push_back(new OSQPStructField(&m_struct->eps_abs, "eps_abs")); + m_structFields.push_back(new OSQPStructField(&m_struct->eps_rel, "eps_rel")); + m_structFields.push_back(new OSQPStructField(&m_struct->eps_prim_inf, "eps_prim_inf")); + m_structFields.push_back(new OSQPStructField(&m_struct->eps_dual_inf, "eps_dual_inf")); + m_structFields.push_back(new OSQPStructField(&m_struct->scaled_termination, "scaled_termination")); + m_structFields.push_back(new OSQPStructField(&m_struct->check_termination, "check_termination")); + m_structFields.push_back(new OSQPStructField(&m_struct->time_limit, "time_limit")); + + // polishing parameters + m_structFields.push_back(new OSQPStructField(&m_struct->delta, "delta")); + m_structFields.push_back(new OSQPStructField(&m_struct->polish_refine_iter, "polish_refine_iter")); +} + + +// Instantiate the OSQPSettings wrapper class +template class OSQPStructWrapper; \ No newline at end of file diff --git a/codegen/files_to_generate/CMakeLists.txt b/codegen/files_to_generate/CMakeLists.txt deleted file mode 100644 index 91035c4..0000000 --- a/codegen/files_to_generate/CMakeLists.txt +++ /dev/null @@ -1,141 +0,0 @@ -# Minimum version required -cmake_minimum_required (VERSION 3.5) - -# Project name -project (osqp) - - -# Set the output folder where your program will be created -set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/out) -set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/out) - - -# Detect operating system -# ---------------------------------------------- -message(STATUS "We are on a ${CMAKE_SYSTEM_NAME} system") -if(${CMAKE_SYSTEM_NAME} STREQUAL "Linux") - set(IS_LINUX ON) -elseif(${CMAKE_SYSTEM_NAME} STREQUAL "Darwin") - set(IS_MAC ON) -elseif(${CMAKE_SYSTEM_NAME} STREQUAL "Windows") - set(IS_WINDOWS ON) -endif() - - -# Set options -# ---------------------------------------------- - - -# Is the code generated for embedded platforms? -# 1 : Yes. Matrix update not allowed. -# 2 : Yes. Matrix update allowed. - -if (NOT DEFINED EMBEDDED) # enable embedded anyway - set (EMBEDDED EMBEDDED_FLAG) -endif() - -message(STATUS "Embedded is ${EMBEDDED}") -message(STATUS "Passing EMBEDDED flag to compiler") - -# Is printing enabled? -option (PRINTING "Enable solver printing" ON) -if (DEFINED EMBEDDED) - message(STATUS "Disabling printing for embedded") - set(PRINTING OFF) -endif() -message(STATUS "Printing is ${PRINTING}") - - -# Is profiling enabled? -option (PROFILING "Enable solver profiling (timing)" ON) -if (DEFINED EMBEDDED) - message(STATUS "Disabling profiling for embedded") - set(PROFILING OFF) -endif() -message(STATUS "Profiling is ${PROFILING}") - -# Use floats instead of integers -option (DFLOAT "Use float numbers instead of doubles" OFF) -message(STATUS "Floats are ${DFLOAT}") - -# Use long integers for indexing -option (DLONG "Use long integers (64bit) for indexing" ON) -if (NOT (CMAKE_SIZEOF_VOID_P EQUAL 8)) - message(STATUS "Disabling long integers (64bit) on 32bit machine") - set(DLONG OFF) -endif() -message(STATUS "Long integers (64bit) are ${DLONG}") - -# Types for QDLDL -# ---------------------------------------------- -if(DFLOAT) - set(QDLDL_FLOAT_TYPE "float") -else() - set(QDLDL_FLOAT_TYPE "double") -endif() - -if(DLONG) - set(QDLDL_INT_TYPE "long long") -else() - set(QDLDL_INT_TYPE "int") -endif() - -# boolean type is always unsigned char for -# now, since _Bool does not exist in C89 -set(QDLDL_BOOL_TYPE "unsigned char") - -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/configure/qdldl_types.h.in - ${CMAKE_CURRENT_SOURCE_DIR}/include/qdldl_types.h - NEWLINE_STYLE LF) - -# Set Compiler flags -# ---------------------------------------------- -set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3") -set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0 -g") -set(CMAKE_POSITION_INDEPENDENT_CODE ON) # -fPIC - -# Include math library if EMBEDDED != 1 -if(NOT (EMBEDDED EQUAL 1)) - set(CMAKE_C_STANDARD_LIBRARIES "${CMAKE_C_STANDARD_LIBRARIES} -lm") -endif() -# Include real time library in linux -if(${CMAKE_SYSTEM_NAME} STREQUAL "Linux") - set(CMAKE_C_STANDARD_LIBRARIES "${CMAKE_C_STANDARD_LIBRARIES} -lrt") -endif() - -# Generate header file with the global options -# --------------------------------------------- -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/configure/osqp_configure.h.in - ${CMAKE_CURRENT_SOURCE_DIR}/include/osqp_configure.h - NEWLINE_STYLE LF) - -# Include header directory -# ---------------------------------------------- -include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include) - -# Set sources -# ---------------------------------------------- -add_subdirectory (src/osqp) -add_subdirectory (include) - -# Append the generated workspace files and qdldl files -list (APPEND - osqp_src - ${CMAKE_CURRENT_SOURCE_DIR}/src/osqp/workspace.c - ${CMAKE_CURRENT_SOURCE_DIR}/src/osqp/qdldl.c - ${CMAKE_CURRENT_SOURCE_DIR}/src/osqp/qdldl_interface.c -) -list (APPEND - osqp_headers - ${CMAKE_CURRENT_SOURCE_DIR}/include/workspace.h - ${CMAKE_CURRENT_SOURCE_DIR}/include/qdldl.h - ${CMAKE_CURRENT_SOURCE_DIR}/include/qdldl_types.h - ${CMAKE_CURRENT_SOURCE_DIR}/include/qdldl_interface.h -) - -# Create static library for embedded solver -add_library (emosqpstatic STATIC ${osqp_src} ${osqp_headers}) - -# Create example executable -add_executable (example ${PROJECT_SOURCE_DIR}/src/example.c) -target_link_libraries (example emosqpstatic) diff --git a/codegen/files_to_generate/emosqp_mex.c b/codegen/files_to_generate/emosqp_mex.c deleted file mode 100644 index f2d6143..0000000 --- a/codegen/files_to_generate/emosqp_mex.c +++ /dev/null @@ -1,492 +0,0 @@ -#include -#include -#include "osqp.h" -#include "workspace.h" - - - -/********************************* - * Timer Structs and Functions * * - *********************************/ - -// Windows -#ifdef IS_WINDOWS - -#include - -typedef struct { - LARGE_INTEGER tic; - LARGE_INTEGER toc; - LARGE_INTEGER freq; -} PyTimer; - -// Mac -#elif defined IS_MAC - -#include - -/* Use MAC OSX mach_time for timing */ -typedef struct { - uint64_t tic; - uint64_t toc; - mach_timebase_info_data_t tinfo; -} PyTimer; - -// Linux -#else - -/* Use POSIX clocl_gettime() for timing on non-Windows machines */ -#include -#include - -typedef struct { - struct timespec tic; - struct timespec toc; -} PyTimer; - -#endif - - -/** - * Timer Methods - */ - -// Windows -#ifdef IS_WINDOWS - -void tic(PyTimer* t) -{ - QueryPerformanceFrequency(&t->freq); - QueryPerformanceCounter(&t->tic); -} - -c_float toc(PyTimer* t) -{ - QueryPerformanceCounter(&t->toc); - return ((t->toc.QuadPart - t->tic.QuadPart) / (c_float)t->freq.QuadPart); -} - -// Mac -#elif defined IS_MAC - -void tic(PyTimer* t) -{ - /* read current clock cycles */ - t->tic = mach_absolute_time(); -} - -c_float toc(PyTimer* t) -{ - uint64_t duration; /* elapsed time in clock cycles*/ - - t->toc = mach_absolute_time(); - duration = t->toc - t->tic; - - /*conversion from clock cycles to nanoseconds*/ - mach_timebase_info(&(t->tinfo)); - duration *= t->tinfo.numer; - duration /= t->tinfo.denom; - - return (c_float)duration / 1e9; -} - - -// Linux -#else - -/* read current time */ -void tic(PyTimer* t) -{ - clock_gettime(CLOCK_MONOTONIC, &t->tic); -} - - -/* return time passed since last call to tic on this timer */ -c_float toc(PyTimer* t) -{ - struct timespec temp; - - clock_gettime(CLOCK_MONOTONIC, &t->toc); - - if ((t->toc.tv_nsec - t->tic.tv_nsec)<0) { - temp.tv_sec = t->toc.tv_sec - t->tic.tv_sec-1; - temp.tv_nsec = 1e9+t->toc.tv_nsec - t->tic.tv_nsec; - } else { - temp.tv_sec = t->toc.tv_sec - t->tic.tv_sec; - temp.tv_nsec = t->toc.tv_nsec - t->tic.tv_nsec; - } - return (c_float)temp.tv_sec + (c_float)temp.tv_nsec / 1e9; -} - - -#endif - -/**************************************** - * END( Timer Structs and Functions ) * * - ****************************************/ - - - -// Internal utility functions -c_float* copyToCfloatVector(double * vecData, c_int numel); -void castToDoubleArr(c_float *arr, double* arr_out, c_int len); -void setToNaN(double* arr_out, c_int len); -#if EMBEDDED != 1 -c_int* copyDoubleToCintVector(double* vecData, c_int numel); -void change1To0Indexing(c_int *vecData, c_int numel); -#endif - - -// Function handler -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - // Get the command string - char cmd[64]; - if (nrhs < 1 || mxGetString(prhs[0], cmd, sizeof(cmd))) - mexErrMsgTxt("First input should be a command string less than 64 characters long."); - - - // SOLVE - if (!strcmp("solve", cmd)) { - - // Allocate timer - double solve_time; - PyTimer * timer; - timer = mxMalloc(sizeof(PyTimer)); - - if (nlhs != 5 || nrhs != 1){ - mexErrMsgTxt("Solve : wrong number of inputs / outputs"); - } - - // solve the problem - tic(timer); // start timer - osqp_solve(&workspace); - solve_time = toc(timer); // stop timer - - - // Allocate space for solution - // primal variables - plhs[0] = mxCreateDoubleMatrix((&workspace)->data->n, 1, mxREAL); - // dual variables - plhs[1] = mxCreateDoubleMatrix((&workspace)->data->m, 1, mxREAL); - // status value - plhs[2] = mxCreateDoubleScalar((&workspace)->info->status_val); - // number of iterations - plhs[3] = mxCreateDoubleScalar((&workspace)->info->iter); - // solve time - plhs[4] = mxCreateDoubleScalar(solve_time); - - - //copy results to mxArray outputs - //assume that three outputs will always - //be returned to matlab-side class wrapper - if (((&workspace)->info->status_val != OSQP_PRIMAL_INFEASIBLE) && - ((&workspace)->info->status_val != OSQP_PRIMAL_INFEASIBLE_INACCURATE) && - ((&workspace)->info->status_val != OSQP_DUAL_INFEASIBLE) && - ((&workspace)->info->status_val != OSQP_DUAL_INFEASIBLE_INACCURATE)){ - - //primal variables - castToDoubleArr((&workspace)->solution->x, mxGetPr(plhs[0]), (&workspace)->data->n); - - //dual variables - castToDoubleArr((&workspace)->solution->y, mxGetPr(plhs[1]), (&workspace)->data->m); - - } else { // Problem is primal or dual infeasible -> NaN values - - // Set primal and dual variables to NaN - setToNaN(mxGetPr(plhs[0]), (&workspace)->data->n); - setToNaN(mxGetPr(plhs[1]), (&workspace)->data->m); - } - - return; - } - - - // update linear cost - if (!strcmp("update_lin_cost", cmd)) { - - // Fill q - const mxArray *q = prhs[1]; - - // Copy vector to ensure it is cast as c_float - c_float *q_vec; - if(!mxIsEmpty(q)){ - q_vec = copyToCfloatVector(mxGetPr(q), (&workspace)->data->n); - } - - if(!mxIsEmpty(q)){ - osqp_update_lin_cost(&workspace, q_vec); - } - - // Free - if(!mxIsEmpty(q)) mxFree(q_vec); - - return; - } - - - // update lower bound - if (!strcmp("update_lower_bound", cmd)) { - - // Fill l - const mxArray *l = prhs[1]; - - // Copy vector to ensure it is cast as c_float - c_float *l_vec; - if(!mxIsEmpty(l)){ - l_vec = copyToCfloatVector(mxGetPr(l), (&workspace)->data->m); - } - - if(!mxIsEmpty(l)){ - osqp_update_lower_bound(&workspace, l_vec); - } - - // Free - if(!mxIsEmpty(l)) mxFree(l_vec); - - return; - } - - - // update upper bound - if (!strcmp("update_upper_bound", cmd)) { - - // Fill l - const mxArray *u = prhs[1]; - - // Copy vector to ensure it is cast as c_float - c_float *u_vec; - if(!mxIsEmpty(u)){ - u_vec = copyToCfloatVector(mxGetPr(u), (&workspace)->data->m); - } - - if(!mxIsEmpty(u)){ - osqp_update_upper_bound(&workspace, u_vec); - } - - // Free - if(!mxIsEmpty(u)) mxFree(u_vec); - - return; - } - - - // update bounds - if (!strcmp("update_bounds", cmd)) { - - // Fill l, u - const mxArray *l = prhs[1]; - const mxArray *u = prhs[2]; - - // Copy vectors to ensure they are cast as c_float - c_float *l_vec; - c_float *u_vec; - if(!mxIsEmpty(l)){ - l_vec = copyToCfloatVector(mxGetPr(l), (&workspace)->data->m); - } - if(!mxIsEmpty(u)){ - u_vec = copyToCfloatVector(mxGetPr(u), (&workspace)->data->m); - } - - if(!mxIsEmpty(u)){ - osqp_update_bounds(&workspace, l_vec, u_vec); - } - - // Free - if(!mxIsEmpty(l)) mxFree(l_vec); - if(!mxIsEmpty(u)) mxFree(u_vec); - - return; - } - - #if EMBEDDED != 1 - // update matrix P - if (!strcmp("update_P", cmd)) { - - // Fill Px and Px_idx - const mxArray *Px = prhs[1]; - const mxArray *Px_idx = prhs[2]; - - int Px_n = mxGetScalar(prhs[3]); - if(Px_n == 0){ - Px_n = (&workspace)->data->P->nzmax; - } - - // Copy vectors to ensure they are cast as c_float and c_int - c_float *Px_vec; - c_int *Px_idx_vec = NULL; - if(!mxIsEmpty(Px)){ - Px_vec = copyToCfloatVector(mxGetPr(Px), Px_n); - } - if(!mxIsEmpty(Px_idx)){ - Px_idx_vec = copyDoubleToCintVector(mxGetPr(Px_idx), Px_n); - - // Change indexing to match C 0-based one - change1To0Indexing(Px_idx_vec, Px_n); - - } - - if(!mxIsEmpty(Px)){ - osqp_update_P((&workspace), Px_vec, Px_idx_vec, Px_n); - } - - // Free - if(!mxIsEmpty(Px)) mxFree(Px_vec); - if(!mxIsEmpty(Px_idx)) mxFree(Px_idx_vec); - - return; - } - - // update matrix A - if (!strcmp("update_A", cmd)) { - - // Fill Ax and Ax_idx - const mxArray *Ax = prhs[1]; - const mxArray *Ax_idx = prhs[2]; - - int Ax_n = mxGetScalar(prhs[3]); - if(Ax_n == 0){ - Ax_n = (&workspace)->data->A->nzmax; - } - - // Copy vectors to ensure they are cast as c_float and c_int - c_float *Ax_vec; - c_int *Ax_idx_vec = NULL; - if(!mxIsEmpty(Ax)){ - Ax_vec = copyToCfloatVector(mxGetPr(Ax), Ax_n); - } - if(!mxIsEmpty(Ax_idx)){ - Ax_idx_vec = copyDoubleToCintVector(mxGetPr(Ax_idx), Ax_n); - - // Change indexing to match C 0-based one - change1To0Indexing(Ax_idx_vec, Ax_n); - } - - if(!mxIsEmpty(Ax)){ - osqp_update_A((&workspace), Ax_vec, Ax_idx_vec, Ax_n); - } - - // Free - if(!mxIsEmpty(Ax)) mxFree(Ax_vec); - if(!mxIsEmpty(Ax_idx)) mxFree(Ax_idx_vec); - - return; - } - - // update matrices P and A - if (!strcmp("update_P_A", cmd)) { - - // Fill vectors - const mxArray *Px = prhs[1]; - const mxArray *Px_idx = prhs[2]; - const mxArray *Ax = prhs[4]; - const mxArray *Ax_idx = prhs[5]; - - int Px_n = mxGetScalar(prhs[3]); - int Ax_n = mxGetScalar(prhs[6]); - - if(Px_n == 0){ - Px_n = (&workspace)->data->P->nzmax; - } - if(Ax_n == 0){ - Ax_n = (&workspace)->data->A->nzmax; - } - - // Copy vectors to ensure they are cast as c_float and c_int - c_float *Px_vec, *Ax_vec; - c_int *Px_idx_vec = NULL; - c_int *Ax_idx_vec = NULL; - if(!mxIsEmpty(Px)){ - Px_vec = copyToCfloatVector(mxGetPr(Px), Px_n); - } - if(!mxIsEmpty(Px_idx)){ - Px_idx_vec = copyDoubleToCintVector(mxGetPr(Px_idx), Px_n); - - // Change indexing to match C 0-based one - change1To0Indexing(Px_idx_vec, Px_n); - } - if(!mxIsEmpty(Ax)){ - Ax_vec = copyToCfloatVector(mxGetPr(Ax), Ax_n); - } - if(!mxIsEmpty(Ax_idx)){ - Ax_idx_vec = copyDoubleToCintVector(mxGetPr(Ax_idx), Ax_n); - - // Change indexing to match C 0-based one - change1To0Indexing(Ax_idx_vec, Ax_n); - } - - if(!mxIsEmpty(Ax) && !mxIsEmpty(Px)){ - osqp_update_P_A((&workspace), Px_vec, Px_idx_vec, Px_n, - Ax_vec, Ax_idx_vec, Ax_n); - } - - // Free - if(!mxIsEmpty(Px)) mxFree(Px_vec); - if(!mxIsEmpty(Px_idx)) mxFree(Px_idx_vec); - if(!mxIsEmpty(Ax)) mxFree(Ax_vec); - if(!mxIsEmpty(Ax_idx)) mxFree(Ax_idx_vec); - - return; - } - #endif // end EMBEDDED - - // Got here, so command not recognized - mexErrMsgTxt("Command not recognized."); -} - -c_float* copyToCfloatVector(double* vecData, c_int numel){ - // This memory needs to be freed! - - c_float* out; - c_int i; - - out = mxMalloc(numel * sizeof(c_float)); - - //copy data - for(i=0; i < numel; i++){ - out[i] = (c_float)vecData[i]; - } - return out; - -} - -void castToDoubleArr(c_float *arr, double* arr_out, c_int len){ - c_int i; - for (i = 0; i < len; i++) { - arr_out[i] = (double)arr[i]; - } -} - -void setToNaN(double* arr_out, c_int len){ - c_int i; - for (i = 0; i < len; i++) { - arr_out[i] = mxGetNaN(); - } -} - -#if EMBEDDED != 1 -c_int* copyDoubleToCintVector(double* vecData, c_int numel){ - // This memory needs to be freed! - - c_int* out; - c_int i; - - out = mxMalloc(numel * sizeof(c_int)); - - //copy data - for(i=0; i < numel; i++){ - out[i] = (c_int)vecData[i]; - } - return out; - -} - -void change1To0Indexing(c_int *vecData, c_int numel){ -c_int i; // Indexing -for(i=0; i < numel; i++){ -vecData[i] -= 1; // Decrease index by 1 -} -} - - -#endif // end EMBEDDED diff --git a/codegen/files_to_generate/example.c b/codegen/files_to_generate/example.c deleted file mode 100644 index beac16e..0000000 --- a/codegen/files_to_generate/example.c +++ /dev/null @@ -1,20 +0,0 @@ -#include "stdio.h" -#include - -#include "workspace.h" -#include "osqp.h" - -int main(int argc, char **argv) { - - // Solve Problem - osqp_solve(&workspace); - - // Print status - printf("Status: %s\n", (&workspace)->info->status); - printf("Number of iterations: %d\n", (int)((&workspace)->info->iter)); - printf("Objective value: %.4e\n", (&workspace)->info->obj_val); - printf("Primal residual: %.4e\n", (&workspace)->info->pri_res); - printf("Dual residual: %.4e\n", (&workspace)->info->dua_res); - - return 0; -} diff --git a/codegen/make_emosqp.m b/codegen/make_emosqp.m deleted file mode 100644 index 5bb87ee..0000000 --- a/codegen/make_emosqp.m +++ /dev/null @@ -1,78 +0,0 @@ -function make_emosqp(target_dir, mex_cfile, EMBEDDED_FLAG, FLOAT_FLAG, LONG_FLAG) -% Matlab MEX makefile for code generated solver. - - -% Get make and mex commands -mex_cmd = sprintf('mex -O -silent'); - -% Add arguments to mex compiler -mexoptflags = '-DMATLAB'; - -% If running on linux, include the c99 option so GCC uses c99 to compile. -% Otherwise it will throw errors about the comment style -if ( ~ismac() && isunix() ) - mexoptflags = sprintf('%s CFLAGS="$CFLAGS -std=c99"', mexoptflags); -end - -% Add embedded flag -cmake_args = sprintf('-DEMBEDDED:INT=%i', EMBEDDED_FLAG); - -% Add float flag -cmake_args = sprintf('%s -DDFLOAT:BOOL=%s', cmake_args, FLOAT_FLAG); - -% Add long flag -cmake_args = sprintf('%s -DDLONG:BOOL=%s', cmake_args, LONG_FLAG); - - -% Generate osqp_configure.h file by running cmake -current_dir = pwd; -build_dir = fullfile(target_dir, 'build'); -cd(target_dir); -if exist(build_dir, 'dir') - rmdir('build', 's'); -end -mkdir('build'); -cd('build'); - - -% Add specific generators for windows linux or mac -if (ispc) - [status, output] = system(sprintf('%s %s -G "MinGW Makefiles" ..', 'cmake', cmake_args)); -else - [status, output] = system(sprintf('%s %s -G "Unix Makefiles" ..', 'cmake', cmake_args)); -end -if(status) - fprintf('\n'); - disp(output); - error('Error generating osqp_configure.h'); -end -cd(current_dir); - -% Set optimizer flag -if (~ispc) - mexoptflags = sprintf('%s %s', mexoptflags, 'COPTIMFLAGS=''-O3'''); -end - -% Include directory -inc_dir = fullfile(sprintf(' -I%s', target_dir), 'include'); - -% Source files -cfiles = ''; -src_files = dir(fullfile(target_dir, 'src', 'osqp', '*c')); -for i = 1 : length(src_files) - cfiles = sprintf('%s %s', cfiles, ... - fullfile(target_dir, 'src', 'osqp', src_files(i).name)); -end - -% Compile interface -fprintf('Compiling and linking osqpmex...'); - -% Compile command -cmd = sprintf('%s %s %s "%s" %s', mex_cmd, mexoptflags, inc_dir, mex_cfile, cfiles); - -% Compile -eval(cmd); -fprintf('\t\t\t\t[done]\n'); - - -end diff --git a/codegen/render_workspace.m b/codegen/render_workspace.m deleted file mode 100644 index e9d9a79..0000000 --- a/codegen/render_workspace.m +++ /dev/null @@ -1,442 +0,0 @@ -function render_workspace( work, hfname, cfname, embedded_flag ) -%RENDER_WORKSPACE Write workspace to header file. - -incFile = fopen(hfname, 'w'); -srcFile = fopen(cfname, 'w'); - -% Add an include-guard statement -[~, fname, ~] = fileparts(hfname); -incGuard = [upper(fname), '_H']; -fprintf(incFile, '#ifndef %s\n', incGuard); -fprintf(incFile, '#define %s\n\n', incGuard); - -% Print comment headers containing the generation time into the files -daystr = datestr(now, 'mmmm dd, yyyy'); -timestr = datestr(now, 'HH:MM:SS'); -fprintf(incFile, '/*\n'); -fprintf(incFile, ' * This file was autogenerated by OSQP-Matlab on %s at %s.\n', daystr, timestr); -fprintf(incFile, ' * \n'); -fprintf(incFile, ' * This file contains the prototypes for all the workspace variables needed\n'); -fprintf(incFile, ' * by OSQP. The actual data is contained inside workspace.c.\n'); -fprintf(incFile, ' */\n\n'); - -fprintf(srcFile, '/*\n'); -fprintf(srcFile, ' * This file was autogenerated by OSQP-Matlab on %s at %s.\n', daystr, timestr); -fprintf(srcFile, ' * \n'); -fprintf(srcFile, ' * This file contains the workspace variables needed by OSQP.\n'); -fprintf(srcFile, ' */\n\n'); - - -% Include types, constants and private header -fprintf(incFile, '#include \"types.h\"\n'); -fprintf(incFile, '#include \"qdldl_interface.h\"\n\n'); - -fprintf(srcFile, '#include \"types.h\"\n'); -fprintf(srcFile, '#include \"qdldl_interface.h\"\n\n'); - -% Write data structure -write_data_src(srcFile, work.data); -write_data_inc(incFile, work.data); - -% Write settings structure -write_settings_src(srcFile, work.settings, embedded_flag); -write_settings_inc(incFile, work.settings, embedded_flag); - -% Write scaling structure -write_scaling_src(srcFile, work.scaling); -write_scaling_inc(incFile, work.scaling); - -% Write linsys_solver structure -write_linsys_solver_src(srcFile, work.linsys_solver, embedded_flag); -write_linsys_solver_inc(incFile, work.linsys_solver, embedded_flag); - -% Define empty solution structure -write_solution_src(srcFile, work.data.n, work.data.m); -write_solution_inc(incFile, work.data.n, work.data.m); - -% Define info structure -write_info_src(srcFile); -write_info_inc(incFile); - -% Define workspace structure -write_workspace_src(srcFile, work.data.n, work.data.m, work.rho_vectors, embedded_flag); -write_workspace_inc(incFile, work.data.n, work.data.m, work.rho_vectors, embedded_flag); - -% The endif for the include-guard -fprintf(incFile, '#endif // ifndef %s\n', incGuard); - -fclose(incFile); -fclose(srcFile); - -end - - - -function write_data_src( f, data ) -%WRITE_DATA_SRC Write data structure to file. - -fprintf(f, '// Define data structure\n'); - -% Define matrix P -write_mat(f, data.P, 'Pdata'); - -% Define matrix A -write_mat(f, data.A, 'Adata'); - -% Define other data vectors -write_vec(f, data.q, 'qdata', 'c_float'); -write_vec(f, data.l, 'ldata', 'c_float'); -write_vec(f, data.u, 'udata', 'c_float'); - -% Define data structure -fprintf(f, 'OSQPData data = {'); -fprintf(f, '%d, ', data.n); -fprintf(f, '%d, ', data.m); -fprintf(f, '&Pdata, &Adata, qdata, ldata, udata'); -fprintf(f, '};\n\n'); - -end - -function write_data_inc( f, data ) -%WRITE_DATA_INC Write data structure prototypes to file. - -fprintf(f, '// Data structure prototypes\n'); - -% Define matrix P -write_mat_extern(f, data.P, 'Pdata'); - -% Define matrix A -write_mat_extern(f, data.A, 'Adata'); - -% Define other data vectors -write_vec_extern(f, data.q, 'qdata', 'c_float'); -write_vec_extern(f, data.l, 'ldata', 'c_float'); -write_vec_extern(f, data.u, 'udata', 'c_float'); - -% Define data structure -fprintf(f, 'extern OSQPData data;\n\n'); - -end - - -function write_settings_src( f, settings, embedded_flag ) -%WRITE_SETTINGS_SRC Write settings structure to file. - -fprintf(f, '// Define settings structure\n'); -fprintf(f, 'OSQPSettings settings = {'); -fprintf(f, '(c_float)%.20f, ', settings.rho); -fprintf(f, '(c_float)%.20f, ', settings.sigma); -fprintf(f, '%d, ', settings.scaling); - -if embedded_flag ~= 1 - fprintf(f, '%d, ', settings.adaptive_rho); - fprintf(f, '%d, ', settings.adaptive_rho_interval); - fprintf(f, '(c_float)%.20f,', settings.adaptive_rho_tolerance); -end - -fprintf(f, '%d, ', settings.max_iter); -fprintf(f, '(c_float)%.20f, ', settings.eps_abs); -fprintf(f, '(c_float)%.20f, ', settings.eps_rel); -fprintf(f, '(c_float)%.20f, ', settings.eps_prim_inf); -fprintf(f, '(c_float)%.20f, ', settings.eps_dual_inf); -fprintf(f, '(c_float)%.20f, ', settings.alpha); -fprintf(f, '(enum linsys_solver_type) LINSYS_SOLVER, '); - - -fprintf(f, '%d, ', settings.scaled_termination); -fprintf(f, '%d, ', settings.check_termination); -fprintf(f, '%d, ', settings.warm_start); - -fprintf(f, '};\n\n'); - -end - -function write_settings_inc( f, settings, embedded_flag ) -%WRITE_SETTINGS_INC Write prototype for settings structure to file. - -fprintf(f, '// Settings structure prototype\n'); -fprintf(f, 'extern OSQPSettings settings;\n\n'); - -end - - -function write_scaling_src( f, scaling ) -%WRITE_SCALING_SRC Write scaling structure to file. - -fprintf(f, '// Define scaling structure\n'); - -if ~isempty(scaling) - write_vec(f, scaling.D, 'Dscaling', 'c_float'); - write_vec(f, scaling.Dinv, 'Dinvscaling', 'c_float'); - write_vec(f, scaling.E, 'Escaling', 'c_float'); - write_vec(f, scaling.Einv, 'Einvscaling', 'c_float'); - fprintf(f, 'OSQPScaling scaling = {'); - fprintf(f, '(c_float)%.20f, ', scaling.c); - fprintf(f, 'Dscaling, Escaling, '); - fprintf(f, '(c_float)%.20f, ', scaling.cinv); - fprintf(f, 'Dinvscaling, Einvscaling};\n\n'); -else - fprintf(f, 'OSQPScaling scaling;\n\n'); -end - -end - - -function write_scaling_inc( f, scaling ) -%WRITE_SCALING_INC Write prototypes for the scaling structure to file. - -fprintf(f, '// Scaling structure prototypes\n'); - -if ~isempty(scaling) - write_vec_extern(f, scaling.D, 'Dscaling', 'c_float'); - write_vec_extern(f, scaling.Dinv, 'Dinvscaling', 'c_float'); - write_vec_extern(f, scaling.E, 'Escaling', 'c_float'); - write_vec_extern(f, scaling.Einv, 'Einvscaling', 'c_float'); -end - -fprintf(f, 'extern OSQPScaling scaling;\n\n'); - - -end - -function write_linsys_solver_src( f, linsys_solver, embedded_flag ) -%WRITE_LINSYS_SOLVER_SRC Write linsys_solver structure to file. - -fprintf(f, '// Define linsys_solver structure\n'); -write_mat(f, linsys_solver.L, 'linsys_solver_L') -write_vec(f, linsys_solver.Dinv, 'linsys_solver_Dinv', 'c_float') -write_vec(f, linsys_solver.P, 'linsys_solver_P', 'c_int') -fprintf(f, 'c_float linsys_solver_bp[%d];\n', length(linsys_solver.bp)); -fprintf(f, 'c_float linsys_solver_sol[%d];\n', length(linsys_solver.sol)); -write_vec(f, linsys_solver.rho_inv_vec, 'linsys_solver_rho_inv_vec', 'c_float') - -if embedded_flag ~= 1 - write_vec(f, linsys_solver.Pdiag_idx, 'linsys_solver_Pdiag_idx', 'c_int'); - write_mat(f, linsys_solver.KKT, 'linsys_solver_KKT'); - write_vec(f, linsys_solver.PtoKKT, 'linsys_solver_PtoKKT', 'c_int'); - write_vec(f, linsys_solver.AtoKKT, 'linsys_solver_AtoKKT', 'c_int'); - write_vec(f, linsys_solver.rhotoKKT, 'linsys_solver_rhotoKKT', 'c_int'); - write_vec(f, linsys_solver.D, 'linsys_solver_D', 'QDLDL_float'); - write_vec(f, linsys_solver.etree, 'linsys_solver_etree', 'QDLDL_int'); - write_vec(f, linsys_solver.Lnz, 'linsys_solver_Lnz', 'QDLDL_int'); - fprintf(f, 'QDLDL_int linsys_solver_iwork[%d];\n', length(linsys_solver.iwork)); - fprintf(f, 'QDLDL_bool linsys_solver_bwork[%d];\n', length(linsys_solver.bwork)); - fprintf(f, 'QDLDL_float linsys_solver_fwork[%d];\n', length(linsys_solver.fwork)); -end - -fprintf(f, 'qdldl_solver linsys_solver = '); -fprintf(f, '{QDLDL_SOLVER, &solve_linsys_qdldl, '); - -if embedded_flag ~= 1 - fprintf(f, '&update_linsys_solver_matrices_qdldl, &update_linsys_solver_rho_vec_qdldl, '); -end - -fprintf(f, '&linsys_solver_L, linsys_solver_Dinv, linsys_solver_P, linsys_solver_bp, linsys_solver_sol, linsys_solver_rho_inv_vec, '); -fprintf(f, '(c_float)%.20f, ', linsys_solver.sigma); -fprintf(f, '%d, ', linsys_solver.n); -fprintf(f, '%d, ', linsys_solver.m); - -if embedded_flag ~= 1 - fprintf(f, 'linsys_solver_Pdiag_idx, '); - fprintf(f, '%d, ', linsys_solver.Pdiag_n); - fprintf(f, ['&linsys_solver_KKT, linsys_solver_PtoKKT, linsys_solver_AtoKKT, linsys_solver_rhotoKKT, ', ... - 'linsys_solver_D, linsys_solver_etree, linsys_solver_Lnz, linsys_solver_iwork, linsys_solver_bwork, linsys_solver_fwork, ']); -end - -fprintf(f, '};\n\n'); - -end - - -function write_linsys_solver_inc( f, linsys_solver, embedded_flag ) -%WRITE_LINSYS_SOLVER_INC Write prototypes for linsys_solver structure to file. - -fprintf(f, '// Prototypes for linsys_solver structure\n'); -write_mat_extern(f, linsys_solver.L, 'linsys_solver_L') -write_vec_extern(f, linsys_solver.Dinv, 'linsys_solver_Dinv', 'c_float') -write_vec_extern(f, linsys_solver.P, 'linsys_solver_P', 'c_int') -fprintf(f, 'extern c_float linsys_solver_bp[%d];\n', length(linsys_solver.bp)); -fprintf(f, 'extern c_float linsys_solver_sol[%d];\n', length(linsys_solver.sol)); -write_vec_extern(f, linsys_solver.rho_inv_vec, 'linsys_solver_rho_inv_vec', 'c_float') - -if embedded_flag ~= 1 - write_vec_extern(f, linsys_solver.Pdiag_idx, 'linsys_solver_Pdiag_idx', 'c_int'); - write_mat_extern(f, linsys_solver.KKT, 'linsys_solver_KKT'); - write_vec_extern(f, linsys_solver.PtoKKT, 'linsys_solver_PtoKKT', 'c_int'); - write_vec_extern(f, linsys_solver.AtoKKT, 'linsys_solver_AtoKKT', 'c_int'); - write_vec_extern(f, linsys_solver.rhotoKKT, 'linsys_solver_rhotoKKT', 'c_int'); - write_vec_extern(f, linsys_solver.D, 'linsys_solver_D', 'QDLDL_float'); - write_vec_extern(f, linsys_solver.etree, 'linsys_solver_etree', 'QDLDL_int'); - write_vec_extern(f, linsys_solver.Lnz, 'linsys_solver_Lnz', 'QDLDL_int'); - fprintf(f, 'extern QDLDL_int linsys_solver_iwork[%d];\n', length(linsys_solver.iwork)); - fprintf(f, 'extern QDLDL_bool linsys_solver_bwork[%d];\n', length(linsys_solver.bwork)); - fprintf(f, 'extern QDLDL_float linsys_solver_fwork[%d];\n', length(linsys_solver.fwork)); -end - -fprintf(f, 'extern qdldl_solver linsys_solver;\n\n'); - -end - - -function write_solution_src( f, n, m ) -%WRITE_SOLUTION_SRC Preallocate solution vectors - -fprintf(f, '// Define solution\n'); -fprintf(f, 'c_float xsolution[%d];\n', n); -fprintf(f, 'c_float ysolution[%d];\n\n', m); -fprintf(f, 'OSQPSolution solution = {xsolution, ysolution};\n\n'); - -end - -function write_solution_inc( f, n, m ) -%WRITE_SOLUTION_INC Prototypes for solution vectors - -fprintf(f, '// Prototypes for solution\n'); -fprintf(f, 'extern c_float xsolution[%d];\n', n); -fprintf(f, 'extern c_float ysolution[%d];\n\n', m); -fprintf(f, 'extern OSQPSolution solution;\n\n'); - -end - - -function write_info_src( f ) -%WRITE_INFO_SRC Preallocate info structure - -fprintf(f, '// Define info\n'); -fprintf(f, 'OSQPInfo info = {0, "Unsolved", OSQP_UNSOLVED, (c_float)0.0, (c_float)0.0, (c_float)0.0};\n\n'); - -end - -function write_info_inc( f ) -%WRITE_INFO_INC Prototype for info structure - -fprintf(f, '// Prototype for info structure\n'); -fprintf(f, 'extern OSQPInfo info;\n\n'); - -end - - -function write_workspace_src( f, n, m, rho_vectors, embedded_flag ) -%WRITE_WORKSPACE_SRC Preallocate workspace structure and populate rho_vectors - -fprintf(f, '// Define workspace\n'); -write_vec(f, rho_vectors.rho_vec, 'work_rho_vec', 'c_float'); -write_vec(f, rho_vectors.rho_inv_vec, 'work_rho_inv_vec', 'c_float'); -if embedded_flag ~= 1 - write_vec(f, rho_vectors.constr_type, 'work_constr_type', 'c_int'); -end -fprintf(f, 'c_float work_x[%d];\n', n); -fprintf(f, 'c_float work_y[%d];\n', m); -fprintf(f, 'c_float work_z[%d];\n', m); -fprintf(f, 'c_float work_xz_tilde[%d];\n', n+m); -fprintf(f, 'c_float work_x_prev[%d];\n', n); -fprintf(f, 'c_float work_z_prev[%d];\n', m); -fprintf(f, 'c_float work_Ax[%d];\n', m); -fprintf(f, 'c_float work_Px[%d];\n', n); -fprintf(f, 'c_float work_Aty[%d];\n', n); -fprintf(f, 'c_float work_delta_y[%d];\n', m); -fprintf(f, 'c_float work_Atdelta_y[%d];\n', n); -fprintf(f, 'c_float work_delta_x[%d];\n', n); -fprintf(f, 'c_float work_Pdelta_x[%d];\n', n); -fprintf(f, 'c_float work_Adelta_x[%d];\n', m); -fprintf(f, 'c_float work_D_temp[%d];\n', n); -fprintf(f, 'c_float work_D_temp_A[%d];\n', n); -fprintf(f, 'c_float work_E_temp[%d];\n\n', m); - -fprintf(f, 'OSQPWorkspace workspace = {\n'); -fprintf(f, '&data, (LinSysSolver *)&linsys_solver,\n'); -fprintf(f, 'work_rho_vec, work_rho_inv_vec,\n'); -if embedded_flag ~= 1 - fprintf(f, 'work_constr_type,\n'); -end -fprintf(f, 'work_x, work_y, work_z, work_xz_tilde,\n'); -fprintf(f, 'work_x_prev, work_z_prev,\n'); -fprintf(f, 'work_Ax, work_Px, work_Aty,\n'); -fprintf(f, 'work_delta_y, work_Atdelta_y,\n'); -fprintf(f, 'work_delta_x, work_Pdelta_x, work_Adelta_x,\n'); -fprintf(f, 'work_D_temp, work_D_temp_A, work_E_temp,\n'); -fprintf(f, '&settings, &scaling, &solution, &info};\n\n'); - -end - -function write_workspace_inc( f, n, m, rho_vectors, embedded_flag ) -%WRITE_WORKSPACE_INC Prototypes for the workspace structure and rho_vectors - -fprintf(f, '// Prototypes for the workspace\n'); -write_vec_extern(f, rho_vectors.rho_vec, 'work_rho_vec', 'c_float'); -write_vec_extern(f, rho_vectors.rho_inv_vec, 'work_rho_inv_vec', 'c_float'); -if embedded_flag ~= 1 - write_vec_extern(f, rho_vectors.constr_type, 'work_constr_type', 'c_int'); -end -fprintf(f, 'extern c_float work_x[%d];\n', n); -fprintf(f, 'extern c_float work_y[%d];\n', m); -fprintf(f, 'extern c_float work_z[%d];\n', m); -fprintf(f, 'extern c_float work_xz_tilde[%d];\n', n+m); -fprintf(f, 'extern c_float work_x_prev[%d];\n', n); -fprintf(f, 'extern c_float work_z_prev[%d];\n', m); -fprintf(f, 'extern c_float work_Ax[%d];\n', m); -fprintf(f, 'extern c_float work_Px[%d];\n', n); -fprintf(f, 'extern c_float work_Aty[%d];\n', n); -fprintf(f, 'extern c_float work_delta_y[%d];\n', m); -fprintf(f, 'extern c_float work_Atdelta_y[%d];\n', n); -fprintf(f, 'extern c_float work_delta_x[%d];\n', n); -fprintf(f, 'extern c_float work_Pdelta_x[%d];\n', n); -fprintf(f, 'extern c_float work_Adelta_x[%d];\n', m); -fprintf(f, 'extern c_float work_D_temp[%d];\n', n); -fprintf(f, 'extern c_float work_D_temp_A[%d];\n', n); -fprintf(f, 'extern c_float work_E_temp[%d];\n\n', m); - -fprintf(f, 'extern OSQPWorkspace workspace;\n\n'); - -end - - -function write_vec(f, vec, name, vec_type) -%WRITE_VEC Write vector to file. - -fprintf(f, '%s %s[%d] = {\n', vec_type, name, length(vec)); - -% Write vector elements -for i = 1 : length(vec) - if strcmp(vec_type, 'c_float') - fprintf(f, '(c_float)%.20f,\n', vec(i)); - else - fprintf(f, '%i,\n', vec(i)); - end -end -fprintf(f, '};\n'); - -end - -function write_vec_extern(f, vec, name, vec_type) -%WRITE_VEC_EXTERN Write vector prototype to file. - -fprintf(f, 'extern %s %s[%d];\n', vec_type, name, length(vec)); - -end - - -function write_mat(f, mat, name) -%WRITE_MAT Write matrix in CSC format to file. - -write_vec(f, mat.i, [name, '_i'], 'c_int'); -write_vec(f, mat.p, [name, '_p'], 'c_int'); -write_vec(f, mat.x, [name, '_x'], 'c_float'); - -fprintf(f, 'csc %s = {', name); -fprintf(f, '%d, ', mat.nzmax); -fprintf(f, '%d, ', mat.m); -fprintf(f, '%d, ', mat.n); -fprintf(f, '%s_p, ', name); -fprintf(f, '%s_i, ', name); -fprintf(f, '%s_x, ', name); -fprintf(f, '%d};\n', mat.nz); - -end - -function write_mat_extern(f, mat, name) -%WRITE_MAT_EXTERN Write matrix the prototype for the matrix. - -fprintf(f, 'extern csc %s;\n', name); - -end diff --git a/examples/easy_lasso_qp.m b/examples/easy_lasso_qp.m index fe13dac..44523d8 100644 --- a/examples/easy_lasso_qp.m +++ b/examples/easy_lasso_qp.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + %% Lasso problem % Problem dimensions and sparity diff --git a/examples/examples_old.m b/examples/examples_old.m index a584c42..1fa160f 100644 --- a/examples/examples_old.m +++ b/examples/examples_old.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + %% OSQPtest diff --git a/examples/examples_osqpmatlab.m b/examples/examples_osqpmatlab.m index 15b4c77..ae2e246 100644 --- a/examples/examples_osqpmatlab.m +++ b/examples/examples_osqpmatlab.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + %% Simple problem m = 50; n = 100; diff --git a/examples/interrupt.m b/examples/interrupt.m index bb81ee8..93e737c 100644 --- a/examples/interrupt.m +++ b/examples/interrupt.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + %% Big lasso problem % Problem dimensions and sparity @@ -27,7 +29,7 @@ solver = osqp; solver.setup(P, q, A, l, u, ... - 'check_termination', 0, 'polish', 0, 'max_iter', 50000); + 'check_termination', 0, 'polishing', 0, 'max_iter', 50000); res = solver.solve(); % PRESS CTRL-C BEFORE REACHING 50K ITERATIONS diff --git a/examples/codegen_test.m b/examples/osqp_codegen.m similarity index 79% rename from examples/codegen_test.m rename to examples/osqp_codegen.m index 5734246..a07acdf 100644 --- a/examples/codegen_test.m +++ b/examples/osqp_codegen.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + %% Simple problem m = 50; n = 100; @@ -11,4 +13,4 @@ m = osqp; m.setup(problem.P, problem.q, problem.A, problem.l, problem.u); -m.codegen('code', 'parameters', 'matrices'); \ No newline at end of file +m.codegen('out/', 'emosqp'); \ No newline at end of file diff --git a/examples/osqp_demo.m b/examples/osqp_demo.m index da9d2eb..d0d9b0e 100644 --- a/examples/osqp_demo.m +++ b/examples/osqp_demo.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + % Demo showing the usage of OSQP from Matlab and the code generation features. % This problem is the same one that is presented in the osqp_demo.c file. @@ -13,10 +15,4 @@ solver.setup(P, q, A, l, u, 'verbose', true) % Solve the problem using the Matlab solver -mat_results = solver.solve() - -% Generate the embedded code with the default options -solver.codegen('osqp_demo') - -% Solve the problem using the generated code -[x, y, status_val, iter, run_time] = emosqp('solve') +results = solver.solve(); \ No newline at end of file diff --git a/examples/osqpmatlab.m b/examples/osqpmatlab.m index dfa2ba0..5a26e16 100644 --- a/examples/osqpmatlab.m +++ b/examples/osqpmatlab.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + function [x, y, cost, status, iter] = osqpmatlab(problem, warm_start, settings) %#codegen % OSQPMATLAB Pure Matlab implementation of the OSQP solver % diff --git a/make_osqp.m b/make_osqp.m deleted file mode 100644 index ce3a271..0000000 --- a/make_osqp.m +++ /dev/null @@ -1,307 +0,0 @@ -function make_osqp(varargin) -% Matlab MEX makefile for OSQP. -% -% MAKE_OSQP(VARARGIN) is a make file for OSQP solver. It -% builds OSQP and its components from source. -% -% WHAT is the last element of VARARGIN and cell array of strings, -% with the following options: -% -% {}, '' (empty string) or 'all': build all components and link. -% -% 'osqp': builds the OSQP solver using CMake -% -% 'osqp_mex': builds the OSQP mex interface and links it to the OSQP -% library -% -% VARARGIN{1:NARGIN-1} specifies the optional flags passed to the compiler -% -% Additional commands: -% -% 'clean': delete all object files (.o and .obj) -% 'purge' : same as above, and also delete the mex files. - - -if( nargin == 0 ) - what = {'all'}; - verbose = false; -elseif ( nargin == 1 && ismember('-verbose', varargin) ) - what = {'all'}; - verbose = true; -else - what = varargin{nargin}; - if(isempty(strfind(what, 'all')) && ... - isempty(strfind(what, 'osqp')) && ... - isempty(strfind(what, 'osqp_mex')) && ... - isempty(strfind(what, 'clean')) && ... - isempty(strfind(what, 'purge'))) - fprintf('No rule to make target "%s", exiting.\n', what); - end - - verbose = ismember('-verbose', varargin); -end - -%% Try to unlock any pre-existing version of osqp_mex - -% this prevents compile errors if a user builds, runs osqp -% and then tries to recompile -if(mislocked('osqp_mex')) - munlock('osqp_mex'); -end - - - -%% Basic compile commands - -% Get make and mex commands -make_cmd = 'cmake --build .'; -mex_cmd = sprintf('mex -O -silent'); -mex_libs = ''; - - -% Add arguments to cmake and mex compiler -cmake_args = '-DMATLAB=ON'; -mexoptflags = '-DMATLAB'; - -% Add specific generators for windows linux or mac -if (ispc) - cmake_args = sprintf('%s %s', cmake_args, '-G "MinGW Makefiles"'); -else - cmake_args = sprintf('%s %s', cmake_args, '-G "Unix Makefiles"'); -end - -% Pass Matlab root to cmake -Matlab_ROOT = strrep(matlabroot, '\', '/'); -cmake_args = sprintf('%s %s%s%s', cmake_args, ... - '-DMatlab_ROOT_DIR="', Matlab_ROOT, '"'); - -% Add parameters options to mex and cmake -% CTRLC -if (ispc) - ut = fullfile(matlabroot, 'extern', 'lib', computer('arch'), ... - 'mingw64', 'libut.lib'); - mex_libs = sprintf('%s "%s"', mex_libs, ut); -else - mex_libs = sprintf('%s %s', mex_libs, '-lut'); -end -% Shared library loading -if (isunix && ~ismac) - mex_libs = sprintf('%s %s', mex_libs, '-ldl'); -end - -% Add large arrays support if computer is 64 bit and a pre-2018 version -% Release R2018a corresponds to Matlab version 9.4 -if (~isempty(strfind(computer, '64')) && verLessThan('matlab', '9.4')) - mexoptflags = sprintf('%s %s', mexoptflags, '-largeArrayDims'); -end - -%Force Matlab to respect old-style usage of mxGetPr in releases after 2018a, -%which use interleaved complex data. Note that the -R2017b flag is badly -%named since it indicates that non-interleaved complex data model is being used; -%it is not really specific to the release year -if ~verLessThan('matlab', '9.4') - mexoptflags = sprintf('%s %s', mexoptflags, '-R2017b'); -end - - -% Set optimizer flag -if (~ispc) - mexoptflags = sprintf('%s %s', mexoptflags, 'COPTIMFLAGS=''-O3'''); -end - -% Set library extension -lib_ext = '.a'; -lib_name = sprintf('libosqp%s', lib_ext); - - -% Set osqp directory and osqp_build directory -current_dir = pwd; -[makefile_path,~,~] = fileparts(which('make_osqp.m')); -osqp_dir = fullfile(makefile_path, 'osqp_sources'); -osqp_build_dir = fullfile(osqp_dir, 'build'); -qdldl_dir = fullfile(osqp_dir, 'lin_sys', 'direct', 'qdldl'); -cg_sources_dir = fullfile('.', 'codegen', 'sources'); - -% Include directory -inc_dir = [ - fullfile(sprintf(' -I%s', osqp_dir), 'include'), ... - sprintf(' -I%s', qdldl_dir), ... - fullfile(sprintf(' -I%s', qdldl_dir), 'qdldl_sources', 'include')]; - - -%% OSQP Solver -if( any(strcmpi(what,'osqp')) || any(strcmpi(what,'all')) ) - fprintf('Compiling OSQP solver...'); - - % Create build directory and go inside - if exist(osqp_build_dir, 'dir') - rmdir(osqp_build_dir, 's'); - end - mkdir(osqp_build_dir); - cd(osqp_build_dir); - - % Extend path for CMake mac (via Homebrew) - PATH = getenv('PATH'); - if ((ismac) && (isempty(strfind(PATH, '/usr/local/bin')))) - setenv('PATH', [PATH ':/usr/local/bin']); - end - - % Compile static library with CMake - [status, output] = system(sprintf('%s %s ..', 'cmake', cmake_args)); - if(status) - fprintf('\n'); - disp(output); - error('Error configuring CMake environment'); - elseif(verbose) - fprintf('\n'); - disp(output); - end - - [status, output] = system(sprintf('%s %s', make_cmd, '--target osqpstatic')); - if (status) - fprintf('\n'); - disp(output); - error('Error compiling OSQP'); - elseif(verbose) - fprintf('\n'); - disp(output); - end - - - % Change directory back to matlab interface - cd(makefile_path); - - % Copy static library to current folder - lib_origin = fullfile(osqp_build_dir, 'out', lib_name); - copyfile(lib_origin, lib_name); - - fprintf('\t\t\t\t\t\t[done]\n'); - -end - -%% osqpmex -if( any(strcmpi(what,'osqp_mex')) || any(strcmpi(what,'all')) ) - % Compile interface - fprintf('Compiling and linking osqpmex...'); - - % Compile command - %cmd = sprintf('%s %s %s %s osqp_mex.cpp', mex_cmd, mexoptflags, inc_dir, lib_name); - cmd = sprintf('%s %s %s %s osqp_mex.cpp %s', ... - mex_cmd, mexoptflags, inc_dir, lib_name, mex_libs); - - % Compile - eval(cmd); - fprintf('\t\t\t\t\t[done]\n'); - -end - - -%% codegen -if( any(strcmpi(what,'codegen')) || any(strcmpi(what,'all')) ) - fprintf('Copying source files for codegen...'); - - % Copy C files - cg_src_dir = fullfile(cg_sources_dir, 'src'); - if ~exist(cg_src_dir, 'dir') - mkdir(cg_src_dir); - end - cdirs = {fullfile(osqp_dir, 'src'),... - fullfile(qdldl_dir),... - fullfile(qdldl_dir, 'qdldl_sources', 'src')}; - for j = 1:length(cdirs) - cfiles = dir(fullfile(cdirs{j},'*.c')); - for i = 1 : length(cfiles) - if ~any(strcmp(cfiles(i).name, {'cs.c', 'ctrlc.c', 'lin_sys.c', 'polish.c'})) - copyfile(fullfile(cdirs{j}, cfiles(i).name), ... - fullfile(cg_src_dir, cfiles(i).name)); - end - end - end - - % Copy H files - cg_include_dir = fullfile(cg_sources_dir, 'include'); - if ~exist(cg_include_dir, 'dir') - mkdir(cg_include_dir); - end - hdirs = {fullfile(osqp_dir, 'include'),... - fullfile(qdldl_dir),... - fullfile(qdldl_dir, 'qdldl_sources', 'include')}; - for j = 1:length(hdirs) - hfiles = dir(fullfile(hdirs{j},'*.h')); - for i = 1 : length(hfiles) - if ~any(strcmp(hfiles(i).name, {'qdldl_types.h', 'osqp_configure.h', ... - 'cs.h', 'ctrlc.h', 'lin_sys.h', 'polish.h'})) - copyfile(fullfile(hdirs{j}, hfiles(i).name), ... - fullfile(cg_include_dir, hfiles(i).name)); - end - end - end - - % Copy configure files - cg_configure_dir = fullfile(cg_sources_dir, 'configure'); - if ~exist(cg_configure_dir, 'dir') - mkdir(cg_configure_dir); - end - configure_dirs = {fullfile(osqp_dir, 'configure'),... - fullfile(qdldl_dir, 'qdldl_sources', 'configure')}; - for j = 1:length(configure_dirs) - configure_files = dir(fullfile(configure_dirs{j},'*.h.in')); - for i = 1 : length(configure_files) - copyfile(fullfile(configure_dirs{j}, configure_files(i).name), ... - fullfile(cg_configure_dir, configure_files(i).name)); - end - end - - % Copy cmake files - copyfile(fullfile(osqp_dir, 'src', 'CMakeLists.txt'), ... - fullfile(cg_src_dir, 'CMakeLists.txt')); - copyfile(fullfile(osqp_dir, 'include', 'CMakeLists.txt'), ... - fullfile(cg_include_dir, 'CMakeLists.txt')); - - fprintf('\t\t\t\t\t[done]\n'); - -end - - -%% clean -if( any(strcmpi(what,'clean')) || any(strcmpi(what,'purge')) ) - fprintf('Cleaning mex files and library...'); - - % Delete mex file - mexfiles = dir(['*.', mexext]); - for i = 1 : length(mexfiles) - delete(mexfiles(i).name); - end - - % Delete static library - lib_full_path = fullfile(makefile_path, lib_name); - if( exist(lib_full_path,'file') ) - delete(lib_full_path); - end - - fprintf('\t\t\t[done]\n'); -end - - -%% purge -if( any(strcmpi(what,'purge')) ) - fprintf('Cleaning OSQP build and codegen directories...'); - - % Delete OSQP build directory - if exist(osqp_build_dir, 'dir') - rmdir(osqp_build_dir, 's'); - end - - % Delete codegen files - if exist(cg_sources_dir, 'dir') - rmdir(cg_sources_dir, 's'); - end - - fprintf('\t\t[done]\n'); -end - - -%% Go back to the original directory -cd(current_dir); - -end diff --git a/osqp.m b/osqp.m deleted file mode 100755 index acc653a..0000000 --- a/osqp.m +++ /dev/null @@ -1,677 +0,0 @@ -classdef osqp < handle - % osqp interface class for OSQP solver - % This class provides a complete interface to the C implementation - % of the OSQP solver. - % - % osqp Properties: - % objectHandle - pointer to the C structure of OSQP solver - % - % osqp Methods: - % - % setup - configure solver with problem data - % solve - solve the QP - % update - modify problem vectors - % warm_start - set warm starting variables x and y - % - % default_settings - create default settings structure - % current_settings - get the current solver settings structure - % update_settings - update the current solver settings structure - % - % get_dimensions - get the number of variables and constraints - % version - return OSQP version - % constant - return a OSQP internal constant - % - % codegen - generate embeddable C code for the problem - - - properties (SetAccess = private, Hidden = true) - objectHandle % Handle to underlying C instance - end - methods(Static) - %% - function out = default_settings() - % DEFAULT_SETTINGS get the default solver settings structure - out = osqp_mex('default_settings', 'static'); - - % Convert linsys solver to string - out.linsys_solver = linsys_solver_to_string(out.linsys_solver); - - end - - %% - function out = constant(constant_name) - % CONSTANT Return solver constant - % C = CONSTANT(CONSTANT_NAME) return constant called CONSTANT_NAME - out = osqp_mex('constant', 'static', constant_name); - end - - %% - function out = version() - % Return OSQP version - out = osqp_mex('version', 'static'); - end - - end - methods - %% Constructor - Create a new solver instance - function this = osqp(varargin) - % Construct OSQP solver class - this.objectHandle = osqp_mex('new', varargin{:}); - end - - %% Destructor - destroy the solver instance - function delete(this) - % Destroy OSQP solver class - osqp_mex('delete', this.objectHandle); - end - - %% - function out = current_settings(this) - % CURRENT_SETTINGS get the current solver settings structure - out = osqp_mex('current_settings', this.objectHandle); - - % Convert linsys solver to string - out.linsys_solver = linsys_solver_to_string(out.linsys_solver); - - end - - %% - function update_settings(this,varargin) - % UPDATE_SETTINGS update the current solver settings structure - - %second input 'false' means that this is *not* a settings - %initialization, so some parameter/values will be disallowed - newSettings = validateSettings(this,false,varargin{:}); - - %write the solver settings. C-mex does not check input - %data or protect against disallowed parameter modifications - osqp_mex('update_settings', this.objectHandle, newSettings); - - end - - %% - function [n,m] = get_dimensions(this) - % GET_DIMENSIONS get the number of variables and constraints - - [n,m] = osqp_mex('get_dimensions', this.objectHandle); - - end - - %% - function update(this,varargin) - % UPDATE modify the linear cost term and/or lower and upper bounds - - %second input 'false' means that this is *not* a settings - %initialization, so some parameter/values will be disallowed - allowedFields = {'q','l','u','Px','Px_idx','Ax','Ax_idx'}; - - if(isempty(varargin)) - return; - elseif(length(varargin) == 1) - if(~isstruct(varargin{1})) - error('Single input should be a structure with new problem data'); - else - newData = varargin{1}; - end - else % param / value style assumed - newData = struct(varargin{:}); - end - - %check for unknown fields - newFields = fieldnames(newData); - badFieldsIdx = find(~ismember(newFields,allowedFields)); - if(~isempty(badFieldsIdx)) - error('Unrecognized input field ''%s'' detected',newFields{badFieldsIdx(1)}); - end - - %get all of the terms. Nonexistent fields will be passed - %as empty mxArrays - try q = double(full(newData.q(:))); catch q = []; end - try l = double(full(newData.l(:))); catch l = []; end - try u = double(full(newData.u(:))); catch u = []; end - try Px = double(full(newData.Px(:))); catch Px = []; end - try Px_idx = double(full(newData.Px_idx(:))); catch Px_idx = []; end - try Ax = double(full(newData.Ax(:))); catch Ax = []; end - try Ax_idx = double(full(newData.Ax_idx(:))); catch Ax_idx = []; end - - [n,m] = get_dimensions(this); - - assert(isempty(q) || length(q) == n, 'input ''q'' is the wrong size'); - assert(isempty(l) || length(l) == m, 'input ''u'' is the wrong size'); - assert(isempty(u) || length(u) == m, 'input ''l'' is the wrong size'); - assert(isempty(Px) || isempty(Px_idx) || length(Px) == length(Px_idx), ... - 'inputs ''Px'' and ''Px_idx'' must be the same size'); - assert(isempty(Ax) || isempty(Ax_idx) || length(Ax) == length(Ax_idx), ... - 'inputs ''Ax'' and ''Ax_idx'' must be the same size'); - - % Adjust index of Px_idx and Ax_idx to match 0-based indexing - % in C - if (~isempty(Px_idx)) - Px_idx = Px_idx - 1; - end - if (~isempty(Ax_idx)) - Ax_idx = Ax_idx - 1; - end - - % Convert infinity values to OSQP_INFTY - if (~isempty(u)) - u = min(u, osqp.constant('OSQP_INFTY')); - end - if (~isempty(l)) - l = max(l, -osqp.constant('OSQP_INFTY')); - end - - %write the new problem data. C-mex does not protect - %against unknown fields, but will handle empty values - osqp_mex('update', this.objectHandle, ... - q, l, u, Px, Px_idx, length(Px), Ax, Ax_idx, length(Ax)); - - end - - %% - function varargout = setup(this, varargin) - % SETUP configure solver with problem data - % - % setup(P,q,A,l,u,options) - - nargin = length(varargin); - - %dimension checks on user data. Mex function does not - %perform any checks on inputs, so check everything here - assert(nargin >= 5, 'incorrect number of inputs'); - [P,q,A,l,u] = deal(varargin{1:5}); - - % - % Get problem dimensions - % - - % Get number of variables n - if (isempty(P)) - if (~isempty(q)) - n = length(q); - else - if (~isempty(A)) - n = size(A, 2); - else - error('The problem does not have any variables'); - end - end - else - n = size(P, 1); - end - - % Get number of constraints m - if (isempty(A)) - m = 0; - else - m = size(A, 1); - assert(size(A, 2) == n, 'Incorrect dimension of A'); - end - - % - % Create sparse matrices and full vectors if they are empty - % - - if (isempty(P)) - P = sparse(n, n); - else - P = sparse(P); - end - if (~istriu(P)) - P = triu(P); - end - if (isempty(q)) - q = zeros(n, 1); - else - q = full(q(:)); - end - - % Create proper constraints if they are not passed - if (isempty(A) && (~isempty(l) || ~isempty(u))) || ... - (~isempty(A) && (isempty(l) && isempty(u))) - error('A must be supplied together with at least one bound l or u'); - end - - if (~isempty(A) && isempty(l)) - l = -Inf(m, 1); - end - - if (~isempty(A) && isempty(u)) - u = Inf(m, 1); - end - - if (isempty(A)) - A = sparse(m, n); - l = -Inf(m, 1); - u = Inf(m, 1); - else - l = full(l(:)); - u = full(u(:)); - A = sparse(A); - end - - - % - % Check vector dimensions (not checked from the C solver) - % - - assert(length(q) == n, 'Incorrect dimension of q'); - assert(length(l) == m, 'Incorrect dimension of l'); - assert(length(u) == m, 'Incorrect dimension of u'); - - % - % Convert infinity values to OSQP_INFINITY - % - u = min(u, osqp.constant('OSQP_INFTY')); - l = max(l, -osqp.constant('OSQP_INFTY')); - - - %make a settings structure from the remainder of the arguments. - %'true' means that this is a settings initialization, so all - %parameter/values are allowed. No extra inputs will result - %in default settings being passed back - theSettings = validateSettings(this,true,varargin{6:end}); - - [varargout{1:nargout}] = osqp_mex('setup', this.objectHandle, n,m,P,q,A,l,u,theSettings); - - end - - - %% - - function warm_start(this, varargin) - % WARM_START warm start primal and/or dual variables - % - % warm_start('x', x, 'y', y) - % - % or warm_start('x', x) - % or warm_start('y', y) - - - % Get problem dimensions - [n, m] = get_dimensions(this); - - % Get data - allowedFields = {'x','y'}; - - if(isempty(varargin)) - return; - elseif(length(varargin) == 1) - if(~isstruct(varargin{1})) - error('Single input should be a structure with new problem data'); - else - newData = varargin{1}; - end - else % param / value style assumed - newData = struct(varargin{:}); - end - - %check for unknown fields - newFields = fieldnames(newData); - badFieldsIdx = find(~ismember(newFields,allowedFields)); - if(~isempty(badFieldsIdx)) - error('Unrecognized input field ''%s'' detected',newFields{badFieldsIdx(1)}); - end - - %get all of the terms. Nonexistent fields will be passed - %as empty mxArrays - try x = double(full(newData.x(:))); catch x = []; end - try y = double(full(newData.y(:))); catch y = []; end - - % Check dimensions - assert(isempty(x) || length(x) == n, 'input ''x'' is the wrong size'); - assert(isempty(y) || length(y) == m, 'input ''y'' is the wrong size'); - - - % Decide which function to call - if (~isempty(x) && isempty(y)) - osqp_mex('warm_start_x', this.objectHandle, x); - return; - end - - if (isempty(x) && ~isempty(y)) - osqp_mex('warm_start_y', this.objectHandle, y); - end - - if (~isempty(x) && ~isempty(y)) - osqp_mex('warm_start', this.objectHandle, x, y); - end - - if (isempty(x) && isempty(y)) - error('Unrecognized fields'); - end - - end - - %% - function varargout = solve(this, varargin) - % SOLVE solve the QP - - nargoutchk(0,1); %either return nothing (but still solve), or a single output structure - [out.x, out.y, out.prim_inf_cert, out.dual_inf_cert, out.info] = osqp_mex('solve', this.objectHandle); - if(nargout) - varargout{1} = out; - end - return; - end - - %% - function codegen(this, target_dir, varargin) - % CODEGEN generate C code for the parametric problem - % - % codegen(target_dir,options) - - % Parse input arguments - p = inputParser; - defaultProject = ''; - expectedProject = {'', 'Makefile', 'MinGW Makefiles', 'Unix Makefiles', 'CodeBlocks', 'Xcode'}; - defaultParams = 'vectors'; - expectedParams = {'vectors', 'matrices'}; - defaultMexname = 'emosqp'; - defaultFloat = false; - defaultLong = true; - defaultFW = false; - - addRequired(p, 'target_dir', @isstr); - addParameter(p, 'project_type', defaultProject, ... - @(x) ischar(validatestring(x, expectedProject))); - addParameter(p, 'parameters', defaultParams, ... - @(x) ischar(validatestring(x, expectedParams))); - addParameter(p, 'mexname', defaultMexname, @isstr); - addParameter(p, 'FLOAT', defaultFloat, @islogical); - addParameter(p, 'LONG', defaultLong, @islogical); - addParameter(p, 'force_rewrite', defaultFW, @islogical); - - parse(p, target_dir, varargin{:}); - - % Set internal variables - if strcmp(p.Results.parameters, 'vectors') - embedded = 1; - else - embedded = 2; - end - if p.Results.FLOAT - float_flag = 'ON'; - else - float_flag = 'OFF'; - end - if p.Results.LONG - long_flag = 'ON'; - else - long_flag = 'OFF'; - end - if strcmp(p.Results.project_type, 'Makefile') - if (ispc) - project_type = 'MinGW Makefiles'; % Windows - elseif (ismac || isunix) - project_type = 'Unix Makefiles'; % Unix - end - else - project_type = p.Results.project_type; - end - - % Check whether the specified directory already exists - if exist(target_dir, 'dir') - if p.Results.force_rewrite - rmdir(target_dir, 's'); - else - while(1) - prompt = sprintf('Directory "%s" already exists. Do you want to replace it? y/n [y]: ', target_dir); - str = input(prompt, 's'); - - if any(strcmpi(str, {'','y'})) - rmdir(target_dir, 's'); - break; - elseif strcmpi(str, 'n') - return; - end - end - end - end - - % Import OSQP path - [osqp_path,~,~] = fileparts(which('osqp.m')); - - % Add codegen directory to path - addpath(fullfile(osqp_path, 'codegen')); - - % Path of osqp module - cg_dir = fullfile(osqp_path, 'codegen'); - files_to_generate_path = fullfile(cg_dir, 'files_to_generate'); - - % Get workspace structure - work = osqp_mex('get_workspace', this.objectHandle); - - % Make target directory - fprintf('Creating target directories...\t\t\t\t\t'); - target_configure_dir = fullfile(target_dir, 'configure'); - target_include_dir = fullfile(target_dir, 'include'); - target_src_dir = fullfile(target_dir, 'src'); - - if ~exist(target_dir, 'dir') - mkdir(target_dir); - end - if ~exist(target_configure_dir, 'dir') - mkdir(target_configure_dir); - end - if ~exist(target_include_dir, 'dir') - mkdir(target_include_dir); - end - if ~exist(target_src_dir, 'dir') - mkdir(fullfile(target_src_dir, 'osqp')); - end - fprintf('[done]\n'); - - % Copy source files to target directory - fprintf('Copying OSQP source files...\t\t\t\t\t'); - cdir = fullfile(cg_dir, 'sources', 'src'); - cfiles = dir(fullfile(cdir, '*.c')); - for i = 1 : length(cfiles) - if embedded == 1 - % Do not copy kkt.c if embedded is 1 - if ~strcmp(cfiles(i).name, 'kkt.c') - copyfile(fullfile(cdir, cfiles(i).name), ... - fullfile(target_src_dir, 'osqp', cfiles(i).name)); - end - else - copyfile(fullfile(cdir, cfiles(i).name), ... - fullfile(target_src_dir, 'osqp', cfiles(i).name)); - end - end - configure_dir = fullfile(cg_dir, 'sources', 'configure'); - configure_files = dir(fullfile(configure_dir, '*.h.in')); - for i = 1 : length(configure_files) - copyfile(fullfile(configure_dir, configure_files(i).name), ... - fullfile(target_configure_dir, configure_files(i).name)); - end - hdir = fullfile(cg_dir, 'sources', 'include'); - hfiles = dir(fullfile(hdir, '*.h')); - for i = 1 : length(hfiles) - if embedded == 1 - % Do not copy kkt.h if embedded is 1 - if ~strcmp(hfiles(i).name, 'kkt.h') - copyfile(fullfile(hdir, hfiles(i).name), ... - fullfile(target_include_dir, hfiles(i).name)); - end - else - copyfile(fullfile(hdir, hfiles(i).name), ... - fullfile(target_include_dir, hfiles(i).name)); - end - end - - % Copy cmake files - copyfile(fullfile(cdir, 'CMakeLists.txt'), ... - fullfile(target_src_dir, 'osqp', 'CMakeLists.txt')); - copyfile(fullfile(hdir, 'CMakeLists.txt'), ... - fullfile(target_include_dir, 'CMakeLists.txt')); - fprintf('[done]\n'); - - % Copy example.c - copyfile(fullfile(files_to_generate_path, 'example.c'), target_src_dir); - - % Render CMakeLists.txt - fidi = fopen(fullfile(files_to_generate_path, 'CMakeLists.txt'),'r'); - fido = fopen(fullfile(target_dir, 'CMakeLists.txt'),'w'); - while ~feof(fidi) - l = fgetl(fidi); % read line - % Replace EMBEDDED_FLAG in CMakeLists.txt by a numerical value - newl = strrep(l, 'EMBEDDED_FLAG', num2str(embedded)); - fprintf(fido, '%s\n', newl); - end - fclose(fidi); - fclose(fido); - - % Render workspace.h and workspace.c - work_hfile = fullfile(target_include_dir, 'workspace.h'); - work_cfile = fullfile(target_src_dir, 'osqp', 'workspace.c'); - fprintf('Generating workspace.h/.c...\t\t\t\t\t\t'); - render_workspace(work, work_hfile, work_cfile, embedded); - fprintf('[done]\n'); - - % Create project - if ~isempty(project_type) - - % Extend path for CMake mac (via Homebrew) - PATH = getenv('PATH'); - if ((ismac) && (isempty(strfind(PATH, '/usr/local/bin')))) - setenv('PATH', [PATH ':/usr/local/bin']); - end - - fprintf('Creating project...\t\t\t\t\t\t\t\t'); - orig_dir = pwd; - cd(target_dir); - mkdir('build') - cd('build'); - cmd = sprintf('cmake -G "%s" ..', project_type); - [status, output] = system(cmd); - if(status) - fprintf('\n'); - fprintf(output); - error('Error configuring CMake environment'); - else - fprintf('[done]\n'); - end - cd(orig_dir); - end - - % Make mex interface to the generated code - mex_cfile = fullfile(files_to_generate_path, 'emosqp_mex.c'); - make_emosqp(target_dir, mex_cfile, embedded, float_flag, long_flag); - - % Rename the mex file - old_mexfile = ['emosqp_mex.', mexext]; - new_mexfile = [p.Results.mexname, '.', mexext]; - movefile(old_mexfile, new_mexfile); - - end - - end -end - - - -function currentSettings = validateSettings(this,isInitialization,varargin) - -%don't allow these fields to be changed -unmodifiableFields = {'scaling', 'linsys_solver'}; - -%get the current settings -if(isInitialization) - currentSettings = osqp_mex('default_settings', this.objectHandle); -else - currentSettings = osqp_mex('current_settings', this.objectHandle); -end - -%no settings passed -> return defaults -if(isempty(varargin)) - return; -end - -%check for structure style input -if(isstruct(varargin{1})) - newSettings = varargin{1}; - assert(length(varargin) == 1, 'too many input arguments'); -else - newSettings = struct(varargin{:}); -end - -%get the osqp settings fields -currentFields = fieldnames(currentSettings); - -%get the requested fields in the update -newFields = fieldnames(newSettings); - -%check for unknown parameters -badFieldsIdx = find(~ismember(newFields,currentFields)); -if(~isempty(badFieldsIdx)) - error('Unrecognized solver setting ''%s'' detected',newFields{badFieldsIdx(1)}); -end - -%convert linsys_solver string to integer -if ismember('linsys_solver',newFields) - if ~ischar(newSettings.linsys_solver) - error('Setting linsys_solver is required to be a string.'); - end - % Convert linsys_solver to number - newSettings.linsys_solver = string_to_linsys_solver(newSettings.linsys_solver); -end - - -%check for disallowed fields if this in not an initialization call -if(~isInitialization) - badFieldsIdx = find(ismember(newFields,unmodifiableFields)); - for i = badFieldsIdx(:)' - if(~isequal(newSettings.(newFields{i}),currentSettings.(newFields{i}))) - error('Solver setting ''%s'' can only be changed at solver initialization.', newFields{i}); - end - end -end - - -%check that everything is a nonnegative scalar (this check is already -%performed in C) -% for i = 1:length(newFields) -% val = double(newSettings.(newFields{i})); -% assert(isscalar(val) & isnumeric(val) & val >= 0, ... -% 'Solver setting ''%s'' not specified as nonnegative scalar', newFields{i}); -% end - -%everything checks out - merge the newSettings into the current ones -for i = 1:length(newFields) - currentSettings.(newFields{i}) = double(newSettings.(newFields{i})); -end - - -end - -function [linsys_solver_string] = linsys_solver_to_string(linsys_solver) -% Convert linear systme solver integer to stringh -switch linsys_solver - case osqp.constant('QDLDL_SOLVER') - linsys_solver_string = 'qdldl'; - case osqp.constant('MKL_PARDISO_SOLVER') - linsys_solver_string = 'mkl pardiso'; - otherwise - error('Unrecognized linear system solver.'); -end -end - - - -function [linsys_solver] = string_to_linsys_solver(linsys_solver_string) -linsys_solver_string = lower(linsys_solver_string); -switch linsys_solver_string - case 'qdldl' - linsys_solver = osqp.constant('QDLDL_SOLVER'); - case 'mkl pardiso' - linsys_solver = osqp.constant('MKL_PARDISO_SOLVER'); - % Default solver: QDLDL - case '' - linsys_solver = osqp.constant('QDLDL_SOLVER'); - otherwise - warning('Linear system solver not recognized. Using default solver QDLDL.') - linsys_solver = osqp.constant('QDLDL_SOLVER'); -end -end - - diff --git a/osqp_mex.cpp b/osqp_mex.cpp deleted file mode 100755 index 5506e1d..0000000 --- a/osqp_mex.cpp +++ /dev/null @@ -1,1198 +0,0 @@ -#include "mex.h" -#include "matrix.h" -#include "osqp_mex.hpp" -#include "osqp.h" -#include "ctrlc.h" // Needed for interrupt -#include "qdldl_interface.h" // To extract workspace for codegen - -// all of the OSQP_INFO fieldnames as strings -const char* OSQP_INFO_FIELDS[] = {"iter", //c_int - "status" , //char* - "status_val" , //c_int - "status_polish", //c_int - "obj_val", //c_float - "pri_res", //c_float - "dua_res", //c_float - "setup_time", //c_float, only used if PROFILING - "solve_time", //c_float, only used if PROFILING - "update_time", //c_float, only used if PROFILING - "polish_time", //c_float, only used if PROFILING - "run_time", //c_float, only used if PROFILING - "rho_updates", //c_int - "rho_estimate"}; //c_float - -const char* OSQP_SETTINGS_FIELDS[] = {"rho", //c_float - "sigma", //c_float - "scaling", //c_int - "adaptive_rho", //c_int - "adaptive_rho_interval", //c_int - "adaptive_rho_tolerance", //c_float - "adaptive_rho_fraction", //c_float - "max_iter", //c_int - "eps_abs", //c_float - "eps_rel", //c_float - "eps_prim_inf", //c_float - "eps_dual_inf", //c_float - "alpha", //c_float - "linsys_solver", //c_int - "delta", //c_float - "polish", //c_int - "polish_refine_iter", //c_int - "verbose", //c_int - "scaled_termination", //c_int - "check_termination", //c_int - "warm_start", //c_int - "time_limit"}; //c_float - -const char* CSC_FIELDS[] = {"nzmax", //c_int - "m", //c_int - "n", //c_int - "p", //c_int* - "i", //c_int* - "x", //c_float* - "nz"}; //c_int - -const char* OSQP_DATA_FIELDS[] = {"n", //c_int - "m", //c_int - "P", //csc - "A", //csc - "q", //c_float* - "l", //c_float* - "u"}; //c_float* - -const char* LINSYS_SOLVER_FIELDS[] = {"L", //csc - "Dinv", //c_float* - "P", //c_int* - "bp", //c_float* - "sol", //c_float* - "rho_inv_vec", //c_float* - "sigma", //c_float - "polish", //c_int - "n", //c_int - "m", //c_int - "Pdiag_idx", //c_int* - "Pdiag_n", //c_int - "KKT", //csc - "PtoKKT", //c_int* - "AtoKKT", //c_int* - "rhotoKKT", //c_int* - "D", //c_float* - "etree", //c_int* - "Lnz", //c_int* - "iwork", //c_int* - "bwork", //c_int* - "fwork"}; //c_float* - -const char* OSQP_SCALING_FIELDS[] = {"c", //c_float - "D", //c_float* - "E", //c_float* - "cinv", //c_float - "Dinv", //c_float* - "Einv"}; //c_float* - -const char* OSQP_RHO_VECTORS_FIELDS[] = {"rho_vec", //c_float* - "rho_inv_vec", //c_float* - "constr_type"}; //c_int* - -const char* OSQP_WORKSPACE_FIELDS[] = {"rho_vectors", - "data", - "linsys_solver", - "scaling", - "settings"}; - - -#define NEW_SETTINGS_TOL (1e-10) - -// wrapper class for all osqp data and settings -class OsqpData -{ -public: - OsqpData(){ - work = NULL; - } - OSQPWorkspace * work; // Workspace -}; - -// internal utility functions -void castToDoubleArr(c_float *arr, double* arr_out, c_int len); -void setToNaN(double* arr_out, c_int len); -void copyMxStructToSettings(const mxArray*, OSQPSettings*); -void copyUpdatedSettingsToWork(const mxArray*, OsqpData*); -void castCintToDoubleArr(c_int *arr, double* arr_out, c_int len); -c_int* copyToCintVector(mwIndex * vecData, c_int numel); -c_int* copyDoubleToCintVector(double* vecData, c_int numel); -c_float* copyToCfloatVector(double * vecData, c_int numel); -mxArray* copyInfoToMxStruct(OSQPInfo* info); -mxArray* copySettingsToMxStruct(OSQPSettings* settings); -mxArray* copyCscMatrixToMxStruct(csc* M); -mxArray* copyDataToMxStruct(OSQPWorkspace* work); -mxArray* copyLinsysSolverToMxStruct(OSQPWorkspace* work); -mxArray* copyScalingToMxStruct(OSQPWorkspace * work); -mxArray* copyWorkToMxStruct(OSQPWorkspace* work); - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ OsqpData* osqpData; // OSQP data identifier - - // Exitflag - c_int exitflag = 0; - - // Static string for static methods - char stat_string[64]; - - // Get the command string - char cmd[64]; - if (nrhs < 1 || mxGetString(prhs[0], cmd, sizeof(cmd))) - mexErrMsgTxt("First input should be a command string less than 64 characters long."); - - // new object - if (!strcmp("new", cmd)) { - // Check parameters - if (nlhs != 1){ - mexErrMsgTxt("New: One output expected."); - } - // Return a handle to a new C++ wrapper instance - plhs[0] = convertPtr2Mat(new OsqpData); - return; - } - - // Check for a second input, which should be the class instance handle or string 'static' - if (nrhs < 2) - mexErrMsgTxt("Second input should be a class instance handle or the string 'static'."); - - if(mxGetString(prhs[1], stat_string, sizeof(stat_string))){ - // If we are dealing with non-static methods, get the class instance pointer from the second input - osqpData = convertMat2Ptr(prhs[1]); - } else { - if (strcmp("static", stat_string)){ - mexErrMsgTxt("Second argument for static functions is string 'static'"); - } - } - - // delete the object and its data - if (!strcmp("delete", cmd)) { - - //clean up the problem workspace - if(osqpData->work){ - osqp_cleanup(osqpData->work); - } - - //clean up the handle object - destroyObject(prhs[1]); - // Warn if other commands were ignored - if (nlhs != 0 || nrhs != 2) - mexWarnMsgTxt("Delete: Unexpected arguments ignored."); - return; - } - - // report the current settings - if (!strcmp("current_settings", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver is uninitialized. No settings have been configured."); - } - //report the current settings - plhs[0] = copySettingsToMxStruct(osqpData->work->settings); - return; - } - - - // return workspace structure - if (!strcmp("get_workspace", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver is uninitialized. No data have been configured."); - } - - //throw an error if linear systems solver is different than qdldl - if(osqpData->work->linsys_solver->type != QDLDL_SOLVER){ - mexErrMsgTxt("Solver setup was not performed using QDLDL! Please perform setup with linsys_solver as QDLDL."); - } - - //return data - plhs[0] = copyWorkToMxStruct(osqpData->work); - return; - } - - // write_settings - if (!strcmp("update_settings", cmd)) { - - //overwrite the current settings. Mex function is responsible - //for disallowing overwrite of selected settings after initialization, - //and for all error checking - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver is uninitialized. No settings have been configured."); - } - - copyUpdatedSettingsToWork(prhs[2],osqpData); - return; - } - - // report the default settings - if (!strcmp("default_settings", cmd)) { - - // Warn if other commands were ignored - if (nrhs > 2) - mexWarnMsgTxt("Default settings: unexpected number of arguments."); - - - //Create a Settings structure in default form and report the results - //Useful for external solver packages (e.g. Yalmip) that want to - //know which solver settings are supported - OSQPSettings* defaults = (OSQPSettings *)mxCalloc(1,sizeof(OSQPSettings)); - osqp_set_default_settings(defaults); - plhs[0] = copySettingsToMxStruct(defaults); - mxFree(defaults); - return; - } - - // setup - if (!strcmp("setup", cmd)) { - - //throw an error if this is called more than once - if(osqpData->work){ - mexErrMsgTxt("Solver is already initialized with problem data."); - } - - //Create data and settings containers - OSQPSettings* settings = (OSQPSettings *)mxCalloc(1,sizeof(OSQPSettings)); - OSQPData* data = (OSQPData *)mxCalloc(1,sizeof(OSQPData)); - - // handle the problem data first. Matlab-side - // class wrapper is responsible for ensuring that - // P and A are sparse matrices, everything - // else is a dense vector and all inputs are - // of compatible dimension - - // Get mxArrays - const mxArray* P = prhs[4]; - const mxArray* q = prhs[5]; - const mxArray* A = prhs[6]; - const mxArray* l = prhs[7]; - const mxArray* u = prhs[8]; - - // Create Data Structure - data->n = (c_int) mxGetScalar(prhs[2]); - data->m = (c_int) mxGetScalar(prhs[3]); - data->q = copyToCfloatVector(mxGetPr(q), data->n); - data->l = copyToCfloatVector(mxGetPr(l), data->m); - data->u = copyToCfloatVector(mxGetPr(u), data->m); - - // Matrix P: nnz = P->p[n] - c_int * Pp = copyToCintVector(mxGetJc(P), data->n + 1); - c_int * Pi = copyToCintVector(mxGetIr(P), Pp[data->n]); - c_float * Px = copyToCfloatVector(mxGetPr(P), Pp[data->n]); - data->P = csc_matrix(data->n, data->n, Pp[data->n], Px, Pi, Pp); - - // Matrix A: nnz = A->p[n] - c_int * Ap = copyToCintVector(mxGetJc(A), data->n + 1); - c_int * Ai = copyToCintVector(mxGetIr(A), Ap[data->n]); - c_float * Ax = copyToCfloatVector(mxGetPr(A), Ap[data->n]); - data->A = csc_matrix(data->m, data->n, Ap[data->n], Ax, Ai, Ap); - - // Create Settings - const mxArray* mxSettings = prhs[9]; - if(mxIsEmpty(mxSettings)){ - // use defaults - osqp_set_default_settings(settings); - } else { - //populate settings structure from mxArray input - copyMxStructToSettings(mxSettings, settings); - } - - // Setup workspace - exitflag = osqp_setup(&(osqpData->work), data, settings); - - //cleanup temporary structures - // Data - if (data->q) c_free(data->q); - if (data->l) c_free(data->l); - if (data->u) c_free(data->u); - if (Px) c_free(Px); - if (Pi) c_free(Pi); - if (Pp) c_free(Pp); - if (data->P) c_free(data->P); - if (Ax) c_free(Ax); - if (Ai) c_free(Ai); - if (Ap) c_free(Ap); - if (data->A) c_free(data->A); - if (data) mxFree(data); - // Settings - if (settings) mxFree(settings); - - // Report error (if any) - if(exitflag){ - mexErrMsgTxt("Invalid problem setup"); - } - - return; - - } - - // get #constraints and variables - if (!strcmp("get_dimensions", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver has not been initialized."); - } - - plhs[0] = mxCreateDoubleScalar(osqpData->work->data->n); - plhs[1] = mxCreateDoubleScalar(osqpData->work->data->m); - - return; - } - - // report the version - if (!strcmp("version", cmd)) { - - plhs[0] = mxCreateString(osqp_version()); - - return; - } - - // update linear cost and bounds - if (!strcmp("update", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver has not been initialized."); - } - - // Fill q, l, u - const mxArray *q = prhs[2]; - const mxArray *l = prhs[3]; - const mxArray *u = prhs[4]; - const mxArray *Px = prhs[5]; - const mxArray *Px_idx = prhs[6]; - const mxArray *Ax = prhs[8]; - const mxArray *Ax_idx = prhs[9]; - - int Px_n, Ax_n; - Px_n = mxGetScalar(prhs[7]); - Ax_n = mxGetScalar(prhs[10]); - - // Copy vectors to ensure they are cast as c_float - c_float *q_vec; - c_float *l_vec; - c_float *u_vec; - c_float *Px_vec; - c_float *Ax_vec; - c_int *Px_idx_vec = NULL; - c_int *Ax_idx_vec = NULL; - if(!mxIsEmpty(q)){ - q_vec = copyToCfloatVector(mxGetPr(q), - osqpData->work->data->n); - } - if(!mxIsEmpty(l)){ - l_vec = copyToCfloatVector(mxGetPr(l), - osqpData->work->data->m); - } - if(!mxIsEmpty(u)){ - u_vec = copyToCfloatVector(mxGetPr(u), - osqpData->work->data->m); - } - if(!mxIsEmpty(Px)){ - Px_vec = copyToCfloatVector(mxGetPr(Px), Px_n); - } - if(!mxIsEmpty(Ax)){ - Ax_vec = copyToCfloatVector(mxGetPr(Ax), Ax_n); - } - if(!mxIsEmpty(Px_idx)){ - Px_idx_vec = copyDoubleToCintVector(mxGetPr(Px_idx), Px_n); - } - if(!mxIsEmpty(Ax_idx)){ - Ax_idx_vec = copyDoubleToCintVector(mxGetPr(Ax_idx), Ax_n); - } - - if(!exitflag && !mxIsEmpty(q)){ - exitflag = osqp_update_lin_cost(osqpData->work, q_vec); - if(exitflag){ - exitflag = 1; - } - } - if(!exitflag && !mxIsEmpty(l) && !mxIsEmpty(u)){ - exitflag = osqp_update_bounds(osqpData->work, l_vec, u_vec); - if(exitflag){ - exitflag = 2; - } - } - else if(!exitflag && !mxIsEmpty(l)){ - exitflag = osqp_update_lower_bound(osqpData->work, l_vec); - if(exitflag){ - exitflag = 3; - } - } - else if(!exitflag && !mxIsEmpty(u)){ - exitflag = osqp_update_upper_bound(osqpData->work, u_vec); - if(exitflag){ - exitflag = 4; - } - } - if(!exitflag && !mxIsEmpty(Px) && !mxIsEmpty(Ax)){ - exitflag = osqp_update_P_A(osqpData->work, Px_vec, Px_idx_vec, Px_n, - Ax_vec, Ax_idx_vec, Ax_n); - if(exitflag){ - exitflag = 5; - } - } - else if(!exitflag && !mxIsEmpty(Px)){ - exitflag = osqp_update_P(osqpData->work, Px_vec, Px_idx_vec, Px_n); - if(exitflag){ - exitflag = 6; - } - } - else if(!exitflag && !mxIsEmpty(Ax)){ - exitflag = osqp_update_A(osqpData->work, Ax_vec, Ax_idx_vec, Ax_n); - if(exitflag){ - exitflag = 7; - } - } - - // Free vectors - if(!mxIsEmpty(q)) c_free(q_vec); - if(!mxIsEmpty(l)) c_free(l_vec); - if(!mxIsEmpty(u)) c_free(u_vec); - if(!mxIsEmpty(Px)) c_free(Px_vec); - if(!mxIsEmpty(Ax)) c_free(Ax_vec); - if(!mxIsEmpty(Px_idx)) c_free(Px_idx_vec); - if(!mxIsEmpty(Ax_idx)) c_free(Ax_idx_vec); - - // Report errors (if any) - switch (exitflag) { - case 1: - mexErrMsgTxt("Linear cost update error!"); - case 2: - mexErrMsgTxt("Bounds update error!"); - case 3: - mexErrMsgTxt("Lower bound update error!"); - case 4: - mexErrMsgTxt("Upper bound update error!"); - case 5: - mexErrMsgTxt("Matrices P and A update error!"); - case 6: - mexErrMsgTxt("Matrix P update error!"); - case 7: - mexErrMsgTxt("Matrix A update error!"); - } - - return; - } - - - // Warm start x and y variables - if (!strcmp("warm_start", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver has not been initialized."); - } - - // Fill x and y - const mxArray *x = prhs[2]; - const mxArray *y = prhs[3]; - - - // Copy vectors to ensure they are cast as c_float - c_float *x_vec; - c_float *y_vec; - - if(!mxIsEmpty(x)){ - x_vec = copyToCfloatVector(mxGetPr(x), - osqpData->work->data->n); - } - if(!mxIsEmpty(y)){ - y_vec = copyToCfloatVector(mxGetPr(y), - osqpData->work->data->m); - } - - // Warm start x and y - osqp_warm_start(osqpData->work, x_vec, y_vec); - - // Free vectors - if(!mxIsEmpty(x)) c_free(x_vec); - if(!mxIsEmpty(y)) c_free(y_vec); - - return; - } - - - // Warm start x variable - if (!strcmp("warm_start_x", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver has not been initialized."); - } - - // Fill x and y - const mxArray *x = prhs[2]; - - - // Copy vectors to ensure they are cast as c_float - c_float *x_vec; - - if(!mxIsEmpty(x)){ - x_vec = copyToCfloatVector(mxGetPr(x), - osqpData->work->data->n); - } - - // Warm start x - osqp_warm_start_x(osqpData->work, x_vec); - - // Free vectors - if(!mxIsEmpty(x)) c_free(x_vec); - - return; - } - - - // Warm start y variable - if (!strcmp("warm_start_y", cmd)) { - - //throw an error if this is called before solver is configured - if(!osqpData->work){ - mexErrMsgTxt("Solver has not been initialized."); - } - - // Fill x and y - const mxArray *y = prhs[2]; - - - // Copy vectors to ensure they are cast as c_float - c_float *y_vec; - - if(!mxIsEmpty(y)){ - y_vec = copyToCfloatVector(mxGetPr(y), - osqpData->work->data->m); - } - - // Warm start x - osqp_warm_start_y(osqpData->work, y_vec); - - // Free vectors - if(!mxIsEmpty(y)) c_free(y_vec); - - return; - } - - - - // SOLVE - if (!strcmp("solve", cmd)) { - if (nlhs != 5 || nrhs != 2){ - mexErrMsgTxt("Solve : wrong number of inputs / outputs"); - } - if(!osqpData->work){ - mexErrMsgTxt("No problem data has been given."); - } - // solve the problem - osqp_solve(osqpData->work); - - - // Allocate space for solution - // primal variables - plhs[0] = mxCreateDoubleMatrix(osqpData->work->data->n,1,mxREAL); - // dual variables - plhs[1] = mxCreateDoubleMatrix(osqpData->work->data->m,1,mxREAL); - // primal infeasibility certificate - plhs[2] = mxCreateDoubleMatrix(osqpData->work->data->m,1,mxREAL); - // dual infeasibility certificate - plhs[3] = mxCreateDoubleMatrix(osqpData->work->data->n,1,mxREAL); - - //copy results to mxArray outputs - //assume that five outputs will always - //be returned to matlab-side class wrapper - if ((osqpData->work->info->status_val != OSQP_PRIMAL_INFEASIBLE) && - (osqpData->work->info->status_val != OSQP_DUAL_INFEASIBLE)){ - - //primal and dual solutions - castToDoubleArr(osqpData->work->solution->x, mxGetPr(plhs[0]), osqpData->work->data->n); - castToDoubleArr(osqpData->work->solution->y, mxGetPr(plhs[1]), osqpData->work->data->m); - - //infeasibility certificates -> NaN values - setToNaN(mxGetPr(plhs[2]), osqpData->work->data->m); - setToNaN(mxGetPr(plhs[3]), osqpData->work->data->n); - - } else if (osqpData->work->info->status_val == OSQP_PRIMAL_INFEASIBLE || - osqpData->work->info->status_val == OSQP_PRIMAL_INFEASIBLE_INACCURATE){ //primal infeasible - - //primal and dual solutions -> NaN values - setToNaN(mxGetPr(plhs[0]), osqpData->work->data->n); - setToNaN(mxGetPr(plhs[1]), osqpData->work->data->m); - - //primal infeasibility certificates - castToDoubleArr(osqpData->work->delta_y, mxGetPr(plhs[2]), osqpData->work->data->m); - - //dual infeasibility certificates -> NaN values - setToNaN(mxGetPr(plhs[3]), osqpData->work->data->n); - - // Set objective value to infinity - osqpData->work->info->obj_val = mxGetInf(); - - } else { //dual infeasible - - //primal and dual solutions -> NaN values - setToNaN(mxGetPr(plhs[0]), osqpData->work->data->n); - setToNaN(mxGetPr(plhs[1]), osqpData->work->data->m); - - //primal infeasibility certificates -> NaN values - setToNaN(mxGetPr(plhs[2]), osqpData->work->data->m); - - //dual infeasibility certificates - castToDoubleArr(osqpData->work->delta_x, mxGetPr(plhs[3]), osqpData->work->data->n); - - // Set objective value to -infinity - osqpData->work->info->obj_val = -mxGetInf(); - } - - if (osqpData->work->info->status_val == OSQP_NON_CVX) { - osqpData->work->info->obj_val = mxGetNaN(); - } - - plhs[4] = copyInfoToMxStruct(osqpData->work->info); // Info structure - - return; - } - - if (!strcmp("constant", cmd)) { // Return solver constants - - char constant[32]; - mxGetString(prhs[2], constant, sizeof(constant)); - - if (!strcmp("OSQP_INFTY", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_INFTY); - return; - } - if (!strcmp("OSQP_NAN", constant)){ - plhs[0] = mxCreateDoubleScalar(mxGetNaN()); - return; - } - - if (!strcmp("OSQP_SOLVED", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_SOLVED); - return; - } - - if (!strcmp("OSQP_SOLVED_INACCURATE", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_SOLVED_INACCURATE); - return; - } - - if (!strcmp("OSQP_UNSOLVED", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_UNSOLVED); - return; - } - - if (!strcmp("OSQP_PRIMAL_INFEASIBLE", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_PRIMAL_INFEASIBLE); - return; - } - - if (!strcmp("OSQP_PRIMAL_INFEASIBLE_INACCURATE", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_PRIMAL_INFEASIBLE_INACCURATE); - return; - } - - if (!strcmp("OSQP_DUAL_INFEASIBLE", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_DUAL_INFEASIBLE); - return; - } - - if (!strcmp("OSQP_DUAL_INFEASIBLE_INACCURATE", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_DUAL_INFEASIBLE_INACCURATE); - return; - } - - if (!strcmp("OSQP_MAX_ITER_REACHED", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_MAX_ITER_REACHED); - return; - } - - if (!strcmp("OSQP_NON_CVX", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_NON_CVX); - return; - } - - if (!strcmp("OSQP_TIME_LIMIT_REACHED", constant)){ - plhs[0] = mxCreateDoubleScalar(OSQP_TIME_LIMIT_REACHED); - return; - } - - // Linear system solvers - if (!strcmp("QDLDL_SOLVER", constant)){ - plhs[0] = mxCreateDoubleScalar(QDLDL_SOLVER); - return; - } - - if (!strcmp("MKL_PARDISO_SOLVER", constant)){ - plhs[0] = mxCreateDoubleScalar(MKL_PARDISO_SOLVER); - return; - } - - - mexErrMsgTxt("Constant not recognized."); - - return; - } - - - // Got here, so command not recognized - mexErrMsgTxt("Command not recognized."); -} - - -c_float* copyToCfloatVector(double* vecData, c_int numel){ - // This memory needs to be freed! - c_float* out = (c_float*)c_malloc(numel * sizeof(c_float)); - - //copy data - for(c_int i=0; i < numel; i++){ - out[i] = (c_float)vecData[i]; - } - return out; - -} - -c_int* copyToCintVector(mwIndex* vecData, c_int numel){ - // This memory needs to be freed! - c_int* out = (c_int*)c_malloc(numel * sizeof(c_int)); - - //copy data - for(c_int i=0; i < numel; i++){ - out[i] = (c_int)vecData[i]; - } - return out; - -} - -c_int* copyDoubleToCintVector(double* vecData, c_int numel){ - // This memory needs to be freed! - c_int* out = (c_int*)c_malloc(numel * sizeof(c_int)); - - //copy data - for(c_int i=0; i < numel; i++){ - out[i] = (c_int)vecData[i]; - } - return out; - -} - -void castCintToDoubleArr(c_int *arr, double* arr_out, c_int len) { - for (c_int i = 0; i < len; i++) { - arr_out[i] = (double)arr[i]; - } -} - -void castToDoubleArr(c_float *arr, double* arr_out, c_int len) { - for (c_int i = 0; i < len; i++) { - arr_out[i] = (double)arr[i]; - } -} - -void setToNaN(double* arr_out, c_int len){ - c_int i; - for (i = 0; i < len; i++) { - arr_out[i] = mxGetNaN(); - } -} - - -mxArray* copyInfoToMxStruct(OSQPInfo* info){ - - //create mxArray with the right number of fields - int nfields = sizeof(OSQP_INFO_FIELDS) / sizeof(OSQP_INFO_FIELDS[0]); - mxArray* mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_INFO_FIELDS); - - //map the OSQP_INFO fields one at a time into mxArrays - //matlab all numeric values as doubles - mxSetField(mxPtr, 0, "iter", mxCreateDoubleScalar(info->iter)); - mxSetField(mxPtr, 0, "status", mxCreateString(info->status)); - mxSetField(mxPtr, 0, "status_val", mxCreateDoubleScalar(info->status_val)); - mxSetField(mxPtr, 0, "status_polish", mxCreateDoubleScalar(info->status_polish)); - mxSetField(mxPtr, 0, "obj_val", mxCreateDoubleScalar(info->obj_val)); - mxSetField(mxPtr, 0, "pri_res", mxCreateDoubleScalar(info->pri_res)); - mxSetField(mxPtr, 0, "dua_res", mxCreateDoubleScalar(info->dua_res)); - - #ifdef PROFILING - //if not profiling, these fields will be empty - mxSetField(mxPtr, 0, "setup_time", mxCreateDoubleScalar(info->setup_time)); - mxSetField(mxPtr, 0, "solve_time", mxCreateDoubleScalar(info->solve_time)); - mxSetField(mxPtr, 0, "update_time", mxCreateDoubleScalar(info->update_time)); - mxSetField(mxPtr, 0, "polish_time", mxCreateDoubleScalar(info->polish_time)); - mxSetField(mxPtr, 0, "run_time", mxCreateDoubleScalar(info->run_time)); - #endif - - mxSetField(mxPtr, 0, "rho_updates", mxCreateDoubleScalar(info->rho_updates)); - mxSetField(mxPtr, 0, "rho_estimate", mxCreateDoubleScalar(info->rho_estimate)); - - - return mxPtr; - -} - - -mxArray* copySettingsToMxStruct(OSQPSettings* settings){ - - int nfields = sizeof(OSQP_SETTINGS_FIELDS) / sizeof(OSQP_SETTINGS_FIELDS[0]); - mxArray* mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_SETTINGS_FIELDS); - - //map the OSQP_SETTINGS fields one at a time into mxArrays - //matlab handles everything as a double - mxSetField(mxPtr, 0, "rho", mxCreateDoubleScalar(settings->rho)); - mxSetField(mxPtr, 0, "sigma", mxCreateDoubleScalar(settings->sigma)); - mxSetField(mxPtr, 0, "scaling", mxCreateDoubleScalar(settings->scaling)); - mxSetField(mxPtr, 0, "adaptive_rho", mxCreateDoubleScalar(settings->adaptive_rho)); - mxSetField(mxPtr, 0, "adaptive_rho_interval", mxCreateDoubleScalar(settings->adaptive_rho_interval)); - mxSetField(mxPtr, 0, "adaptive_rho_tolerance", mxCreateDoubleScalar(settings->adaptive_rho_tolerance)); - mxSetField(mxPtr, 0, "adaptive_rho_fraction", mxCreateDoubleScalar(settings->adaptive_rho_fraction)); - mxSetField(mxPtr, 0, "max_iter", mxCreateDoubleScalar(settings->max_iter)); - mxSetField(mxPtr, 0, "eps_abs", mxCreateDoubleScalar(settings->eps_abs)); - mxSetField(mxPtr, 0, "eps_rel", mxCreateDoubleScalar(settings->eps_rel)); - mxSetField(mxPtr, 0, "eps_prim_inf", mxCreateDoubleScalar(settings->eps_prim_inf)); - mxSetField(mxPtr, 0, "eps_dual_inf", mxCreateDoubleScalar(settings->eps_dual_inf)); - mxSetField(mxPtr, 0, "alpha", mxCreateDoubleScalar(settings->alpha)); - mxSetField(mxPtr, 0, "linsys_solver", mxCreateDoubleScalar(settings->linsys_solver)); - mxSetField(mxPtr, 0, "delta", mxCreateDoubleScalar(settings->delta)); - mxSetField(mxPtr, 0, "polish", mxCreateDoubleScalar(settings->polish)); - mxSetField(mxPtr, 0, "polish_refine_iter", mxCreateDoubleScalar(settings->polish_refine_iter)); - mxSetField(mxPtr, 0, "verbose", mxCreateDoubleScalar(settings->verbose)); - mxSetField(mxPtr, 0, "scaled_termination", mxCreateDoubleScalar(settings->scaled_termination)); - mxSetField(mxPtr, 0, "check_termination", mxCreateDoubleScalar(settings->check_termination)); - mxSetField(mxPtr, 0, "warm_start", mxCreateDoubleScalar(settings->warm_start)); - mxSetField(mxPtr, 0, "time_limit", mxCreateDoubleScalar(settings->time_limit)); - - return mxPtr; -} - - -// ====================================================================== -mxArray* copyCscMatrixToMxStruct(csc* M){ - int nnzM; - int nfields = sizeof(CSC_FIELDS) / sizeof(CSC_FIELDS[0]); - mxArray* mxPtr = mxCreateStructMatrix(1,1,nfields,CSC_FIELDS); - - // Get number of nonzeros - nnzM = M->p[M->n]; - - // Create vectors - mxArray* p = mxCreateDoubleMatrix((M->n)+1,1,mxREAL); - mxArray* i = mxCreateDoubleMatrix(nnzM,1,mxREAL); - mxArray* x = mxCreateDoubleMatrix(nnzM,1,mxREAL); - - // Populate vectors - castCintToDoubleArr(M->p, mxGetPr(p), (M->n)+1); - castCintToDoubleArr(M->i, mxGetPr(i), nnzM); - castToDoubleArr(M->x, mxGetPr(x), nnzM); - - //map the CSC fields one at a time into mxArrays - //matlab handles everything as a double - mxSetField(mxPtr, 0, "nzmax", mxCreateDoubleScalar(M->nzmax)); - mxSetField(mxPtr, 0, "m", mxCreateDoubleScalar(M->m)); - mxSetField(mxPtr, 0, "n", mxCreateDoubleScalar(M->n)); - mxSetField(mxPtr, 0, "p", p); - mxSetField(mxPtr, 0, "i", i); - mxSetField(mxPtr, 0, "x", x); - mxSetField(mxPtr, 0, "nz", mxCreateDoubleScalar(M->nz)); - - return mxPtr; -} - -mxArray* copyDataToMxStruct(OSQPWorkspace* work){ - - int nfields = sizeof(OSQP_DATA_FIELDS) / sizeof(OSQP_DATA_FIELDS[0]); - mxArray* mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_DATA_FIELDS); - - // Create vectors - mxArray* q = mxCreateDoubleMatrix(work->data->n,1,mxREAL); - mxArray* l = mxCreateDoubleMatrix(work->data->m,1,mxREAL); - mxArray* u = mxCreateDoubleMatrix(work->data->m,1,mxREAL); - - // Populate vectors - castToDoubleArr(work->data->q, mxGetPr(q), work->data->n); - castToDoubleArr(work->data->l, mxGetPr(l), work->data->m); - castToDoubleArr(work->data->u, mxGetPr(u), work->data->m); - - // Create matrices - mxArray* P = copyCscMatrixToMxStruct(work->data->P); - mxArray* A = copyCscMatrixToMxStruct(work->data->A); - - //map the OSQP_DATA fields one at a time into mxArrays - //matlab handles everything as a double - mxSetField(mxPtr, 0, "n", mxCreateDoubleScalar(work->data->n)); - mxSetField(mxPtr, 0, "m", mxCreateDoubleScalar(work->data->m)); - mxSetField(mxPtr, 0, "P", P); - mxSetField(mxPtr, 0, "A", A); - mxSetField(mxPtr, 0, "q", q); - mxSetField(mxPtr, 0, "l", l); - mxSetField(mxPtr, 0, "u", u); - - return mxPtr; -} - -mxArray* copyLinsysSolverToMxStruct(OSQPWorkspace * work){ - - int nfields; - mxArray* mxPtr; - OSQPData * data; - qdldl_solver * linsys_solver; - - nfields = sizeof(LINSYS_SOLVER_FIELDS) / sizeof(LINSYS_SOLVER_FIELDS[0]); - mxPtr = mxCreateStructMatrix(1,1,nfields,LINSYS_SOLVER_FIELDS); - - data = work->data; - linsys_solver = (qdldl_solver *) work->linsys_solver; - - // Dimensions - int n = linsys_solver->L->n; - int Pdiag_n = linsys_solver->Pdiag_n; - int nnzP = data->P->p[data->P->n]; - int nnzA = data->A->p[data->A->n]; - - // Create vectors - mxArray* Dinv = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* P = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* bp = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* sol = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* rho_inv_vec = mxCreateDoubleMatrix(data->m,1,mxREAL); - mxArray* Pdiag_idx = mxCreateDoubleMatrix(Pdiag_n,1,mxREAL); - mxArray* PtoKKT = mxCreateDoubleMatrix(nnzP,1,mxREAL); - mxArray* AtoKKT = mxCreateDoubleMatrix(nnzA,1,mxREAL); - mxArray* rhotoKKT = mxCreateDoubleMatrix(data->m,1,mxREAL); - mxArray* D = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* etree = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* Lnz = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* iwork = mxCreateDoubleMatrix(3*n,1,mxREAL); - mxArray* bwork = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* fwork = mxCreateDoubleMatrix(n,1,mxREAL); - - // Populate vectors - castToDoubleArr(linsys_solver->Dinv, mxGetPr(Dinv), n); - castCintToDoubleArr(linsys_solver->P, mxGetPr(P), n); - castToDoubleArr(linsys_solver->bp, mxGetPr(bp), n); - castToDoubleArr(linsys_solver->sol, mxGetPr(sol), n); - castToDoubleArr(linsys_solver->rho_inv_vec, mxGetPr(rho_inv_vec), data->m); - castCintToDoubleArr(linsys_solver->Pdiag_idx, mxGetPr(Pdiag_idx), Pdiag_n); - castCintToDoubleArr(linsys_solver->PtoKKT, mxGetPr(PtoKKT), nnzP); - castCintToDoubleArr(linsys_solver->AtoKKT, mxGetPr(AtoKKT), nnzA); - castCintToDoubleArr(linsys_solver->rhotoKKT, mxGetPr(rhotoKKT), data->m); - castToDoubleArr(linsys_solver->D, mxGetPr(D), n); - castCintToDoubleArr(linsys_solver->etree, mxGetPr(etree), n); - castCintToDoubleArr(linsys_solver->Lnz, mxGetPr(Lnz), n); - castCintToDoubleArr(linsys_solver->iwork, mxGetPr(iwork), 3*n); - castCintToDoubleArr((c_int *)linsys_solver->bwork, mxGetPr(bwork), n); - castToDoubleArr(linsys_solver->fwork, mxGetPr(fwork), n); - - // Create matrices - mxArray* L = copyCscMatrixToMxStruct(linsys_solver->L); - mxArray* KKT = copyCscMatrixToMxStruct(linsys_solver->KKT); - - //map the PRIV fields one at a time into mxArrays - mxSetField(mxPtr, 0, "L", L); - mxSetField(mxPtr, 0, "Dinv", Dinv); - mxSetField(mxPtr, 0, "P", P); - mxSetField(mxPtr, 0, "bp", bp); - mxSetField(mxPtr, 0, "sol", sol); - mxSetField(mxPtr, 0, "rho_inv_vec", rho_inv_vec); - mxSetField(mxPtr, 0, "sigma", mxCreateDoubleScalar(linsys_solver->sigma)); - mxSetField(mxPtr, 0, "polish", mxCreateDoubleScalar(linsys_solver->polish)); - mxSetField(mxPtr, 0, "n", mxCreateDoubleScalar(linsys_solver->n)); - mxSetField(mxPtr, 0, "m", mxCreateDoubleScalar(linsys_solver->m)); - mxSetField(mxPtr, 0, "Pdiag_idx", Pdiag_idx); - mxSetField(mxPtr, 0, "Pdiag_n", mxCreateDoubleScalar(Pdiag_n)); - mxSetField(mxPtr, 0, "KKT", KKT); - mxSetField(mxPtr, 0, "PtoKKT", PtoKKT); - mxSetField(mxPtr, 0, "AtoKKT", AtoKKT); - mxSetField(mxPtr, 0, "rhotoKKT", rhotoKKT); - mxSetField(mxPtr, 0, "D", D); - mxSetField(mxPtr, 0, "etree", etree); - mxSetField(mxPtr, 0, "Lnz", Lnz); - mxSetField(mxPtr, 0, "iwork", iwork); - mxSetField(mxPtr, 0, "bwork", bwork); - mxSetField(mxPtr, 0, "fwork", fwork); - - return mxPtr; -} - -mxArray* copyScalingToMxStruct(OSQPWorkspace *work){ - - int n, m, nfields; - mxArray* mxPtr; - - - if (work->settings->scaling){ // Scaling performed - n = work->data->n; - m = work->data->m; - - nfields = sizeof(OSQP_SCALING_FIELDS) / sizeof(OSQP_SCALING_FIELDS[0]); - mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_SCALING_FIELDS); - - // Create vectors - mxArray* D = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* E = mxCreateDoubleMatrix(m,1,mxREAL); - mxArray* Dinv = mxCreateDoubleMatrix(n,1,mxREAL); - mxArray* Einv = mxCreateDoubleMatrix(m,1,mxREAL); - - // Populate vectors - castToDoubleArr(work->scaling->D, mxGetPr(D), n); - castToDoubleArr(work->scaling->E, mxGetPr(E), m); - castToDoubleArr(work->scaling->Dinv, mxGetPr(Dinv), n); - castToDoubleArr(work->scaling->Einv, mxGetPr(Einv), m); - - //map the SCALING fields one at a time - mxSetField(mxPtr, 0, "c", mxCreateDoubleScalar(work->scaling->c)); - mxSetField(mxPtr, 0, "D", D); - mxSetField(mxPtr, 0, "E", E); - mxSetField(mxPtr, 0, "cinv", mxCreateDoubleScalar(work->scaling->cinv)); - mxSetField(mxPtr, 0, "Dinv", Dinv); - mxSetField(mxPtr, 0, "Einv", Einv); - - } else { - mxPtr = mxCreateDoubleMatrix(0, 0, mxREAL); - } - - return mxPtr; -} - -mxArray* copyRhoVectorsToMxStruct(OSQPWorkspace *work){ - - int m, nfields; - mxArray* mxPtr; - - m = work->data->m; - - nfields = sizeof(OSQP_RHO_VECTORS_FIELDS) / sizeof(OSQP_RHO_VECTORS_FIELDS[0]); - mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_RHO_VECTORS_FIELDS); - - // Create vectors - mxArray* rho_vec = mxCreateDoubleMatrix(m,1,mxREAL); - mxArray* rho_inv_vec = mxCreateDoubleMatrix(m,1,mxREAL); - mxArray* constr_type = mxCreateDoubleMatrix(m,1,mxREAL); - - // Populate vectors - castToDoubleArr(work->rho_vec, mxGetPr(rho_vec), m); - castToDoubleArr(work->rho_inv_vec, mxGetPr(rho_inv_vec), m); - castCintToDoubleArr(work->constr_type, mxGetPr(constr_type), m); - - //map the RHO_VECTORS fields one at a time into mxArrays - mxSetField(mxPtr, 0, "rho_vec", rho_vec); - mxSetField(mxPtr, 0, "rho_inv_vec", rho_inv_vec); - mxSetField(mxPtr, 0, "constr_type", constr_type); - - return mxPtr; -} - - -mxArray* copyWorkToMxStruct(OSQPWorkspace* work){ - - int nfields = sizeof(OSQP_WORKSPACE_FIELDS) / sizeof(OSQP_WORKSPACE_FIELDS[0]); - mxArray* mxPtr = mxCreateStructMatrix(1,1,nfields,OSQP_WORKSPACE_FIELDS); - - // Create workspace substructures - mxArray* rho_vectors = copyRhoVectorsToMxStruct(work); - mxArray* data = copyDataToMxStruct(work); - mxArray* linsys_solver = copyLinsysSolverToMxStruct(work); - mxArray* scaling = copyScalingToMxStruct(work); - mxArray* settings = copySettingsToMxStruct(work->settings); - - //map the WORKSPACE fields one at a time into mxArrays - mxSetField(mxPtr, 0, "rho_vectors", rho_vectors); - mxSetField(mxPtr, 0, "data", data); - mxSetField(mxPtr, 0, "linsys_solver", linsys_solver); - mxSetField(mxPtr, 0, "scaling", scaling); - mxSetField(mxPtr, 0, "settings", settings); - - return mxPtr; -} - -// ====================================================================== - - -void copyMxStructToSettings(const mxArray* mxPtr, OSQPSettings* settings){ - - //this function assumes that only a complete and validated structure - //will be passed. matlab mex-side function is responsible for checking - //structure validity - - //map the OSQP_SETTINGS fields one at a time into mxArrays - //matlab handles everything as a double - settings->rho = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "rho")); - settings->sigma = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "sigma")); - settings->scaling = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "scaling")); - settings->adaptive_rho = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho")); - settings->adaptive_rho_interval = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho_interval")); - settings->adaptive_rho_tolerance = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho_tolerance")); - settings->adaptive_rho_fraction = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "adaptive_rho_fraction")); - settings->max_iter = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "max_iter")); - settings->eps_abs = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_abs")); - settings->eps_rel = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_rel")); - settings->eps_prim_inf = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_dual_inf")); - settings->eps_dual_inf = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_dual_inf")); - settings->alpha = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "alpha")); - settings->linsys_solver = (enum linsys_solver_type) (c_int) mxGetScalar(mxGetField(mxPtr, 0, "linsys_solver")); - settings->delta = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "delta")); - settings->polish = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "polish")); - settings->polish_refine_iter = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "polish_refine_iter")); - settings->verbose = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "verbose")); - settings->scaled_termination = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "scaled_termination")); - settings->check_termination = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "check_termination")); - settings->warm_start = (c_int)mxGetScalar(mxGetField(mxPtr, 0, "warm_start")); - settings->time_limit = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "time_limit")); - -} - -void copyUpdatedSettingsToWork(const mxArray* mxPtr ,OsqpData* osqpData){ - - c_int exitflag; - - //This does basically the same job as copyMxStructToSettings which was used - //during setup, but uses the provided update functions in osqp.h to update - //settings in the osqp workspace. Protects against bad parameter writes - //or future modifications to updated settings handling - - osqp_update_max_iter(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "max_iter"))); - osqp_update_eps_abs(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_abs"))); - osqp_update_eps_rel(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_rel"))); - osqp_update_eps_prim_inf(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_prim_inf"))); - osqp_update_eps_dual_inf(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "eps_dual_inf"))); - osqp_update_alpha(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "alpha"))); - osqp_update_delta(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "delta"))); - osqp_update_polish(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "polish"))); - osqp_update_polish_refine_iter(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "polish_refine_iter"))); - osqp_update_verbose(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "verbose"))); - osqp_update_scaled_termination(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "scaled_termination"))); - osqp_update_check_termination(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "check_termination"))); - osqp_update_warm_start(osqpData->work, - (c_int)mxGetScalar(mxGetField(mxPtr, 0, "warm_start"))); - osqp_update_time_limit(osqpData->work, - (c_float)mxGetScalar(mxGetField(mxPtr, 0, "time_limit"))); - - - // Check for settings that need special update - // Update them only if they are different than already set values - c_float rho_new = (c_float)mxGetScalar(mxGetField(mxPtr, 0, "rho")); - // Check if it has changed - if (c_absval(rho_new - osqpData->work->settings->rho) > NEW_SETTINGS_TOL){ - exitflag = osqp_update_rho(osqpData->work, rho_new); - if (exitflag){ - mexErrMsgTxt("rho update error!"); - } - } - - -} diff --git a/osqp_sources b/osqp_sources deleted file mode 160000 index f9fc23d..0000000 --- a/osqp_sources +++ /dev/null @@ -1 +0,0 @@ -Subproject commit f9fc23d3436e4b17dd2cb95f70cfa1f37d122c24 diff --git a/package/install_osqp.m b/package/install_osqp.m index 7df3add..4fc6109 100644 --- a/package/install_osqp.m +++ b/package/install_osqp.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + function install_osqp % Install the OSQP solver Matlab interface diff --git a/package/package_osqp.m b/package/package_osqp.m index b0c39b2..cd730f3 100644 --- a/package/package_osqp.m +++ b/package/package_osqp.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + function package_osqp(version) % Create OSQP matlab interface package diff --git a/run_osqp_tests.m b/run_osqp_tests.m index 29d9d18..da1973a 100644 --- a/run_osqp_tests.m +++ b/run_osqp_tests.m @@ -1,11 +1,19 @@ +% SPDX-License-Identifier: Apache-2.0 + import matlab.unittest.TestSuite; -[osqp_path,~,~] = fileparts(which('osqp.m')); -unittest_dir = fullfile(osqp_path, 'unittests'); +[osqp_classpath,~,~] = fileparts( mfilename( 'fullpath' ) ); +unittest_dir = fullfile(osqp_classpath, 'unittests'); suiteFolder = TestSuite.fromFolder(unittest_dir); % Solve individual test file %suiteFolder = TestSuite.fromFile('unittests/dual_infeasibility_tests.m'); -% Run all suite -result = run(suiteFolder); +% Run all tests +results = run(suiteFolder); + +if all([results.Failed] == 0) + exit(0); +else + exit(1); +end \ No newline at end of file diff --git a/unittests/basic_tests.m b/unittests/basic_tests.m index 269efe0..4d21ad1 100644 --- a/unittests/basic_tests.m +++ b/unittests/basic_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef basic_tests < matlab.unittest.TestCase %TEST_BASIC_QP Solve Basic QP Problem @@ -124,6 +126,9 @@ function test_update_max_iter(testCase) opts.max_iter = 30; testCase.solver.update_settings(opts); + set = testCase.solver.current_settings(); + testCase.verifyEqual(set.max_iter, 30) + % Solve again results = testCase.solver.solve(); @@ -139,6 +144,9 @@ function test_update_early_termination(testCase) opts.check_termination = 0; testCase.solver.update_settings(opts); + set = testCase.solver.current_settings(); + testCase.verifyEqual(set.check_termination, 0) + % Solve again results = testCase.solver.solve(); @@ -170,7 +178,7 @@ function test_update_rho(testCase) end function test_update_time_limit(testCase) - testCase.verifyEqual(testCase.options.time_limit, 0) + testCase.verifyEqual(testCase.options.time_limit, 1e10) results = testCase.solver.solve(); testCase.verifyEqual(results.info.status_val, ... @@ -184,6 +192,11 @@ function test_update_time_limit(testCase) 'max_iter', 2e9,... 'check_termination', 0); + set = testCase.solver.current_settings(); + testCase.verifyEqual(set.check_termination, 0) + testCase.verifyEqual(set.time_limit, 1e-6) + testCase.verifyEqual(set.max_iter, 2e9) + results = testCase.solver.solve(); testCase.verifyEqual(results.info.status_val, ... testCase.solver.constant('OSQP_TIME_LIMIT_REACHED')) diff --git a/unittests/codegen_mat_tests.m b/unittests/codegen_mat_tests.m index 31be98d..5087676 100644 --- a/unittests/codegen_mat_tests.m +++ b/unittests/codegen_mat_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef codegen_mat_tests < matlab.unittest.TestCase %TEST_BASIC_QP Solve Basic QP Problem @@ -17,6 +19,7 @@ methods(TestMethodSetup) function setup_problem(testCase) % Create Problem + assumeFail(testCase); testCase.P = sparse([11 0; 0, 0.1]); testCase.q = [3; 4]; testCase.A = sparse([-1. 0; 0 -1; -1 -3; 2 5; 3 4]); diff --git a/unittests/codegen_vec_tests.m b/unittests/codegen_vec_tests.m index 972c11f..d9c6843 100644 --- a/unittests/codegen_vec_tests.m +++ b/unittests/codegen_vec_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef codegen_vec_tests < matlab.unittest.TestCase %TEST_BASIC_QP Solve Basic QP Problem @@ -17,6 +19,7 @@ methods(TestMethodSetup) function setup_problem(testCase) % Create Problem + assumeFail(testCase); testCase.P = sparse([11 0; 0, 0]); testCase.q = [3; 4]; testCase.A = sparse([-1. 0; 0 -1; -1 -3; 2 5; 3 4]); diff --git a/unittests/dual_infeasibility_tests.m b/unittests/dual_infeasibility_tests.m index edeb680..3cc861c 100644 --- a/unittests/dual_infeasibility_tests.m +++ b/unittests/dual_infeasibility_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef dual_infeasibility_tests < matlab.unittest.TestCase properties diff --git a/unittests/feasibility_tests.m b/unittests/feasibility_tests.m index cc15f6d..b2e2801 100644 --- a/unittests/feasibility_tests.m +++ b/unittests/feasibility_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef feasibility_tests < matlab.unittest.TestCase %FEASIBILITY_TESTS Solve equality constrained feasibility problem diff --git a/unittests/non_cvx_tests.m b/unittests/non_cvx_tests.m index 53df8e9..ab9eab2 100644 --- a/unittests/non_cvx_tests.m +++ b/unittests/non_cvx_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef non_cvx_tests < matlab.unittest.TestCase %NON_CVX_TESTS Try to solve a non-convex QP diff --git a/unittests/primal_infeasibility_tests.m b/unittests/primal_infeasibility_tests.m index fdea926..a9568a8 100644 --- a/unittests/primal_infeasibility_tests.m +++ b/unittests/primal_infeasibility_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef primal_infeasibility_tests < matlab.unittest.TestCase properties diff --git a/unittests/unconstrained_tests.m b/unittests/unconstrained_tests.m index fd4c1c1..a92b1f5 100644 --- a/unittests/unconstrained_tests.m +++ b/unittests/unconstrained_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef unconstrained_tests < matlab.unittest.TestCase %UNCONSTRAINED_TESTS Solve unconstrained quadratic program diff --git a/unittests/update_matrices_tests.m b/unittests/update_matrices_tests.m index 8948822..e52e615 100644 --- a/unittests/update_matrices_tests.m +++ b/unittests/update_matrices_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef update_matrices_tests < matlab.unittest.TestCase %TEST_BASIC_QP Solve Basic QP Problem @@ -41,7 +43,7 @@ function setup_problem(testCase) % Setup solver testCase.solver = osqp; testCase.solver.setup(testCase.P, testCase.q, testCase.A, ... - testCase.l, testCase.u, 'verbose', false, 'eps_rel', 1e-7, 'eps_abs', 1e-07, 'polish', true); + testCase.l, testCase.u, 'verbose', false, 'eps_rel', 1e-7, 'eps_abs', 1e-07, 'polishing', true); % Setup tolerance testCase.tol = 1e-04; diff --git a/unittests/warm_start_tests.m b/unittests/warm_start_tests.m index f6ab098..2b33375 100644 --- a/unittests/warm_start_tests.m +++ b/unittests/warm_start_tests.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + classdef warm_start_tests < matlab.unittest.TestCase %WARM_START_TESTS Warm Start problems solution @@ -33,8 +35,7 @@ function setup_problem(testCase) end methods (Test) - function test_warm_start(testCase) - + function test_warm_start_zeros(testCase) % big example rng(4) testCase.n = 100; @@ -63,12 +64,99 @@ function test_warm_start(testCase) testCase.solver.warm_start('x', zeros(testCase.n, 1), 'y', zeros(testCase.m, 1)); results = testCase.solver.solve(); testCase.verifyEqual(results.info.iter, tot_iter, 'AbsTol', testCase.tol) + end + + function test_warm_start_optimal(testCase) + % big example + rng(4) + testCase.n = 100; + testCase.m = 200; + Pt = sprandn(testCase.n, testCase.n, 0.6); + testCase.P = Pt' * Pt; + testCase.q = randn(testCase.n, 1); + testCase.A = sprandn(testCase.m, testCase.n, 0.8); + testCase.u = 2*rand(testCase.m, 1); + testCase.l = -2*rand(testCase.m, 1); + + % Setup solver + testCase.solver = osqp; + testCase.solver.setup(testCase.P, testCase.q, ... + testCase.A, testCase.l, testCase.u, testCase.options); + + % Solve with OSQP + results = testCase.solver.solve(); + + % Store optimal values + x_opt = results.x; + y_opt = results.y; + tot_iter = results.info.iter; % Warm start with optimal values and check that number of iterations is < 10 testCase.solver.warm_start('x', x_opt, 'y', y_opt); results = testCase.solver.solve(); testCase.verifyThat(results.info.iter, matlab.unittest.constraints.IsLessThan(10)); + end + + function test_warm_start_duals(testCase) + % big example + rng(4) + testCase.n = 100; + testCase.m = 200; + Pt = sprandn(testCase.n, testCase.n, 0.6); + testCase.P = Pt' * Pt; + testCase.q = randn(testCase.n, 1); + testCase.A = sprandn(testCase.m, testCase.n, 0.8); + testCase.u = 2*rand(testCase.m, 1); + testCase.l = -2*rand(testCase.m, 1); + + % Setup solver + testCase.solver = osqp; + testCase.solver.setup(testCase.P, testCase.q, ... + testCase.A, testCase.l, testCase.u, testCase.options); + + % Solve with OSQP + results = testCase.solver.solve(); + + % Store optimal values + x_opt = results.x; + y_opt = results.y; + tot_iter = results.info.iter; + + % Warm start with zeros for dual variables + testCase.solver.warm_start('y', zeros(testCase.m, 1)); + results = testCase.solver.solve(); + testCase.verifyEqual(results.y, y_opt, 'AbsTol', testCase.tol) + end + + function test_warm_start_primal(testCase) + % big example + rng(4) + testCase.n = 100; + testCase.m = 200; + Pt = sprandn(testCase.n, testCase.n, 0.6); + testCase.P = Pt' * Pt; + testCase.q = randn(testCase.n, 1); + testCase.A = sprandn(testCase.m, testCase.n, 0.8); + testCase.u = 2*rand(testCase.m, 1); + testCase.l = -2*rand(testCase.m, 1); + + % Setup solver + testCase.solver = osqp; + testCase.solver.setup(testCase.P, testCase.q, ... + testCase.A, testCase.l, testCase.u, testCase.options); + + % Solve with OSQP + results = testCase.solver.solve(); + + % Store optimal values + x_opt = results.x; + y_opt = results.y; + tot_iter = results.info.iter; + % Warm start with zeros for primal variables + testCase.solver.warm_start('x', zeros(testCase.n, 1)); + results = testCase.solver.solve(); + testCase.verifyEqual(results.x, x_opt, 'AbsTol', testCase.tol) end end diff --git a/utils/convertProblemToMat.m b/utils/convertProblemToMat.m index 614e6c0..cc9eece 100644 --- a/utils/convertProblemToMat.m +++ b/utils/convertProblemToMat.m @@ -1,3 +1,5 @@ +% SPDX-License-Identifier: Apache-2.0 + function convertProblemToMat(filename, Pdata, qdata, Adata, ldata, udata) Ptrip = load(Pdata); P = spconvert(Ptrip);