Added Matlab interface for sparse linear algebra.

SQProblemSchur now works with Matlab's shipped MA57 (libmwma57, libmwlapack, libmwblas).

.gitattributes

@@ -85,6 +85,8 @@ interfaces/matlab/qpOASES_matlab_utils.hpp -text
 interfaces/matlab/qpOASES_options.m -text
 interfaces/matlab/qpOASES_sequence.cpp -text
 interfaces/matlab/qpOASES_sequence.m -text
+interfaces/matlab/testQPset.m -text
+interfaces/matlab/testSchur.m -text
 interfaces/octave/clean -text
 interfaces/octave/clean.sh -text
 interfaces/octave/example1.mat -text
@@ -187,6 +189,7 @@ testing/cpp/test_qrecipeSchur.cpp -text
 testing/cpp/test_qrecipe_data.hpp -text
 testing/cpp/test_runAllOqpExamples.cpp -text
 testing/cpp/test_sebastien1.cpp -text
+testing/cpp/test_smallSchur.cpp -text
 testing/cpp/test_vanBarelsUnboundedQP.cpp -text
 testing/matlab/auxFiles/generateExample.m -text
 testing/matlab/auxFiles/generateRandomQp.m -text

include/qpOASES.hpp

@@ -38,8 +38,10 @@
 #include <SubjectTo.cpp>
 #include <Bounds.cpp>
 #include <Constraints.cpp>
+#ifndef __MATLAB__
 #include <BLASReplacement.cpp>
 #include <LAPACKReplacement.cpp>
+#endif
 #include <Matrices.cpp>
 #include <Options.cpp>
 #include <QProblemB.cpp>

include/qpOASES/Types.hpp

@@ -174,7 +174,11 @@ typedef unsigned int uint_t;
 
 
 /** typedef for Fortran INTEGER type. Might be platform dependent! */
+#ifdef __USE_LONG_FINTS__
+typedef long fint;
+#else
 typedef int fint;
+#endif /* __USE_LONG_FINTS__ */
 
 
 /**

interfaces/matlab/make.m

@@ -62,10 +62,10 @@
     QPOASESPATH = '../../';
 
     DEBUGFLAGS = ' ';
-    %DEBUGFLAGS = ' -g CXXDEBUGFLAGS=''$CXXDEBUGFLAGS -Wall -pedantic -Wshadow'' ';
+    %DEBUGFLAGS = ' -v -g CXXDEBUGFLAGS=''$CXXDEBUGFLAGS -Wall -pedantic -Wshadow'' ';
 
     IFLAGS = [ '-I. -I',QPOASESPATH,'include',' -I',QPOASESPATH,'src',' ' ];
-    CPPFLAGS = [ IFLAGS, DEBUGFLAGS, '-largeArrayDims -D__cpluplus -D__MATLAB__ -D__SINGLE_OBJECT__',' ' ];
+    CPPFLAGS = [ IFLAGS, DEBUGFLAGS, '-largeArrayDims -D__cpluplus -D__MATLAB__ -D__SINGLE_OBJECT__ -lmwblas',' ' ];
     defaultFlags = '-O -D__NO_COPYRIGHT__ '; %% -D__SUPPRESSANYOUTPUT__
 
     if ( ispc == 0 )
@@ -80,6 +80,18 @@
         CPPFLAGS = [ CPPFLAGS, userFlags,' ' ];
     end
 
+	%% determine if MA57 is available for sparse linear algebra
+	isoctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
+	if isoctave
+		warning('Sparse linear algebra is currently not available for qpOASES in Octave. Passing sparse matrices works but will likely be slow.')
+		SPARSEFLAGS = '';
+	elseif verLessThan('matlab', '7.8')
+		warning('Sparse linear algebra is currently available for qpOASES only for Matlab versions 7.8 and later. Passing sparse matrices works but will likely be slow.')
+		SPARSEFLAGS = '';
+	else
+		SPARSEFLAGS = '-largeArrayDims -D__USE_LONG_INTEGERS__ -D__USE_LONG_FINTS__ -DWITH_SPARSE_LA -DSOLVER_MA57 -lmwma57 -lmwlapack ';
+	end
+
     mexExt = eval('mexext');
 
 
@@ -108,7 +120,7 @@
     %% call mex compiler
     for ii=1:length(fcnNames)
 
-        cmd = [ 'mex -output ', fcnNames{ii}, ' ', CPPFLAGS, [fcnNames{ii},'.cpp'] ];
+        cmd = [ 'mex -output ', fcnNames{ii}, ' ', CPPFLAGS, SPARSEFLAGS, [fcnNames{ii},'.cpp'] ];
 
         if ( exist( [fcnNames{ii},'.',mexExt],'file' ) == 0 )
 

interfaces/matlab/qpOASES.cpp

@@ -60,15 +60,20 @@ int_t QProblem_qpOASES(	int_t nV, int_t nC, HessianType hessianType, int_t nP,
 						const double* const x0, Options* options,
 						int_t nOutputs, mxArray* plhs[],
 						const double* const guessedBounds, const double* const guessedConstraints,
-						const double* const _R
+						const double* const _R,
+						BooleanType isSparse
 						)
 {
 	int_t nWSRout;
 	real_t maxCpuTimeOut;
 
 	/* 1) Setup initial QP. */
