diff --git a/.gitattributes b/.gitattributes index f233d00..7b013b1 100644 --- a/.gitattributes +++ b/.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 diff --git a/include/qpOASES.hpp b/include/qpOASES.hpp index 9d51cf9..5726e14 100644 --- a/include/qpOASES.hpp +++ b/include/qpOASES.hpp @@ -38,8 +38,10 @@ #include #include #include +#ifndef __MATLAB__ #include #include +#endif #include #include #include diff --git a/include/qpOASES/Types.hpp b/include/qpOASES/Types.hpp index 9654e3a..9b4f1f5 100644 --- a/include/qpOASES/Types.hpp +++ b/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__ */ /** diff --git a/interfaces/matlab/make.m b/interfaces/matlab/make.m index 1ba7768..727d87d 100644 --- a/interfaces/matlab/make.m +++ b/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 ) diff --git a/interfaces/matlab/qpOASES.cpp b/interfaces/matlab/qpOASES.cpp index 957636f..a1eebd1 100644 --- a/interfaces/matlab/qpOASES.cpp +++ b/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; diff --git a/interfaces/matlab/qpOASES_matlab_utils.cpp b/interfaces/matlab/qpOASES_matlab_utils.cpp index ee96ee2..c2c3fe9 100644 --- a/interfaces/matlab/qpOASES_matlab_utils.cpp +++ b/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 ); diff --git a/interfaces/matlab/qpOASES_matlab_utils.hpp b/interfaces/matlab/qpOASES_matlab_utils.hpp index 135cdb9..99756c0 100644 --- a/interfaces/matlab/qpOASES_matlab_utils.hpp +++ b/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; diff --git a/interfaces/matlab/qpOASES_sequence.cpp b/interfaces/matlab/qpOASES_sequence.cpp index 95b3038..8c2992b 100644 --- a/interfaces/matlab/qpOASES_sequence.cpp +++ b/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) */ diff --git a/interfaces/matlab/testQPset.m b/interfaces/matlab/testQPset.m new file mode 100644 index 0000000..40981f7 --- /dev/null +++ b/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 + diff --git a/interfaces/matlab/testSchur.m b/interfaces/matlab/testSchur.m new file mode 100644 index 0000000..8ac27c1 --- /dev/null +++ b/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); + diff --git a/make_linux.mk b/make_linux.mk index 0670f11..cbc1fbe 100644 --- a/make_linux.mk +++ b/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,7 +107,7 @@ 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 @@ -109,8 +115,8 @@ 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} diff --git a/src/LAPACKReplacement.cpp b/src/LAPACKReplacement.cpp index 61a4a98..f64cc71 100644 --- a/src/LAPACKReplacement.cpp +++ b/src/LAPACKReplacement.cpp @@ -123,7 +123,7 @@ extern "C" void dtrtrs_( const char *UPLO, const char *TRANS, const char *DIAG, double *A, const unsigned long *LDA, double *B, const unsigned long *LDB, long *INFO ) { - ; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ + INFO[0] = 0xDEADBEEF; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ } extern "C" void strtrs_( const char *UPLO, const char *TRANS, const char *DIAG, @@ -131,7 +131,7 @@ extern "C" void strtrs_( const char *UPLO, const char *TRANS, const char *DIAG, float *A, const unsigned long *LDA, float *B, const unsigned long *LDB, long *INFO ) { - ; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ + INFO[0] = 0xDEADBEEF; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ } extern "C" void dtrcon_( const char *NORM, const char *UPLO, const char *DIAG, @@ -139,7 +139,7 @@ extern "C" void dtrcon_( const char *NORM, const char *UPLO, const char *DIAG, double *RCOND, double *WORK, const unsigned long *IWORK, long *INFO ) { - ; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ + INFO[0] = 0xDEADBEEF; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ } extern "C" void strcon_( const char *NORM, const char *UPLO, const char *DIAG, @@ -147,5 +147,6 @@ extern "C" void strcon_( const char *NORM, const char *UPLO, const char *DIAG, float *RCOND, float *WORK, const unsigned long *IWORK, long *INFO ) { - ; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ + INFO[0] = 0xDEADBEEF; /* Dummy. If SQProblemSchur is to be used, system LAPACK must be used */ } + diff --git a/src/SQProblemSchur.cpp b/src/SQProblemSchur.cpp index e70aa4b..ffcbd15 100644 --- a/src/SQProblemSchur.cpp +++ b/src/SQProblemSchur.cpp @@ -2456,6 +2456,8 @@ returnValue SQProblemSchur::backsolveSchurQR( int_t dimS, const real_t* const rh if ( INFO != 0 ) { MyPrintf("TRTRS returns INFO = %d\n", INFO); + if ( INFO == 0xDEADBEEF ) + MyPrintf( "If SQProblemSchur is to be used, system LAPACK must be used instead of the qpOASES LAPACK replacement" ); return RET_QR_FACTORISATION_FAILED; } diff --git a/src/SparseSolver.cpp b/src/SparseSolver.cpp index defba5e..84ab424 100644 --- a/src/SparseSolver.cpp +++ b/src/SparseSolver.cpp @@ -781,6 +781,9 @@ returnValue Ma57SparseSolver::factorize( ) fint *keep_ma57 = new fint[lkeep_ma57]; fint *iwork_ma57 = new fint[5*dim]; + /* Set initial pivot sequence. */ + for (fint i = 0; i < lkeep_ma57; i++) keep_ma57[i] = i+1; + fint info_ma57[40]; double rinfo_ma57[20]; diff --git a/testing/cpp/Makefile b/testing/cpp/Makefile index 6f62d1c..2e6dc49 100644 --- a/testing/cpp/Makefile +++ b/testing/cpp/Makefile @@ -56,6 +56,7 @@ QPOASES_TEST_EXES = \ ${BINDIR}/test_exampleLP${EXE} \ ${BINDIR}/test_qrecipe${EXE} \ ${BINDIR}/test_qrecipeSchur${EXE} \ + ${BINDIR}/test_smallSchur${EXE} \ ${BINDIR}/test_infeasible1${EXE} \ ${BINDIR}/test_hs268${EXE} \ ${BINDIR}/test_gradientShift${EXE} \ diff --git a/testing/cpp/test_qrecipeSchur.cpp b/testing/cpp/test_qrecipeSchur.cpp index e5f9ff2..0e3630f 100644 --- a/testing/cpp/test_qrecipeSchur.cpp +++ b/testing/cpp/test_qrecipeSchur.cpp @@ -56,6 +56,10 @@ int main( ) real_t *x3 = new real_t[180]; real_t *y3 = new real_t[271]; + Options options; + options.setToDefault(); + options.enableEqualities = BT_TRUE; + /* create sparse matrices */ SymSparseMat *H = new SymSparseMat(180, 180, H_ir, H_jc, H_val); SparseMatrix *A = new SparseMatrix(91, 180, A_ir, A_jc, A_val); @@ -71,6 +75,7 @@ int main( ) /* solve with dense matrices */ nWSR = 1000; QProblem qrecipeD(180, 91); + qrecipeD.setOptions(options); tic = getCPUtime(); qrecipeD.init(Hd, g, Ad, lb, ub, lbA, ubA, nWSR, 0); toc = getCPUtime(); @@ -82,6 +87,7 @@ int main( ) /* solve with sparse matrices (nullspace factorization) */ nWSR = 1000; QProblem qrecipeS(180, 91); + qrecipeS.setOptions(options); tic = getCPUtime(); qrecipeS.init(H, g, A, lb, ub, lbA, ubA, nWSR, 0); toc = getCPUtime(); @@ -94,6 +100,7 @@ int main( ) #ifndef SOLVER_NONE nWSR = 1000; SQProblemSchur qrecipeSchur(180, 91); + qrecipeSchur.setOptions(options); tic = getCPUtime(); qrecipeSchur.init(H, g, A, lb, ub, lbA, ubA, nWSR, 0); toc = getCPUtime(); diff --git a/testing/cpp/test_smallSchur.cpp b/testing/cpp/test_smallSchur.cpp new file mode 100644 index 0000000..7dba6b9 --- /dev/null +++ b/testing/cpp/test_smallSchur.cpp @@ -0,0 +1,140 @@ +/* + * This file is part of qpOASES. + * + * qpOASES -- An Implementation of the Online Active Set Strategy. + * Copyright (C) 2007-2015 by Hans Joachim Ferreau, Andreas Potschka, + * Christian Kirches et al. All rights reserved. + * + * qpOASES is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * qpOASES is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * See the GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with qpOASES; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + * + */ + + +/** + * \file testing/cpp/test_smallSchur.cpp + * \author Dennis Janka + * \version 3.2 + * \date 2007-2015 + * + */ + + + +#include +#include + +int main( ) +{ + USING_NAMESPACE_QPOASES + + int_t i; + const int_t m = 200; + const int_t n = 2*m + 1; + + /* problem data */ + sparse_int_t H_jc[n+1], H_ir[n], A_jc[n+1], A_ir[m]; + real_t H_val[n], A_val[m], g[n], lb[n], ub[n], lbA[m], ubA[m]; + for (i = 0; i < n+1; i++) H_jc[i] = i; + for (i = 0; i < n; i++) H_ir[i] = i; + for (i = 0; i < n; i++) H_val[i] = 2.0; + for (i = 0; i < m; i++) A_jc[i] = i; + for (i = m; i < n+1; i++) A_jc[i] = m; + for (i = 0; i < m; i++) A_ir[i] = i; + for (i = 0; i < m; i++) A_val[i] = 1.0; + for (i = 0; i < n; i++) g[i] = 0.0; + for (i = 0; i < n; i++) lb[i] = -1.0; + for (i = 0; i < n; i++) ub[i] = 1.0; + for (i = 0; i < m; i++) lbA[i] = 0.5; + for (i = 0; i < m; i++) ubA[i] = 0.5; + + int_t nWSR; + real_t errP, errD, tic, toc; + real_t *xref = new real_t[n]; + real_t *yref = new real_t[n+m]; + real_t *x = new real_t[n]; + real_t *y = new real_t[n+m]; + + Options options; + options.setToDefault(); + options.printLevel = PL_TABULAR; + options.initialStatusBounds = ST_UPPER; + + /* create sparse matrices */ + SymSparseMat *H = new SymSparseMat(n, n, H_ir, H_jc, H_val); + SparseMatrix *A = new SparseMatrix(m, n, A_ir, A_jc, A_val); + + H->createDiagInfo(); + + /* solve with dense matrices */ + nWSR = 1000; + QProblem qpD(n, m); + qpD.setOptions(options); + tic = getCPUtime(); + qpD.init(H, g, A, lb, ub, lbA, ubA, nWSR, 0); + toc = getCPUtime(); + qpD.getPrimalSolution(xref); + qpD.getDualSolution(yref); + + fprintf(stdFile, "Solved problem with dense LA in %d iterations, %.3f seconds.\n", (int)nWSR, toc-tic); + + /* solve with sparse matrices (Schur complement) */ + #ifndef SOLVER_NONE + nWSR = 1000; + SQProblemSchur qp(n, m); + qp.setOptions(options); + tic = getCPUtime(); + qp.init(H, g, A, lb, ub, lbA, ubA, nWSR, 0); + toc = getCPUtime(); + qp.getPrimalSolution(x); + qp.getDualSolution(y); + + fprintf(stdFile, "Solved problem with Schur complement approach in %d iterations, %.3f seconds.\n", (int)nWSR, toc-tic); + #endif /* SOLVER_NONE */ + + /* check distance of solutions */ + errP = 0.0; + #ifndef SOLVER_NONE + for (i = 0; i < n; i++) + if (getAbs(x[i] - xref[i]) > errP) + errP = getAbs(x[i] - xref[i]); + #endif /* SOLVER_NONE */ + fprintf(stdFile, "Primal error: %9.2e\n", errP); + + errD = 0.0; + #ifndef SOLVER_NONE + for (i = 0; i < n+m; i++) + if (getAbs(y[i] - yref[i]) > errD) + errD = getAbs(y[i] - yref[i]); + #endif /* SOLVER_NONE */ + fprintf(stdFile, "Dual error: %9.2e\n", errD); + + delete H; + delete A; + + delete[] y; + delete[] x; + delete[] yref; + delete[] xref; + + QPOASES_TEST_FOR_TOL( errP,1e-12 ); + QPOASES_TEST_FOR_TOL( errD,1e-12 ); + + return TEST_PASSED; +} + + +/* + * end of file + */ diff --git a/testing/runUnitTests b/testing/runUnitTests index 3a1203e..20298c4 100755 --- a/testing/runUnitTests +++ b/testing/runUnitTests @@ -78,6 +78,7 @@ runTest $counter ../bin/test_example5; runTest $counter ../bin/test_exampleLP; runTest $counter ../bin/test_qrecipe; runTest $counter ../bin/test_qrecipeSchur; +runTest $counter ../bin/test_smallSchur; runTest $counter ../bin/test_infeasible1; runTest $counter ../bin/test_hs268;