-	QProblem QP( nV,nC,hessianType );
-	QP.setOptions( *options );
+	QProblem *QP;
+	if ( isSparse )
+		QP = new SQProblemSchur ( nV,nC,hessianType );
+	else
+		QP = new QProblem ( nV,nC,hessianType );
+	QP->setOptions( *options );
 
 	/* 2) Solve initial QP. */
 	returnValue returnvalue;
@@ -114,7 +119,7 @@ int_t QProblem_qpOASES(	int_t nV, int_t nC, HessianType hessianType, int_t nP,
 	nWSRout = nWSRin;
 	maxCpuTimeOut = (maxCpuTimeIn >= 0.0) ? maxCpuTimeIn : INFTY;
 
-	returnvalue = QP.init(	H,g,A,lb,ub,lbA,ubA,
+	returnvalue = QP->init(	H,g,A,lb,ub,lbA,ubA,
 							nWSRout,&maxCpuTimeOut,
 							x0,0,
 							(guessedBounds != 0) ? &bounds : 0, (guessedConstraints != 0) ? &constraints : 0,
@@ -147,16 +152,17 @@ int_t QProblem_qpOASES(	int_t nV, int_t nC, HessianType hessianType, int_t nP,
 
 			nWSRout = nWSRin;
 			maxCpuTimeOut = (maxCpuTimeIn >= 0.0) ? maxCpuTimeIn : INFTY;
-			returnvalue = QP.hotstart( g_current,lb_current,ub_current,lbA_current,ubA_current, nWSRout,&maxCpuTimeOut );
+			returnvalue = QP->hotstart( g_current,lb_current,ub_current,lbA_current,ubA_current, nWSRout,&maxCpuTimeOut );
 		}
 
 		/* write results into output vectors */
-		obtainOutputs(	k,&QP,returnvalue,nWSRout,maxCpuTimeOut,
+		obtainOutputs(	k,QP,returnvalue,nWSRout,maxCpuTimeOut,
 						nOutputs,plhs,nV,nC );
 	}
 
-	//QP.writeQpDataIntoMatFile( "qpDataMat0.mat" );
+	//QP->writeQpDataIntoMatFile( "qpDataMat0.mat" );
 
+	delete QP;
 	return 0;
 }
 
@@ -277,6 +283,11 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
 	/* dimensions */
 	uint_t nV=0, nC=0, nP=0;
 	BooleanType isSimplyBoundedQp = BT_FALSE;
+	#ifdef WITH_SPARSE_LA
+	BooleanType isSparse = BT_TRUE; // This will be set to BT_FALSE later if a dense matrix is encountered.
+	#else
+	BooleanType isSparse = BT_FALSE;
+	#endif
 
 	/* sparse matrix indices and values */
 	sparse_int_t *Hir=0, *Hjc=0, *Air=0, *Ajc=0;
@@ -529,6 +540,12 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
 
 	allocateOutputs( nlhs,plhs,nV,nC,nP );
 
+	/* check if QP is sparse */
+	if ( H_idx >= 0 && !mxIsSparse( prhs[H_idx] ) )
+		isSparse = BT_FALSE;
+	if ( nC > 0 && A_idx >= 0 && !mxIsSparse( prhs[A_idx] ) )
+		isSparse = BT_FALSE;
+
 	if ( nC == 0 )
 	{
 		/* Call qpOASES (using QProblemB class). */
@@ -563,7 +580,8 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
 							nWSRin,maxCpuTimeIn,
 							x0,&options,
 							nlhs,plhs,
-							guessedBounds,guessedConstraints,R
+							guessedBounds,guessedConstraints,R,
+							isSparse
 							);
 
 		if (R != 0) delete R;

interfaces/matlab/qpOASES_matlab_utils.cpp

@@ -36,8 +36,9 @@
 
 
 QPInstance::QPInstance(	uint_t _nV, uint_t _nC, HessianType _hessianType,
-						BooleanType _isSimplyBounded
+						BooleanType _isSimplyBounded, BooleanType _sparseLA
 						)
+	: sparseLA(_sparseLA)
 {
 	handle = s_nexthandle++;
 
@@ -46,16 +47,21 @@ QPInstance::QPInstance(	uint_t _nV, uint_t _nC, HessianType _hessianType,
 	else
 		isSimplyBounded = _isSimplyBounded;
 
-	if ( isSimplyBounded == BT_TRUE )
+	if ( isSimplyBounded == BT_TRUE && sparseLA == BT_FALSE )
 	{
 		sqp = 0;
 		qpb = new QProblemB( _nV,_hessianType );
 	}
-	else
+	else if ( sparseLA == BT_FALSE )
 	{
 		sqp = new SQProblem( _nV,_nC,_hessianType );
 		qpb = 0;
 	}
+	else
+	{
+		sqp = new SQProblemSchur( _nV,_nC,_hessianType );
+		qpb = 0;
+	}
 
 	H = 0;
 	A = 0;
@@ -179,10 +185,10 @@ bool mxIsScalar( const mxArray *pm )
  *	a l l o c a t e Q P r o b l e m I n s t a n c e
  */
 int_t allocateQPInstance(	int_t nV, int_t nC, HessianType hessianType,
-							BooleanType isSimplyBounded, const Options* options
+							BooleanType isSimplyBounded, BooleanType isSparse, const Options* options
 							)
 {
-	QPInstance* inst = new QPInstance( nV,nC,hessianType, isSimplyBounded );
+	QPInstance* inst = new QPInstance( nV,nC,hessianType, isSimplyBounded, isSparse );
 
 	if ( ( inst->sqp != 0 ) && ( options != 0 ) )
 		inst->sqp->setOptions( *options );

interfaces/matlab/qpOASES_matlab_utils.hpp

@@ -61,7 +61,8 @@ class QPInstance
 		QPInstance(	uint_t _nV = 0,
 					uint_t _nC = 0,
 					HessianType _hessianType = HST_UNKNOWN,
-					BooleanType _isSimplyBounded = BT_FALSE
+					BooleanType _isSimplyBounded = BT_FALSE,
+					BooleanType _sparseLA = BT_FALSE
 					);
 
 		~QPInstance( );
@@ -76,6 +77,7 @@ class QPInstance
 		SQProblem* sqp;
 		QProblemB* qpb;
 		BooleanType isSimplyBounded;
+		BooleanType sparseLA;
 
 		SymmetricMatrix* H;
 		Matrix* A;

interfaces/matlab/qpOASES_sequence.cpp

@@ -347,6 +347,11 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
 	int_t x0_idx=-1, auxInput_idx=-1;
 
 	BooleanType isSimplyBoundedQp = BT_FALSE;
+	#ifdef WITH_SPARSE_LA
+	BooleanType isSparse = BT_TRUE; // This will be set to BT_FALSE later if a dense matrix is encountered.
+	#else
+	BooleanType isSparse = BT_FALSE;
+	#endif
 
 	Options options;
 	options.printLevel = PL_LOW;
@@ -609,9 +614,15 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] )
 			R = new real_t[nV*nV];
 			convertFortranToC( R_for, nV,nV, R );
 		}
+
+		/* check if QP is sparse */
+		if ( H_idx >= 0 && !mxIsSparse( prhs[H_idx] ) )
+			isSparse = BT_FALSE;
+		if ( nC > 0 && A_idx >= 0 && !mxIsSparse( prhs[A_idx] ) )
+			isSparse = BT_FALSE;
 
 		/* allocate instance */
-		handle = allocateQPInstance( nV,nC,hessianType, isSimplyBoundedQp,&options );	
+		handle = allocateQPInstance( nV,nC,hessianType, isSimplyBoundedQp, isSparse, &options );	
 		globalQP = getQPInstance( handle );
 
 		/* make a deep-copy of the user-specified Hessian matrix (possibly sparse) */

interfaces/matlab/testQPset.m

@@ -0,0 +1,43 @@
+% Requires Lukas Schork's Matlab repository of the Maros Meszaros test set
+QPsetDIR = '~/git/QPset/maros';
+
+files = dir(QPsetDIR);
+
+options = qpOASES_options('default', ...
+	'printLevel', 0, 'maxIter', 1e8);%, 'initialStatusBounds', 0);
+nQP = 0;
+for i = 1:length(files)
+	clear mex
+	try
+		load([QPsetDIR, '/', files(i).name]);
+		H = sparse(H);
+		A = sparse(A);
+	catch
+		continue
+	end
+
+	if max(size(A)) > 2000
+		continue
+	end
+
+	startTime = tic;
+	try
+		[x,fval,exitflag,iter,lambda] = qpOASES(H,g,A,xl,xu,al,au,options);
+	catch
+		exitflag = 666;
+		iter = 666;
+		fval = 666;
+	end
+	elapsedTime = toc(startTime);
+
+	if mod(nQP, 20) == 0
+		fprintf('\n%-10s %6s %6s %5s %7s %15s %9s\n', ...
+			'name', 'nvar', 'ncon', 'eflag', 'iter', 'fval', 'time [s]');
+	end
+	fprintf('%-10s %6d %6d %5d %7d %15g %9.3e\n', ...
+		files(i).name(1:findstr(files(i).name, '.')-1), ...
+		size(A,2), size(A,1), exitflag, iter, fval, elapsedTime);
+
+	nQP = nQP + 1;
+end
+

interfaces/matlab/testSchur.m

@@ -0,0 +1,16 @@
+m = 2000;
+n = 2*m+1;
+H = 2*speye(n);
+A = 0.5*H(1:m,:);
+g = zeros(n,1);
+lb = -ones(n,1);
+ub = -lb;
+lbA = 0.5*ones(m,1);
+ubA = 0.5*ones(m,1);
+
+options = qpOASES_options('default', 'printLevel', -1, 'maxIter', n+m+1, ...
+	'enableEqualities', 0, 'initialStatusBounds', 1);
+
+%[x,fval,exitflag,iter,lambda] = qpOASES(full(H),g,full(A),lb,ub,lbA,ubA,options);
+[x,fval,exitflag,iter,lambda] = qpOASES(H,g,A,lb,ub,lbA,ubA,options);
+

make_linux.mk

@@ -49,23 +49,29 @@ ifeq ($(REPLACE_LINALG), 1)
 	LIB_BLAS =   ${SRCDIR}/BLASReplacement.o
 	LIB_LAPACK = ${SRCDIR}/LAPACKReplacement.o
 else
-	LIB_BLAS =   /usr/lib/libblas.so
-	LIB_LAPACK = /usr/lib/liblapack.so
+	LIB_BLAS =   /usr/lib/libblas.so.3gf
+	LIB_LAPACK = /usr/lib/liblapack.so.3gf
+#	LIB_BLAS = ${MATLAB_LIBDIR}/libmwblas.so
+#	LIB_LAPACK = ${MATLAB_LIBDIR}/libmwlapack.so
 endif
 
 # choice of sparse solver: NONE, MA27, or MA57
 # If choice is not 'NONE', BLAS and LAPACK replacements must not be used
 USE_SOLVER = NONE
+#USE_SOLVER = MA57
 
 ifeq ($(USE_SOLVER), MA57)
-	LIB_SOLVER = /usr/local/lib/libhsl_ma57.a /usr/local/lib/libfakemetis.a
+	LIB_SOLVER = ${MATLAB_LIBDIR}/libmwma57.so
 	DEF_SOLVER = SOLVER_MA57
+	LINKHSL = -Wl,-rpath=${MATLAB_LIBDIR}
 else ifeq ($(USE_SOLVER), MA27)
 	LIB_SOLVER = /usr/local/lib/libhsl_ma27.a
 	DEF_SOLVER = SOLVER_MA27
+	LINKHSL =
 else
 	LIB_SOLVER =
 	DEF_SOLVER = SOLVER_NONE
+	LINKHSL =
 endif
 
 ################################################################################
@@ -101,16 +107,16 @@ endif
 
 
 
-CPPFLAGS = -Wall -pedantic -Wshadow -Wfloat-equal -O3 -Wconversion -Wsign-conversion -finline-functions -fPIC -DLINUX -D${DEF_SOLVER} -D__NO_COPYRIGHT__
+CPPFLAGS = -Wall -pedantic -Wshadow -Wfloat-equal -O3 -Wconversion -Wsign-conversion -finline-functions -fPIC -DLINUX -D__USE_LONG_INTEGERS__ -D__USE_LONG_FINTS__ -D${DEF_SOLVER} -D__NO_COPYRIGHT__
 #          -g -D__DEBUG__ -D__NO_COPYRIGHT__ -D__SUPPRESSANYOUTPUT__ -D__USE_SINGLE_PRECISION__
 
 # libraries to link against when building qpOASES .so files
 LINK_LIBRARIES = ${LIB_LAPACK} ${LIB_BLAS} -lm ${LIB_SOLVER}
 LINK_LIBRARIES_WRAPPER = -lm ${LIB_SOLVER} -lstdc++
 
 # how to link against the qpOASES shared library
-QPOASES_LINK = -L${BINDIR} -Wl,-rpath=${BINDIR} -lqpOASES
-QPOASES_LINK_WRAPPER = -L${BINDIR} -Wl,-rpath=${BINDIR} -lqpOASES_wrapper
+QPOASES_LINK = -L${BINDIR} -Wl,-rpath=${BINDIR} ${LINKHSL} -lqpOASES
+QPOASES_LINK_WRAPPER = -L${BINDIR} -Wl,-rpath=${BINDIR} ${LINKHSL} -lqpOASES_wrapper
 
 # link dependencies when creating executables
 LINK_DEPENDS = ${LIB_LAPACK} ${LIB_BLAS} ${BINDIR}/libqpOASES.${LIBEXT} ${BINDIR}/libqpOASES.${DLLEXT}