Changeset 261


Ignore:
Timestamp:
Aug 9, 2017 11:58:08 AM (22 months ago)
Author:
dcebulla
Message:

Changes in Makefiles (presolver can be turned off). IMPORTANT: Interfaces must also be adapted!

Location:
branches/presolver
Files:
16 added
20 edited

Legend:

Unmodified
Added
Removed
  • branches/presolver/Makefile

    r260 r261  
    3838##
    3939
    40 
    4140all: qpPresolver src examples
    4241#src_aw testing
    4342
     43# We only do not compile the presolver if the user tells us!
    4444qpPresolver:
     45ifeq (${NO_PRESOLVER}, 0)
    4546        @cd $@; ${MAKE} -s
     47endif
    4648
    4749src:
  • branches/presolver/examples/Makefile

    r260 r261  
    3737
    3838IFLAGS      =  -I. \
    39                -I${IDIR} \
    40                -I${PRESOLVER_IDIR}
     39               -I${IDIR}
     40
     41ifeq (${NO_PRESOLVER}, 0)
     42IFLAGS += \
     43    -I${PRESOLVER_IDIR} \
     44    -I${SUITESPARSE_IDIR}
     45endif
    4146
    4247QPOASES_EXES = \
     
    5156        ${BINDIR}/exampleLP${EXE} \
    5257        ${BINDIR}/qrecipe${EXE} \
    53         ${BINDIR}/qrecipeSchur${EXE} \
    54         ${BINDIR}/qrecipePresolverDense${EXE} \
     58        ${BINDIR}/qrecipeSchur${EXE}
     59
     60ifeq (${NO_PRESOLVER}, 0)
     61QPOASES_EXES += \
     62    ${BINDIR}/qrecipePresolverDense${EXE} \
    5563        ${BINDIR}/qrecipePresolverSparse${EXE}
    56 
     64endif
    5765
    5866##
  • branches/presolver/include/qpOASES.hpp

    r260 r261  
    4949#include <SparseSolver.cpp>
    5050#include <SQProblemSchur.cpp>
     51
     52#ifndef NO_PRESOLVER
    5153#include <Presolver.cpp>
     54#endif
    5255
    5356#ifndef __C_WRAPPER__
     
    6467#include <qpOASES/extras/OQPinterface.hpp>
    6568#include <qpOASES/extras/SolutionAnalysis.hpp>
     69
     70#ifndef NO_PRESOLVER
    6671#include <qpOASES/Presolver.hpp>
     72#endif
    6773
    6874#endif
  • branches/presolver/include/qpOASES/Presolver.hpp

    r260 r261  
    5858#define MST_COLUMN_WISE         QPP_MST_COLUMN_WISE
    5959
    60 /** \brief Provides functionalities to presolve a given QP (primarily to reduce its size).
     60/** \brief Provides functionalities to presolve QPs (primarily to reduce their size).
    6161 *
    6262 *      The presolver must be used before solving the QP (via one of its presolve() methods)
     
    6767 *      CAVEAT: The impact of preprocessing depends on the QP data, i.e. on the matrices,
    6868 *      gradient vector and bounds. Care must be taken if a sequence of QPs shall be presolved.
    69  *      We only recommend to presolve the initial QP, then solving it, afterwards postsolving
    70  *      it (i.e. retrieving the optimal working set and primal-dual solution of the
    71  *      original QP). The working set can be used for a warm-start.
     69 *      We only recommend to presolve the initial QP, solving it, then postsolving it
     70 *      (i.e. retrieving the optimal working set and primal-dual solution of the
     71 *      original QP). Finally, the computed working set may be used as a warm-start for
     72 *      the sequence of QPs.
    7273 */
    7374class Presolver
     
    7980         *      Allocates (most of the) memory the presolver requires (e.g. for storing the QP).
    8081         *      If the number of nonzero elements of the matrices are not known, the presolver
    81          *      reallocates memory later.
     82         *      reallocates memory for the matrices later.
    8283         *
    8384         *      \param nV Number of variables of the QP.
    8485         *      \param nC Number of constraints of the QP.
    8586         *      \param nzH Number of nonzero elements of the Hessian matrix. Default value: 0.
    86          *              Note that the presolver only needs access to the lower triangular part of
    87          *              the Hessian matrix and hence the number of nonzero elements of the lower
    88          *              triangular part of the Hessian matrix must be passed.
     87         *              Note that the presolver only accesses the lower triangular part of
     88         *              the Hessian matrix and hence only the number of nonzero elements of the lower
     89         *              triangular part of the Hessian matrix must be passed to the presolver.
    8990         *      \param nzA Number of nonzero elements of the constraint matrix. Default value: 0.
    9091         *      \param presolveStackSize Initial size (= number of elements) of the presolve stack.
     
    103104
    104105
     106        /** \brief Presolves the given QP with matrices stored in CSC (Harwell-Boeing) scheme.
     107         *
     108         *      All input parameters are going to be modified; they will contain the presolved QP.
     109         *      A description of all parameters can be found in \a presolveCoordMat(), except
     110         *      for those which are specifically introduced here.
     111         *
     112         *      \param Hcp Column pointers of the Hessian matrix in CSC scheme. Will be overwritten
     113         *              with the presolved Hessian matrix. Is allowed to be a null pointer if
     114         *              the Hessian matrix is empty. \n
     115         *              CAVEAT: Only the lower triangular part of the Hessian matrix will be returned!
     116         *      \param Acp Column pointers of the constraint matrix in CSC scheme. Will be
     117         *              overwritten with the presolved constraint matrix. Is allowed to be a
     118         *              null pointer if the constraint matrix is empty.
     119         *
     120         *      \return See \a presolveCoordMat()
     121         */
     122        returnValue presolve(int_t* nV,
     123                                                 int_t* nC,
     124                                                 int_t* const Hirn,
     125                                                 int_t* const Hcp,
     126                                                 real_t* const Hx,
     127                                                 int_t* nzH,
     128                                                 real_t* const g,
     129                                                 int_t* const Airn,
     130                                                 int_t* const Acp,
     131                                                 real_t* const Ax,
     132                                                 int_t* nzA,
     133                                                 real_t* const lb,
     134                                                 real_t* const ub,
     135                                                 real_t* const lbA,
     136                                                 real_t* const ubA);
     137
     138
     139        /** \brief Presolves the given QP with matrices stored in \p SparseMatrix type.
     140         *
     141         *      All parameters, except for the given matrices \p H and \p A are going to be modified;
     142         *      they contain the presolved QP. A description of all parameters can be found in \a
     143         *      presolveCoordMat(), except for those which are specifically introduced here.
     144         *
     145         *      \param H Hessian matrix. (Internal) Data is not altered. User must provide diagonal
     146         *              information via a previous call of H->createDiagInfo(). Is allowed to be a
     147         *              null pointer if the Hessian matrix is empty.
     148         *      \param Hprs Lower triangular part (!) of the presolved Hessian matrix. Memory is
     149         *              allocated in this method and must be deallocated by the user.
     150         *              Diagonal information of the matrix is not provided; the user must do this.
     151         *              If \p H is a null pointer, then \p Hprs will also be a null pointer.
     152         *      \param A Constraint matrix. (Internal) Data is not altered. Is allowed to be a
     153         *              null pointer if the constraint matrix is empty.
     154         *      \param Aprs Contains the presolved constraint matrix. Memory is allocated by
     155         *              this method and must be deallocated by the user. If \p A
     156         *              is a null pointer, then \p Aprs will also be a null pointer.
     157         *
     158         *      \return See \a presolveCoordMat().
     159         */
     160        returnValue presolve(int_t* nV,
     161                                                 int_t* nC,
     162                                                 const SymSparseMat* const H,
     163                                                 SymSparseMat** Hprs,
     164                                                 real_t* const g,
     165                                                 const SparseMatrix* const A,
     166                                                 SparseMatrix** Aprs,
     167                                                 real_t* const lb,
     168                                                 real_t* const ub,
     169                                                 real_t* const lbA,
     170                                                 real_t* const ubA);
     171
     172
     173        /**     \brief Presolves the given QP with matrices given in dense (row-major) format.
     174         *
     175         *      All input parameters are going to be modified; they contain the presolved QP.
     176         *      A description of all parameters can be found in \a presolveCoordMat(), except
     177         *      for those which are specifically introduced here.
     178         *
     179         *      \param H Dense Hessian matrix in row-major format. Will contain the (full!)
     180         *              presolved Hessian matrix. Is allowed to be a null pointer if the original
     181         *              Hessian matrix is empty.
     182         *      \param A Dense constraint matrix in row-major format. Will contain the
     183         *              presolved constraint matrix. Is allowed to be a null pointer if the original
     184         *              constraint matrix is empty. \n
     185         *
     186         *      \return See \a presolveCoordMat().
     187         */
     188        returnValue presolve(int_t* nV,
     189                                                 int_t* nC,
     190                                                 real_t* const H,
     191                                                 real_t* const g,
     192                                                 real_t* const A,
     193                                                 real_t* const lb,
     194                                                 real_t* const ub,
     195                                                 real_t* const lbA,
     196                                                 real_t* const ubA);
     197
     198
     199        /** \brief Presolves the given QP with matrices stored in \p DenseMatrix type.
     200         *
     201         *      All parameters, except for the given matrices \p H and \p A are going to be modified;
     202         *      they contain the presolved QP. A description of all parameters can be found in \a
     203         *      presolveCoordMat(), except for those which are specifically introduced here.
     204         *
     205         *      \param H Dense Hessian matrix. (Internal) Data is not altered. Is allowed to be
     206         *              a null pointer if the Hessian matrix is empty.
     207         *      \param Hprs Will point to the presolved Hessian matrix. Both the lower and
     208         *              upper triangular part are stored. Memory is allocated in this method and
     209         *              must be deallocated by the user. If \p H is a null pointer, then
     210         *              \p Hprs will also be a null pointer.
     211         *      \param A Dense constraint matrix. (Internal) Data is not altered. Is allowed to
     212         *              be a null pointer if the constraint matrix is empty.
     213         *      \param Aprs Will point to the presolved constraint matrix. Memory is allocated in
     214         *              this method and must be deallocated by the user.
     215         *              If \p A is a null pointer, then \p Aprs will also be a null pointer.
     216         *
     217         *      \return See \a presolveCoordMat().
     218         */
     219        returnValue presolve(int_t* nV,
     220                                                 int_t* nC,
     221                                                 const SymDenseMat* const H,
     222                                                 SymDenseMat** Hprs,
     223                                                 real_t* const g,
     224                                                 const DenseMatrix* const A,
     225                                                 DenseMatrix** Aprs,
     226                                                 real_t* const lb,
     227                                                 real_t* const ub,
     228                                                 real_t* const lbA,
     229                                                 real_t* const ubA);
     230
     231
     232        /**     \brief Postsolves the QP.
     233         *
     234         *      Computes a primal-dual optimal solution of the original QP on the basis of
     235         *      the primal-dual solution of the presolved QP. Tries also to compute optimal
     236         *      working sets of the original QP on the basis of the working sets of the presolved
     237         *      QP if desired.
     238         *
     239         *      \param x Primal solution of the presolved QP. Will contain the primal solution of original QP.
     240         *      \param y Dual solution of the presolved QP. Will contain the dual solution of original QP.
     241         *      \param bounds Working set w.r.t. the bound constraints. If not a null pointer,
     242         *              then this parameter will contain the optimal working set of the original QP afterwards.
     243         *      \param constraints Working set w.r.t. the linear constraints. If not a null pointer,
     244         *              then this parameter will contain the optimal working set of the original QP afterwards.
     245         *
     246         *      \return SUCCESSFUL_RETURN \n RET_ERROR_UNDEFINED \n RET_INVALID_ARGUMENTS
     247         */
     248        returnValue postsolve(real_t* const x,
     249                                                  real_t* const y,
     250                                                  Bounds* const bounds = 0,
     251                                                  Constraints* const constraints = 0);
     252
     253
     254        /** Sets the internal options based on the given options.
     255         *
     256         *      This method identifies invalid option values and sets those values to default values.
     257         *
     258         *      \param opt Options.
     259         *
     260         *      \return SUCCESSFUL_RETURN
     261         */
     262        inline returnValue setOptions(const PresolverOptions& opt);
     263
     264
     265        /** Return the current options of the presolver. */
     266        inline PresolverOptions getOptions() const;
     267
     268
     269        /** Print the current options of the presolver. */
     270        returnValue printOptions() const;
     271
     272private:
     273        qpp_data_t* data;       /**< Struct containing all necessary data for qpPresolver. */
     274
     275
    105276        /** \brief Presolves the given QP with matrices stored in coordinate scheme.
    106277         *
    107          *      All input parameters are going to be modified; they contain the presolved QP afterwards.
    108          *      Furthermore, all other presolve() methods call this function internally.
     278         *      All input parameters are going to be modified; they contain the presolved QP
     279         *      afterwards. Furthermore, all other presolve() methods call this function internally.
    109280         *
    110281         *      \param nV Will contain the number of variables of the presolved QP.
     
    143314         *      \param ubA Upper bounds on the constraints. Will be overwritten with presolved bounds.
    144315         *      \param sortType MST_ROW_WISE: Matrices of presolved QP are sorted row-wise. \n
    145          *                                      MST_COLUMN_WISE: Matrices of presolved QP are sorted column-wise.
     316         *                                      MST_COLUMN_WISE: Matrices of presolved QP are sorted column-wise (default).
    146317         *
    147318         *      \return SUCCESSFUL_RETURN \n RET_INVALID_ARGUMENTS \n RET_QP_INFEASIBLE
     
    163334                                                                 real_t* const lbA,
    164335                                                                 real_t* const ubA,
    165                                                                  const MatrixSortType& sortType);
    166 
    167 
    168         /** \brief Presolves the given QP with matrices stored in CSC (Harwell-Boeing) scheme.
    169          *
    170          *      All input parameters are going to be modified; they will contain the presolved QP.
    171          *      A description of all parameters can be found in \a presolveCoordMat(), except
    172          *      for those which are specifically introduced here.
    173          *
    174          *      \param Hcp Column pointers of the Hessian matrix in CSC scheme. Will be overwritten
    175          *              with the presolved Hessian matrix. Is allowed to be a null pointer if
    176          *              the Hessian matrix is empty. \n
    177          *              CAVEAT: Only the lower triangular part of the Hessian matrix will be returned!
    178          *      \param Acp Column pointers of the constraint matrix in CSC scheme. Will be
    179          *              overwritten with the presolved constraint matrix. Is allowed to be a
    180          *              null pointer if the constraint matrix is empty.
    181          *
    182          *      \return See \a presolveCoordMat()
    183          */
    184         returnValue presolve(int_t* nV,
    185                                                  int_t* nC,
    186                                                  int_t* const Hirn,
    187                                                  int_t* const Hcp,
    188                                                  real_t* const Hx,
    189                                                  int_t* nzH,
    190                                                  real_t* const g,
    191                                                  int_t* const Airn,
    192                                                  int_t* const Acp,
    193                                                  real_t* const Ax,
    194                                                  int_t* nzA,
    195                                                  real_t* const lb,
    196                                                  real_t* const ub,
    197                                                  real_t* const lbA,
    198                                                  real_t* const ubA);
    199 
    200 
    201         /** \brief Presolves the given QP with matrices stored in \p SparseMatrix type.
    202          *
    203          *      All parameters, except for the given matrices \p H and \p A are going to be modified;
    204          *      they contain the presolved QP. A description of all parameters can be found in \a
    205          *      presolveCoordMat(), except for those which are specifically introduced here.
    206          *
    207          *      \param H Hessian matrix. (Internal) Data is not altered. User must provide diagonal
    208          *              information via a previous call of H->createDiagInfo(). Is allowed to be a
    209          *              null pointer if the Hessian matrix is empty.
    210          *      \param Hprs Lower triangular part (!) of the presolved Hessian matrix. Memory is
    211          *              allocated in this method and must be deallocated by the user.
    212          *              Diagonal information of the matrix is not provided; the user must do this.
    213          *              If \p H is a null pointer, then \p Hprs will also be a null pointer.
    214          *      \param A Constraint matrix. (Internal) Data is not altered. Is allowed to be a
    215          *              null pointer if the constraint matrix is empty.
    216          *      \param Aprs Contains the presolved constraint matrix. Memory is allocated by
    217          *              this method and must be deallocated by the user. If \p A
    218          *              is a null pointer, then \p Aprs will also be a null pointer.
    219          *
    220          *      \return See \a presolveCoordMat().
    221          */
    222         returnValue presolve(int_t* nV,
    223                                                  int_t* nC,
    224                                                  const SymSparseMat* const H,
    225                                                  SymSparseMat** Hprs,
    226                                                  real_t* const g,
    227                                                  const SparseMatrix* const A,
    228                                                  SparseMatrix** Aprs,
    229                                                  real_t* const lb,
    230                                                  real_t* const ub,
    231                                                  real_t* const lbA,
    232                                                  real_t* const ubA);
    233 
    234 
    235         /**     \brief Presolves the given QP with matrices given in dense (row-major) format.
    236          *
    237          *      All input parameters are going to be modified; they contain the presolved QP.
    238          *      A description of all parameters can be found in \a presolveCoordMat(), except
    239          *      for those which are specifically introduced here.
    240          *
    241          *      \param H Dense Hessian matrix in row-major format. Will contain the (full!)
    242          *              presolved Hessian matrix. Is allowed to be a null pointer if the original
    243          *              Hessian matrix is empty.
    244          *      \param A Dense constraint matrix in row-major format. Will contain the
    245          *              presolved constraint matrix. Is allowed to be a null pointer if the original
    246          *              constraint matrix is empty. \n
    247          *
    248          *      \return See \a presolveCoordMat().
    249          */
    250         returnValue presolve(int_t* nV,
    251                                                  int_t* nC,
    252                                                  real_t* const H,
    253                                                  real_t* const g,
    254                                                  real_t* const A,
    255                                                  real_t* const lb,
    256                                                  real_t* const ub,
    257                                                  real_t* const lbA,
    258                                                  real_t* const ubA);
    259 
    260 
    261         /** \brief Presolves the given QP with matrices stored in \p DenseMatrix type.
    262          *
    263          *      All parameters, except for the given matrices \p H and \p A are going to be modified;
    264          *      they contain the presolved QP. A description of all parameters can be found in \a
    265          *      presolveCoordMat(), except for those which are specifically introduced here.
    266          *
    267          *      \param H Dense Hessian matrix. (Internal) Data is not altered. Is allowed to be
    268          *              a null pointer if the Hessian matrix is empty.
    269          *      \param Hprs Will point to the presolved Hessian matrix. Both the lower and
    270          *              upper triangular part are stored. Memory is allocated in this method and
    271          *              must be deallocated by the user. If \p H is a null pointer, then
    272          *              \p Hprs will also be a null pointer.
    273          *      \param A Dense constraint matrix. (Internal) Data is not altered. Is allowed to
    274          *              be a null pointer if the constraint matrix is empty.
    275          *      \param Aprs Will point to the presolved constraint matrix. Memory is allocated in
    276          *              this method and must be deallocated by the user.
    277          *              If \p A is a null pointer, then \p Aprs will also be a null pointer.
    278          *
    279          *      \return See \a presolveCoordMat().
    280          */
    281         returnValue presolve(int_t* nV,
    282                                                  int_t* nC,
    283                                                  const SymDenseMat* const H,
    284                                                  SymDenseMat** Hprs,
    285                                                  real_t* const g,
    286                                                  const DenseMatrix* const A,
    287                                                  DenseMatrix** Aprs,
    288                                                  real_t* const lb,
    289                                                  real_t* const ub,
    290                                                  real_t* const lbA,
    291                                                  real_t* const ubA);
    292 
    293 
    294         /**     \brief Postsolves the QP.
    295          *
    296          *      Computes a primal-dual optimal solution of the original QP on the basis of
    297          *      the primal-dual solution of the presolved QP. Tries also to compute optimal
    298          *      working sets of the original QP on the basis of the working sets of the presolved
    299          *      QP if desired.
    300          *
    301          *      \param x Primal solution of the presolved QP. Will contain the primal solution of original QP.
    302          *      \param y Dual solution of the presolved QP. Will contain the dual solution of original QP.
    303          *      \param bounds Working set w.r.t. the bound constraints. If not a null pointer,
    304          *              then this parameter will contain the optimal working set of the original QP afterwards.
    305          *      \param constraints Working set w.r.t. the linear constraints. If not a null pointer,
    306          *              then this parameter will contain the optimal working set of the original QP afterwards.
    307          *
    308          *      \return SUCCESSFUL_RETURN \n RET_ERROR_UNDEFINED \n RET_INVALID_ARGUMENTS
    309          */
    310         returnValue postsolve(real_t* const x,
    311                                                   real_t* const y,
    312                                                   Bounds* const bounds = 0,
    313                                                   Constraints* const constraints = 0);
    314 
    315 
    316         /** Sets the internal options based on the given options.
    317          *
    318          *      This method identifies invalid option values and sets those values to default values.
    319          *
    320          *      \param opt Options.
    321          *
    322          *      \return SUCCESSFUL_RETURN
    323          */
    324         inline returnValue setOptions(const PresolverOptions& opt);
    325 
    326 
    327         /** Return the current options of the presolver. */
    328         inline PresolverOptions getOptions() const;
    329 
    330 
    331         /** Print the current options of the presolver. */
    332         returnValue printOptions() const;
    333 
    334 private:
    335         qpp_data_t* data;       /**< Struct containing all necessary data for qpPresolver. */
    336 
     336                                                                 const MatrixSortType& sortType = MST_COLUMN_WISE);
    337337
    338338        /**     Converts error code from qpPresolver to error code from qpOASES (if possible).
     
    342342         *      \return SUCCESSFUL_RETURN \n RET_INVALID_ARGUMENTS \n RET_QP_INFEASIBLE
    343343         *              \n RET_QP_UNBOUNDED \n RET_ERROR_UNDEFINED \n RET_OPTIONS_ADJUSTED
     344         *              \n RET_UNABLE_TO_OPEN_FILE \n RET_UNABLE_TO_READ_FILE
     345         *              \n RET_UNABLE_TO_WRITE_FILE
    344346         */
    345347        returnValue convertErrorCode(const qpp_return_value_t err) const;
  • branches/presolver/include/qpOASES/PresolverOptions.hpp

    r260 r261  
    3737#define QPOASES_PRESOLVEROPTIONS_HPP
    3838
    39 
    4039extern "C"
    4140{
     
    5453class Presolver;
    5554
    56 /** Redefinition of bound types for qpPresolver in order to match the qpOASES
    57     naming convention. */
     55/* Redefinition of bound types for qpPresolver in order to match the qpOASES
     56   naming convention. */
    5857typedef qpp_bound_mode_t        PresolverBoundType;
    5958#define PBT_TIGHTEST            QPP_BM_TIGHTEST_BOUNDS
     
    6463 *
    6564 *      This class manages all user-specified options used for preprocessing
    66  *      quadratic programs. The options can not be changed once the QP has been presolved!
     65 *      quadratic programs. The options can not be changed once the QP has been presolved.
    6766 */
    6867class PresolverOptions : private qpp_options_t
     
    138137     *
    139138     *  Note that the passed options are not checked for consistency. This will later
    140      *  happen within the Presolver class (via \a Presolver::setOptions())
     139     *  be done within the Presolver class (via \a Presolver::setOptions()).
    141140     */
    142141    /**@{*/
     142
    143143    /** Defines which bounds are passed to the user.
    144144     *
     
    160160        inline void setMaxIter(const int_t maxIter);
    161161
    162         /** Sets the amount of information written to the logfile (value between 0 (= none) and
    163      *      5 (detailed information). We refer to the qpPresolver documentation for
    164      *      further details). */
     162        /** Sets the amount of information written to the logfile.
     163         *
     164         *  Value must be between 0 (= no information) and 5 (detailed information).
     165         *  We refer to the qpPresolver documentation for further details. */
    165166        inline void setLogfileLevel(const int_t logLevel);
    166167        /**@}*/
     
    168169        /** \name Methods for enabling/disabling certain preprocessing methods.
    169170         *
    170          *  A description can be found e.g. in: Gould, N.I.M. and Toint, P.L.
     171         *  A description of the methods can be found e.g. in: Gould, N.I.M. and Toint, P.L.
    171172         *  "Preprocessing for quadratic programming." In: Mathematical Programming 100.1
    172173         *  (2004), pp. 95–132. */
    173174        /**@{*/
     175
    174176        /** Enables bound tightening on the primal variables. */
    175177        inline void enableBoundTightening();
  • branches/presolver/interfaces/matlab/Makefile

    r239 r261  
    2929##      Date:      2007-2017
    3030##
     31
    3132
    3233include ../../make.mk
  • branches/presolver/interfaces/matlab/make.m

    r255 r261  
    6565    %DEBUGFLAGS = ' -v -g CXXDEBUGFLAGS=''$CXXDEBUGFLAGS -Wall -pedantic -Wshadow'' ';
    6666
    67     IFLAGS = [ '-I. -I',QPOASESPATH,'include',' -I',QPOASESPATH,'src',' ' ];
     67    IFLAGS = [ '-I. -I',QPOASESPATH,'include',' -I',QPOASESPATH,'src',' -I', QPOASESPATH,'qpPresolver/include', ' ' ];
    6868    CPPFLAGS = [ IFLAGS, DEBUGFLAGS, '-largeArrayDims -D__cpluplus -D__MATLAB__ -D__AVOID_LA_NAMING_CONFLICTS__ -D__SINGLE_OBJECT__',' ' ];
    6969    defaultFlags = '-O -D__NO_COPYRIGHT__ '; %% -D__SUPPRESSANYOUTPUT__
  • branches/presolver/make_linux.mk

    r260 r261  
    3333# user configuration
    3434
     35# Presolver is only compiled if it is not turned off via: make NO_PRESOLVER=1
     36ifndef NO_PRESOLVER
     37NO_PRESOLVER = 0
     38endif
     39
     40# suitesparse suite (for umfpack, cholmod)
     41SUITESPARSE_IDIR   = /usr/include/suitesparse
     42SUITESPARSE_LIBDIR = /usr/lib/x86_64-linux-gnu
     43LIB_UMFPACK        = ${SUITESPARSE_LIBDIR}/libumfpack.so
     44LIB_CHOLMOD        = ${SUITESPARSE_LIBDIR}/libcholmod.so
     45LIB_SS_CONFIG      = ${SUITESPARSE_LIBDIR}/libsuitesparseconfig.so
     46
    3547# include directories for presolver, relative
    36 PRESOLVER_IDIR ${TOP}/qpPresolver/include
     48PRESOLVER_IDIR   = ${TOP}/qpPresolver/include
    3749PRESOLVER_SRCDIR = ${TOP}/qpPresolver/src
    3850PRESOLVER_LIBDIR = ${TOP}/qpPresolver/lib
     51
     52# library of presolver
     53LIB_PRESOLVER = ${PRESOLVER_LIBDIR}/libqpPresolver.a
    3954
    4055# include directories, relative
     
    4257SRCDIR = ${TOP}/src
    4358BINDIR = ${TOP}/bin
    44 
    45 # library of presolver
    46 LIB_PRESOLVER = ${PRESOLVER_LIBDIR}/libqpPresolver.a
    4759
    4860# Matlab include directory (ADAPT TO YOUR LOCAL SETTINGS!)
     
    113125#          -g -D__DEBUG__ -D__NO_COPYRIGHT__ -D__SUPPRESSANYOUTPUT__ -D__USE_SINGLE_PRECISION__
    114126
     127ifneq (${NO_PRESOLVER}, 0)
     128CPPFLAGS += -DNO_PRESOLVER
     129endif
     130
    115131# libraries to link against when building qpOASES .so files
    116 LINK_LIBRARIES = ${LIB_LAPACK} ${LIB_BLAS} ${LIB_PRESOLVER} -lm ${LIB_SOLVER}
    117 LINK_LIBRARIES_WRAPPER = -lm ${LIB_SOLVER} -lstdc++
     132ifeq (${NO_PRESOLVER}, 0)
     133LINK_LIBRARIES = ${LIB_UMFPACK} ${LIB_CHOLMOD} ${LIB_SS_CONFIG} ${LIB_PRESOLVER}
     134endif
     135LINK_LIBRARIES += ${LIB_LAPACK} ${LIB_BLAS} -lm ${LIB_SOLVER}
     136
     137# CHECK!
     138ifeq (${NO_PRESOLVER}, 0)
     139LINK_LIBRARIES_WRAPPER = ${LIB_UMFPACK} ${LIB_CHOLMOD} ${LIB_SS_CONFIG} ${LIB_PRESOLVER}
     140endif
     141LINK_LIBRARIES_WRAPPER += -lm ${LIB_SOLVER} -lstdc++
    118142
    119143# how to link against the qpOASES shared library
  • branches/presolver/qpPresolver/Makefile

    r260 r261  
    1 # Makefile for compiling the qpPresolver library and its dependencies.
    2 # Use only this Makefile or the installation scripts!!!
    3 
    4 # ==============================================================================
    5 # integer type
    6 # ==============================================================================
    7 #
    8 # Per default QPP_USE_INT64 is not defined. Call
    9 #
    10 # $ make DEFINES="-DQPP_USE_INT64"
    11 #
    12 # to compile with 64 bit integers.
    13 
    14 # ==============================================================================
    15 # single precision
    16 # ==============================================================================
    17 #
    18 # Per default QPP_USE_SINGLE_PRECISION is not defined. Call
    19 #
    20 # $ make DEFINES="-DQPP_USE_SINGLE_PRECISION"
    21 #
    22 # to use single precision floating point numbers.
     1##
     2##      This file is part of qpPresolver.
     3##
     4##      qpPresolver -- An implementation of presolving (= preprocessing) techniques
     5##  for Quadratic Programming.
     6##      Copyright (C) 2017 by Dominik Cebulla et al. All rights reserved.
     7##
     8##      qpPresolver is free software; you can redistribute it and/or
     9##      modify it under the terms of the GNU Lesser General Public
     10##      License as published by the Free Software Foundation; either
     11##      version 2.1 of the License, or (at your option) any later version.
     12##
     13##      qpPresolver is distributed in the hope that it will be useful,
     14##      but WITHOUT ANY WARRANTY; without even the implied warranty of
     15##      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
     16##      See the GNU Lesser General Public License for more details.
     17##
     18##      You should have received a copy of the GNU Lesser General Public
     19##      License along with qpPresolver; if not, write to the Free Software
     20##      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
     21##
    2322
    2423
    25 # ==============================================================================
    26 # C compiler and flags
    27 #
    28 # CFLAGS are additional compiler flags not given to external code
    29 # ==============================================================================
     24##
     25##      Filename:  Makefile
     26##      Author:    Dominik Cebulla
     27##      Version:   1.0 Beta
     28##      Date:      2017
     29##
    3030
    31 CC = gcc
    32 CFLAGS += -Wall -Wextra -std=c99 -pedantic -Wno-unused-variable -Wshadow \
    33           -Wno-unused-function -finline-functions -Wno-float-equal \
    34           -Wno-sign-conversion -Wno-conversion -fPIC
     31include make.mk
    3532
    36 export CC
    37 export CFLAGS
     33# Several compiler options can be applied. Call
     34#     make DEFINES=
     35
     36# -DQPP_USE_INT64
     37#   to use integers with at least 8 Byte.
     38
     39# -DQPP_USE_SINGLE_PRECISION
     40#   to use single precision arithmetics (e.g. float instead of double).
     41
     42# -DQPP_WRITE_LOGFILE
     43#   to use a logfile.
     44
     45# The whole statment should be enclosed in " ", e.g. make DEFINES="-DQPP_USE_INT64"
     46
    3847export DEFINES
     48
    3949
    4050.PHONY: all mmio src clean purge
    4151
    4252
    43 # ==============================================================================
    44 # compile libraries
    45 # ==============================================================================
    46 all: src mmio
    47        
     53## TARGETS:
     54
     55all: src
     56
    4857src: mmio
    49         @cd $@; $(MAKE) -s;
     58        @${CD} $@; ${MAKE} -s;
     59
    5060
    5161mmio:
    52         @cd extern/$@; $(MAKE) -s;
    53 
    54 # ==============================================================================
    55 # clean removes the object files but keeps the libraries
    56 # ==============================================================================
    57 clean:
    58         @cd extern/mmio; $(MAKE) -s clean;
    59         @cd src; $(MAKE) -s clean;
     62        @${CD} extern/$@; ${MAKE} -s;
    6063
    6164
    62 # ==============================================================================
    63 # purge removes object files and library files
    64 # ==============================================================================
     65clean:
     66        @${CD} extern/mmio; ${MAKE} -s clean;
     67        @${CD} src; ${MAKE} -s clean;
     68
     69
    6570purge:
    66         @cd extern/mmio; $(MAKE) -s purge;
    67         @cd src; $(MAKE) -s purge;
     71        @${CD} extern/mmio; ${MAKE} -s purge;
     72        @${CD} src; ${MAKE} -s purge;
  • branches/presolver/qpPresolver/extern/mmio/Makefile

    r260 r261  
    1 # Makefile for compiling the mmio code. No (static or dynamic) library will be
    2 # created, only the object file.
     1##
     2##      This file is part of qpPresolver.
     3##
     4##      qpPresolver -- An implementation of presolving (= preprocessing) techniques
     5##  for Quadratic Programming.
     6##      Copyright (C) 2017 by Dominik Cebulla et al. All rights reserved.
     7##
     8##      qpPresolver is free software; you can redistribute it and/or
     9##      modify it under the terms of the GNU Lesser General Public
     10##      License as published by the Free Software Foundation; either
     11##      version 2.1 of the License, or (at your option) any later version.
     12##
     13##      qpPresolver is distributed in the hope that it will be useful,
     14##      but WITHOUT ANY WARRANTY; without even the implied warranty of
     15##      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
     16##      See the GNU Lesser General Public License for more details.
     17##
     18##      You should have received a copy of the GNU Lesser General Public
     19##      License along with qpPresolver; if not, write to the Free Software
     20##      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
     21##
    322
    4 # ==============================================================================
    5 # header files for dependencies
    6 # ==============================================================================
    7 HDR = include/mmio.h
     23##
     24##      Filename:  extern/mmio/Makefile
     25##      Author:    Dominik Cebulla
     26##      Version:   1.0 Beta
     27##      Date:      2017
     28##
    829
    9 # ==============================================================================
    10 # object files
    11 # ==============================================================================
    12 OBJ = lib/mmio.o
     30include ../../make.mk
    1331
    14 # ==============================================================================
    15 # library files
    16 # ==============================================================================
    17 # LIB = lib/libmmio.a
    1832
    19 # ==============================================================================
    20 # compiler command for generating objects
    21 # ==============================================================================
    22 C = $(CC) $(CFLAGS) $(DEFINES) -Iinclude
     33MMIO_DEPENDS  = include/mmio.h
     34MMIO_OBJECTS  = src/mmio.o
     35MMIO_IFLAGS   = -Iinclude
     36MMIO_LIB_DIR  = lib
     37MMIO_LIB      = ${MMIO_LIB_DIR}/libmmio.a
    2338
    24 # ==============================================================================
    25 # command for creating lib directory if not existent
    26 # ==============================================================================
    27 MKDIR_P = mkdir -p
    2839
    29 # ==============================================================================
    30 # build library
    31 # ==============================================================================
    32 # $(LIB): $(OBJ)
    33 #       @$(AR) rs $(LIB) $(OBJ)
     40C = ${CC} ${CFLAGS} ${DEFINES} ${MMIO_IFLAGS}
    3441
    35 # ==============================================================================
    36 # compile all object files
    37 # ==============================================================================
    38 lib/mmio.o: src/mmio.c $(HDR)
    39         @$(MKDIR_P) -p lib
    40         @$(C) -c $< -o $@
    4142
    42 # ==============================================================================
    43 # remove object files
    44 # ==============================================================================
     43${MMIO_LIB}: ${MMIO_OBJECTS}
     44        @${MKDIR} ${MMIO_LIB_DIR}
     45        @${AR} rs ${MMIO_LIB} ${MMIO_OBJECTS}
     46
     47
     48${MMIO_OBJECTS}: src/mmio.c
     49        @${C} -c $< ${DEF_TARGET}
     50
     51
    4552clean:
    46         @$(RM) $(OBJ)
     53        @${RM} ${MMIO_OBJECTS}
    4754
    48 # ==============================================================================
    49 # remove object and library files
    50 # ==============================================================================
     55
    5156purge: clean
    52 #       @$(RM) $(LIB)
     57        @${RM} ${MMIO_LIB}
  • branches/presolver/qpPresolver/include/qpPresolver.h

    r260 r261  
    3838
    3939/* qpOASES defines: */
     40
    4041#ifdef __USE_SINGLE_PRECISION__
    4142    #define QPP_USE_SINGLE_PRECISION
     
    5051#endif
    5152
    52 
     53#include <qpPresolver/cholmod_config.h>
    5354#include <qpPresolver/constants.h>
    5455#include <qpPresolver/ecrmatrix.h>
  • branches/presolver/qpPresolver/include/qpPresolver/presolver.h

    r260 r261  
    4545#include <time.h>
    4646
     47#include <qpPresolver/cholmod_config.h>
    4748#include <qpPresolver/constants.h>
    4849#include <qpPresolver/ecrmatrix.h>
  • branches/presolver/qpPresolver/include/qpPresolver/types.h

    r260 r261  
    4949
    5050#ifdef QPP_USE_INT64
    51     //typedef int_least64_t qpp_int_t;    /**< We use integers with at least 64 bit length. */
    5251    typedef long int qpp_int_t;             /**< Integer type. */
    5352    typedef unsigned long int qpp_uint_t;   /**< Unsigned integer type. */
    5453#else
    55     //typedef int_least32_t qpp_int_t;    /**< We use integers with at least 32 bit length. */
    5654    typedef int qpp_int_t;                  /**< Integer type. */
    5755    typedef unsigned int qpp_uint_t;        /**< Unsigned integer type. */
  • branches/presolver/qpPresolver/src/Makefile

    r260 r261  
    1 # Makefile for creating the qpPresolver library.
     1##
     2##      This file is part of qpPresolver.
     3##
     4##      qpPresolver -- An implementation of presolving (= preprocessing) techniques
     5##  for Quadratic Programming.
     6##      Copyright (C) 2017 by Dominik Cebulla et al. All rights reserved.
     7##
     8##      qpPresolver is free software; you can redistribute it and/or
     9##      modify it under the terms of the GNU Lesser General Public
     10##      License as published by the Free Software Foundation; either
     11##      version 2.1 of the License, or (at your option) any later version.
     12##
     13##      qpPresolver is distributed in the hope that it will be useful,
     14##      but WITHOUT ANY WARRANTY; without even the implied warranty of
     15##      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
     16##      See the GNU Lesser General Public License for more details.
     17##
     18##      You should have received a copy of the GNU Lesser General Public
     19##      License along with qpPresolver; if not, write to the Free Software
     20##      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
     21##
    222
    3 EXTDIR = ../extern
     23##
     24##      Filename:  src/Makefile
     25##      Author:    Dominik Cebulla
     26##      Version:   1.0 Beta
     27##      Date:      2017
     28##
    429
    5 LIBDIR = ../lib
    6 QPPDIR = ../include/qpPresolver
     30include ../make.mk
    731
    8 MKDIR_P = mkdir -p
     32QPPRESOLVER_EXTERN_DIR   = ../extern
     33QPPRESOLVER_LIB_DIR      = ../lib
     34QPPRESOLVER_INCLUDE_DIR  = ../include
    935
    10 # ==============================================================================
    11 # header files for dependencies
    12 # ==============================================================================
    13 HDR = ../include/qpPresolver.h \
    14       $(QPPDIR)/constants.h \
    15       $(QPPDIR)/ecrmatrix.h \
    16       $(QPPDIR)/errorcodes.h \
    17       $(QPPDIR)/minheap.h \
    18       $(QPPDIR)/options.h \
    19       $(QPPDIR)/presolver.h \
    20       $(QPPDIR)/stack.h \
    21       $(QPPDIR)/types.h \
    22       $(QPPDIR)/utility.h \
    23       $(EXTDIR)/mmio/include/mmio.h
    24      
    25 INCLUDES = $(EXTDIR)/mmio/include
     36SUITESPARSE_INCLUDE_DIR  = /usr/include/suitesparse
    2637
    27 # ==============================================================================
    28 # object files
    29 # ==============================================================================
    30 OBJ = $(LIBDIR)/ecrmatrix.o \
    31           $(LIBDIR)/minheap.o \
    32           $(LIBDIR)/options.o \
    33           $(LIBDIR)/presolver.o \
    34           $(LIBDIR)/stack.o \
    35           $(LIBDIR)/utility.o
     38QPPRESOLVER_IFLAGS       = -I../include \
     39                           -I${QPPRESOLVER_EXTERN_DIR}/mmio/include \
     40                           -I${SUITESPARSE_INCLUDE_DIR}
    3641
    37 EXTERN_OBJ = $(EXTDIR)/mmio/lib/mmio.o
    3842
    39 # ==============================================================================
    40 # library files
    41 # ==============================================================================
    42 LIB = $(LIBDIR)/libqpPresolver.a
    43 SLIB = $(LIBDIR)/libqpPresolver.so
     43QPPRESOLVER_DEPENDS = ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver.h \
     44                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/cholmod_config.h \
     45                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/constants.h \
     46                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/ecrmatrix.h \
     47                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/errorcodes.h \
     48                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/minheap.h \
     49                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/options.h \
     50                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/presolver.h \
     51                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/stack.h \
     52                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/types.h \
     53                      ${QPPRESOLVER_INCLUDE_DIR}/qpPresolver/utility.h \
     54                      ${QPPRESOLVER_EXTERN_DIR}/mmio/include/mmio.h
    4455
    45 # ==============================================================================
     56
     57QPPRESOLVER_OBJECTS = ecrmatrix.${OBJEXT} \
     58                          minheap.${OBJEXT} \
     59                          options.${OBJEXT} \
     60                          presolver.${OBJEXT} \
     61                          stack.${OBJEXT} \
     62                          utility.${OBJEXT}
     63
     64QPPRESOLVER_EXTERN_OBJECTS = ${QPPRESOLVER_EXTERN_DIR}/mmio/src/mmio.${OBJEXT}
     65
     66
     67QPPRESOLVER_LIB         = ${QPPRESOLVER_LIB_DIR}/libqpPresolver.${LIBEXT}
     68QPPRESOLVER_SHARED_LIB  = ${QPPRESOLVER_LIB_DIR}/libqpPresolver.${DLLEXT}
     69
    4670# compiler command for generating objects
    47 # ==============================================================================
    48 C = $(CC) $(CFLAGS) $(DEFINES) -I../include -I$(INCLUDES)
     71C = ${CC} ${CFLAGS} ${DEFINES} ${QPPRESOLVER_IFLAGS}
    4972
    50 all: $(LIB) $(SLIB)
     73
     74all: create_lib_dir ${QPPRESOLVER_LIB} ${QPPRESOLVER_SHARED_LIB}
    5175   
    52 # ==============================================================================
    53 # compile all object files
    54 # ==============================================================================
    55 $(LIBDIR)/%.o: %.c
    56         @$(MKDIR_P) $(LIBDIR)
    57         @$(C) -c $< -o $@
     76create_lib_dir:
     77        ${MKDIR} ${QPPRESOLVER_LIB_DIR}
    5878
    59 # ==============================================================================
    60 # build static library
    61 # ==============================================================================
    62 $(LIB): $(OBJ) $(EXTERN_OBJ)
    63         @$(AR) rs $@ $^
    6479
    65 # ==============================================================================
    66 # build shared library
    67 # ==============================================================================
    68 $(SLIB): $(OBJ) $(EXTERN_OBJ)
    69         @$(C) -o $@ -shared $^
     80%.${OBJEXT}: %.c
     81        @${C} -c $< ${DEF_TARGET}
    7082
    71 # ==============================================================================
    72 # remove object files
    73 # ==============================================================================
     83
     84${QPPRESOLVER_LIB}: ${QPPRESOLVER_OBJECTS} ${QPPRESOLVER_EXTERN_OBJECTS}
     85        @${ECHO} "Creating static qpPresolver lib."
     86        @${AR} rs $@ $^
     87
     88
     89${QPPRESOLVER_SHARED_LIB}: ${QPPRESOLVER_OBJECTS} ${QPPRESOLVER_EXTERN_OBJECTS}
     90        @${ECHO} "Creating shared qpPresolver lib."
     91        @${C} ${DEF_TARGET} ${SHARED} $^ ${LINK_LIBRARIES}
     92
     93
    7494clean:
    75 #       @echo "Cleaning up (lib/*.o)"
    76         @$(RM) $(OBJ)
     95        @${ECHO} "Removing all qpPresolver/src/*.o."
     96        @${RM} ${QPPRESOLVER_OBJECTS}
    7797
    7898# ==============================================================================
     
    80100# ==============================================================================
    81101purge: clean
    82 #       @echo "Removing (shared and static) libraries"
    83         @$(RM) $(LIB) $(SLIB)
     102        @${ECHO} "Removing (shared and static) qpPresolver libraries."
     103        @${RM} ${QPPRESOLVER_LIB} ${QPPRESOLVER_SHARED_LIB}
  • branches/presolver/qpPresolver/src/presolver.c

    r260 r261  
    3737#include <qpPresolver/presolver.h>
    3838
    39 
    4039/*  ======================================================================================
    4140    Forward declaration of static (= private) functions
     
    171170                                       qpp_int_t *ub_length);
    172171
     172/** \name KKT system (linear dependency and curvature test). */
     173/**@{*/
     174
     175/** \brief Solves a KKT system given by matrices A and H. */
     176static cholmod_dense *solveKKTSystem(const qpp_ecrmatrix_t *A,
     177                                     const qpp_ecrmatrix_t *H,
     178                                     cholmod_dense *rhs,
     179                                     const qpp_int_t *av,
     180                                     const qpp_int_t num_av,
     181                                     const qpp_int_t *wi,
     182                                     const qpp_int_t num_actv_constr,
     183                                     const qpp_int_t *wj,
     184                                     const qpp_int_t num_actv_bnds,
     185                                     qpp_int_t *storage,
     186                                     qpp_int_t *vec_indices,
     187                                     qpp_real_t *vec_values);
     188
     189/** \brief Performs a curvature test */
     190static qpp_return_value_t curvatureTest(qpp_real_t *result,
     191                                        const qpp_int_t index,
     192                                        const qpp_ecrmatrix_t *A,
     193                                        const qpp_ecrmatrix_t *H,
     194                                        const qpp_int_t *av,
     195                                        const qpp_int_t num_av,
     196                                        const qpp_int_t *wi,
     197                                        const qpp_int_t num_actv_constr,
     198                                        const qpp_int_t *wj,
     199                                        const qpp_int_t num_actv_bnds,
     200                                        qpp_int_t *storage,
     201                                        qpp_int_t *vec_indices,
     202                                        qpp_real_t *vec_values);
     203
     204/** \brief Checks, if the vector given by rhs_* is linearly dependent on the constraint
     205 *      matrix (w.r.t. the working set).
     206 */
     207static qpp_return_value_t linearDependencyTest(qpp_real_t *result,
     208                                               const qpp_ecrmatrix_t *A,
     209                                               const qpp_ecrmatrix_t *H,
     210                                               const qpp_int_t *rhs_indices,
     211                                               const qpp_real_t *rhs_values,
     212                                               const qpp_int_t rhs_length,
     213                                               const qpp_int_t *av,
     214                                               const qpp_int_t num_av,
     215                                               const qpp_int_t *wi,
     216                                               const qpp_int_t num_actv_constr,
     217                                               const qpp_int_t *wj,
     218                                               const qpp_int_t num_actv_bnds,
     219                                               qpp_int_t *storage,
     220                                               qpp_int_t *vec_indices,
     221                                               qpp_real_t *vec_values);
     222
     223/** \brief Pushes a multiplier to 0 such that the working set remains linearly independent. */
     224static qpp_return_value_t ldPush(const qpp_int_t index,
     225                                 const qpp_int_t sign,
     226                                 const qpp_int_t ws,
     227                                 qpp_int_t m,
     228                                 qpp_int_t n,
     229                                 qpp_int_t *wi,
     230                                 qpp_int_t *num_wi,
     231                                 qpp_int_t *wj,
     232                                 qpp_int_t *num_wj,
     233                                 qpp_real_t *y,
     234                                 qpp_real_t *z,
     235                                 qpp_real_t *scalars,
     236                                 qpp_real_t tol);
     237
     238/**@}*/
    173239
    174240/*  ======================================================================================
     
    223289    memsize += 2*m * (qpp_int_t) sizeof(qpp_real_t);    // al, au
    224290    memsize += 3*n * (qpp_int_t) sizeof(qpp_real_t);    // zl, zu, H_diag
    225     memsize += (m+n) * (qpp_int_t) sizeof(qpp_real_t);  // mem_mpn_r
    226     memsize += (m+n) * (qpp_int_t) sizeof(qpp_int_t);   // mem_mpn_i
     291    memsize += (m+2*n) * (qpp_int_t) sizeof(qpp_real_t);  // mem_mpn_r
     292    memsize += (m+2*n) * (qpp_int_t) sizeof(qpp_int_t);   // mem_mpn_i
    227293
    228294    iptr = (qpp_int_t*) malloc((size_t) memsize);
     
    268334    data->A_nnz_rows = iptr; iptr += m;
    269335    data->H_nnz_columns = iptr; iptr += n;
    270     data->mem_mpn_i = iptr; iptr += m+n;
     336    data->mem_mpn_i = iptr; iptr += m+2*n;
    271337
    272338    rptr = (void*) iptr;
     
    283349    data->zu = rptr; rptr += n;
    284350    data->H_diag = rptr; rptr += n;
    285     data->mem_mpn_r = rptr; rptr += m+n;
     351    data->mem_mpn_r = rptr; rptr += m+2*n;
    286352
    287353    data->xl = NULL;
     
    700766
    701767
    702         //err = checkForFixVar(data);
     768        err = checkForFixVar(data);
    703769        if (err != QPP_OK)
    704770        {
     
    734800        }
    735801
     802        if (opt->enable_sparsification_method)
     803        {
     804            err = sparsification(data);
     805            if (err != QPP_OK)
     806            {
     807                break;
     808            }
     809        }
     810
    736811        if (opt->enable_empty_columns_method)
    737812        {
     
    776851        }
    777852
    778         if (opt->enable_sparsification_method)
    779         {
    780             err = sparsification(data);
    781             if (err != QPP_OK)
    782             {
    783                 break;
    784             }
    785         }
     853
    786854
    787855        /* In very rare occurrences: slow. */
     
    12401308                                const qpp_int_t use_copies)
    12411309{
    1242     qpp_int_t m, n, m_ps, n_ps, i, ii, j, jj, k, kk, it, sign, length1,
    1243         length2, vec_length_1, vec_length_2, stack_ptr;
     1310    qpp_int_t m, n, m_ps, n_ps, i, j, jj, k, kk, it, sign, length1,
     1311        length2, vec_length_1, vec_length_2, stack_ptr, num_actv_bnds, num_actv_constr;
    12441312    qpp_int_t *av, *ac, *indices1, *indices2, *vec_indices_1, *vec_indices_2, *iptr;
    1245     qpp_real_t a_ij, value, alpha, new_y, f, eq_tol;
     1313    qpp_real_t a_ij, value, alpha, new_y, f, eq_tol, feas_tol;
    12461314    qpp_real_t *xl, *xu, *R, *C, *values1, *H_diag, *g, *vec_values_1, *vec_values_2, *rptr;
    1247 
     1315    qpp_return_value_t ret_val;
    12481316    #ifdef QPP_WRITE_LOGFILE
    12491317    qpp_int_t log_level;
     
    13151383    m_ps = data->m_ps;
    13161384    n_ps = data->n_ps;
     1385    num_actv_bnds = 0;
     1386    num_actv_constr = 0;
    13171387
    13181388    if ( ((x == NULL) || (z == NULL)) && (n_ps > 0) )
     
    13981468    f = data->f_ps;
    13991469    eq_tol = opt->eq_tol;
     1470    feas_tol = opt->feas_tol;
    14001471
    14011472    /*A_nnz_rows = data->A_nnz_rows;
     
    14171488    vec_values_2  = &vec_values_1[data->mem_length];
    14181489
     1490
    14191491    /*  The elements x[0],...,x[n_ps-1] are the solution of the presolved problem.
    14201492        -> Extract solution! */
     
    14371509    if (wj != NULL)
    14381510    {
     1511        num_actv_bnds = 0;
    14391512        it = n_ps - 1;
    14401513        for (j = n-1; j >= 0; --j)
     
    14441517                wj[j] = wj[it];
    14451518                --it;
     1519                if (wj[j] != QPP_AT_INACTIVE)
     1520                {
     1521                    ++num_actv_bnds;
     1522                }
    14461523            }
    14471524            else
     
    14511528        }
    14521529    }
    1453 
    14541530
    14551531    it = m_ps - 1;
     
    14691545    if (wi != NULL)
    14701546    {
     1547        num_actv_constr = 0;
    14711548        it = m_ps - 1;
    14721549        for (i = m-1; i >= 0; --i)
     
    14761553                wi[i] = wi[it];
    14771554                --it;
     1555                if (wi[i] != QPP_AT_INACTIVE)
     1556                {
     1557                    ++num_actv_constr;
     1558                }
    14781559            }
    14791560            else
     
    14981579            y[i] = 0.0;
    14991580            ac[i] = 1;
     1581            ++m_ps;
    15001582
    15011583            /* Determine optimal working set */
     
    15261608            if ( (wj != NULL) && (wi != NULL) )
    15271609            {
     1610                assert(wj[j] != QPP_AT_EQ_CONSTRAINT);
     1611
    15281612                if ( ((tbt == QPP_TBT_BOTH) && (wj[j] != QPP_AT_INACTIVE))  ||
    15291613                     ((tbt == QPP_TBT_LOWER) && (wj[j] == QPP_AT_LOWER_BOUND)) ||
     
    15311615                {
    15321616                    y[i] = z[j] / a_ij;
    1533 
    1534                     /*if (y[i] > 0.0)
    1535                     {
    1536                         wi[i] = QPP_AT_LOWER_BOUND;
    1537                     }
    1538                     else if (y[i] <  0.0)
    1539                     {
    1540                         wi[i] = QPP_AT_UPPER_BOUND;
     1617                    z[j] = 0.0;
     1618
     1619                    if (a_ij > 0)
     1620                    {
     1621                        wi[i] = wj[j];
    15411622                    }
    15421623                    else
    1543                     {*/
    1544                         /* z_j was zero, but x_j was active -> linear constraint is active! */
    1545                         if (a_ij > 0)
    1546                         {
    1547                             wi[i] = wj[j];
    1548                         }
    1549                         else
    1550                         {
    1551                             wi[i] = -wj[j];
    1552                         }
    1553                     //}
    1554 
    1555                     z[j] = 0.0;
     1624                    {
     1625                        wi[i] = -wj[j];
     1626                    }
     1627                    ++num_actv_constr;
     1628
    15561629                    wj[j] = QPP_AT_INACTIVE;
     1630                    --num_actv_bnds;
    15571631                }
    15581632                else
     
    15641638            else    /* No working sets. We only consider the sign of the multiplier z_j. */
    15651639            {
    1566                 if ( ((tbt == QPP_TBT_BOTH) && (fabs(z[j]) > eq_tol)) ||
    1567                      ((tbt == QPP_TBT_LOWER) && (z[j] > eq_tol)) ||
    1568                      ((tbt == QPP_TBT_UPPER) && (z[j] < -eq_tol)) )
     1640                if ( ((tbt == QPP_TBT_BOTH) && (fabs(z[j]) > feas_tol)) ||
     1641                     ((tbt == QPP_TBT_LOWER) && (z[j] > feas_tol)) ||
     1642                     ((tbt == QPP_TBT_UPPER) && (z[j] < -feas_tol)) )
    15691643                {
    15701644                    y[i] = z[j] / a_ij;
     
    15841658
    15851659            ac[i] = 1;
     1660            ++m_ps;
    15861661
    15871662            /*free(sr);*/
     
    15941669
    15951670            getColH(H, j, av, vec_indices_1, vec_values_1, &vec_length_1);
    1596             for (ii = 0; ii < vec_length_1; ++ii)
    1597             {
    1598                 i = vec_indices_1[ii];
    1599                 /*++H_nnz_columns[i];*/
    1600                 g[i] -= vec_values_1[ii] * value;
    1601             }
     1671            for (kk = 0; kk < vec_length_1; ++kk)
     1672            {
     1673                k = vec_indices_1[kk];
     1674                /*++H_nnz_columns[k];*/
     1675                g[k] -= vec_values_1[kk] * value;
     1676            }
     1677            f -= (g[j] * value + 0.5 * value * value * H_diag[j]);
     1678
    16021679            av[j] = 1;
    1603             f -= (g[j] * value + 0.5 * value * value * H_diag[j]);
     1680            ++n_ps;
    16041681
    16051682            getColA(A, j, ac, vec_indices_2, vec_values_2, &vec_length_2);
     
    16261703                else
    16271704                {
    1628                     if (z[j] >= 0.0)
     1705                    if (z[j] >= -1.0e3 * eq_tol)
    16291706                    {
    16301707                        wj[j] = QPP_AT_LOWER_BOUND;
     
    16351712                    }
    16361713                }
     1714                ++num_actv_bnds;
    16371715            }
    16381716
     
    16401718            break;
    16411719
    1642         case QPP_SE_DUPLICATE_COLUMN:   /*  In some cases it is impossible to derive
    1643                                             an optimal working set! (QBRANDY)*/
     1720        case QPP_SE_DUPLICATE_COLUMN:   /*  With curvature test, works on test problems. */
    16441721            dc = se->elements.dc;
    16451722            j = dc->j;
     
    16711748
    16721749            av[j] = 1;
     1750            ++n_ps;
    16731751            value = x[k];
    16741752
     
    17011779                }
    17021780            }
     1781
    17031782            x[j] = value/alpha - x[k]/alpha;
    17041783            z[j] = alpha * z[k];
     
    17071786            if (wj != NULL)
    17081787            {
    1709                 //printf("wj[k] = %" QPP_PRID "\n", wj[k]);
    1710                 if (alpha > 0.0)
    1711                 {
    1712                     wj[j] = wj[k];
     1788                /*  Note that the bounds on x_j and x_k are correct! Variable j has been
     1789                    removed, hence its bounds were not changed afterwards. And the bounds
     1790                    on x_k have been restored in this (postprocessing) procedure. */
     1791
     1792                if (wj[k] == QPP_AT_INACTIVE)
     1793                {
     1794                    /*  Here we have to check whether we must add a bound constraint to the
     1795                        working set in order to guarantee positive definiteness of the
     1796                        reduced Hessian. */
     1797                    if (qppIsEqual(x[k], xl[k], eq_tol))
     1798                    {
     1799                        wj[k] = QPP_AT_LOWER_BOUND;
     1800                        ++num_actv_bnds;
     1801                    }
     1802                    else if (qppIsEqual(x[k], xu[k], eq_tol))
     1803                    {
     1804                        wj[k] = QPP_AT_UPPER_BOUND;
     1805                        ++num_actv_bnds;
     1806                    }
     1807
     1808                    if (wj[k] == QPP_AT_INACTIVE)
     1809                    {
     1810                        /*  x_k is not fixed at one of its bounds. Now we try for x_j. */
     1811                        if (qppIsEqual(x[j], xl[j], eq_tol))
     1812                        {
     1813                            wj[j] = QPP_AT_LOWER_BOUND;
     1814                            ++num_actv_bnds;
     1815                        }
     1816                        else if (qppIsEqual(x[j], xu[j], eq_tol))
     1817                        {
     1818                            wj[j] = QPP_AT_UPPER_BOUND;
     1819                            ++num_actv_bnds;
     1820                        }
     1821
     1822                        if (wj[j] != QPP_AT_INACTIVE)
     1823                        {
     1824                            /*  x_j is fixed at a bound -> perform curvature test
     1825                                (cf. Schork, Master's Thesis, Lemma 3.1). */
     1826                            values1 = data->mem_mpn_r;
     1827                            ret_val = curvatureTest(values1, m+j, A, H, av, n_ps, wi, num_actv_constr,
     1828                                                    wj, num_actv_bnds, data->mem_mpn_i,
     1829                                                    data->mem_for_bounds_int, data->mem_for_bounds_real);
     1830                            if (ret_val != QPP_OK)
     1831                            {
     1832                                return ret_val;
     1833                            }
     1834
     1835                            value = values1[n+m+j];
     1836
     1837                            /* Bound constraint j must stay in working set iff value <= 0. */
     1838                            if (value > feas_tol)
     1839                            {
     1840                                wj[j] = QPP_AT_INACTIVE;
     1841                                --num_actv_bnds;
     1842                            }
     1843                        }
     1844                    }
     1845                    else
     1846                    {
     1847                        /*  x_k is fixed at a bound -> perform curvature test (cf. Schork,
     1848                            Master's Thesis, Lemma 3.1). */
     1849                        values1 = data->mem_mpn_r;
     1850                        ret_val = curvatureTest(values1, m+k, A, H, av, n_ps, wi, num_actv_constr,
     1851                                                wj, num_actv_bnds, data->mem_mpn_i,
     1852                                                data->mem_for_bounds_int, data->mem_for_bounds_real);
     1853                        if (ret_val != QPP_OK)
     1854                        {
     1855                            return ret_val;
     1856                        }
     1857
     1858                        value = values1[n+m+k];
     1859
     1860                        /* Bound constraint k must stay in working set iff value <= 0. */
     1861                        if (value > feas_tol)
     1862                        {
     1863                            wj[k] = QPP_AT_INACTIVE;
     1864                            --num_actv_bnds;
     1865                        }
     1866                    }
    17131867                }
    17141868                else
    17151869                {
    1716                     wj[j] = -wj[k];
    1717                 }
    1718             }
    1719 
    1720             //printf("k = %" QPP_PRID ", j = %" QPP_PRID "\n", k, j);
     1870                    /*  Bound constr. k is in working set -> also constraint j must be
     1871                        added to the working set (both x_k and x_j are active!) */
     1872                    if (alpha > 0.0)
     1873                    {
     1874                        wj[j] = wj[k];
     1875                    }
     1876                    else
     1877                    {
     1878                        wj[j] = -wj[k];
     1879                    }
     1880                    ++num_actv_bnds;
     1881                }
     1882            }
    17211883
    17221884            /*free(dc);*/
     
    17281890
    17291891            ac[i] = 1;
     1892            ++m_ps;
    17301893            y[i] = 0.0;
    17311894
     
    18171980                for (kk = 0; kk < length1; ++kk)
    18181981                {
    1819                     int sup = 0;
    18201982                    k = indices1[kk];
    18211983
     
    18321994                        {
    18331995                            j = vec_indices_1[jj];
    1834                             if (wj[j] != QPP_AT_INACTIVE)
    1835                                 ++sup;
    1836 
     1996
     1997                            if (j == k) continue;
     1998
     1999                            /* Determine, where the bound is active. */
    18372000                            if (alpha > 0.0)
    18382001                            {
    18392002                                if (vec_values_1[jj] > 0.0)
    18402003                                {
    1841                                     wj[j] = QPP_AT_UPPER_BOUND;
     2004                                    it = QPP_AT_UPPER_BOUND;
    18422005                                }
    18432006                                else
    18442007                                {
    1845                                     wj[j] = QPP_AT_LOWER_BOUND;
     2008                                    it = QPP_AT_LOWER_BOUND;
    18462009                                }
    18472010                            }
     
    18502013                                if (vec_values_1[jj] > 0.0)
    18512014                                {
    1852                                     wj[j] = QPP_AT_LOWER_BOUND;
     2015                                    it = QPP_AT_LOWER_BOUND;
    18532016                                }
    18542017                                else
    18552018                                {
    1856                                     wj[j] = QPP_AT_UPPER_BOUND;
     2019                                    it = QPP_AT_UPPER_BOUND;
    18572020                                }
    18582021                            }
     2022
     2023                            if (wj[j] == QPP_AT_INACTIVE)
     2024                            {
     2025                                /*  We know that bound j is active. It is not part of the
     2026                                    working set because of linear dependency. */
     2027                                vec_values_2[0] = 1.0;
     2028                                vec_indices_2[0] = j;
     2029                                vec_length_2 = 1;
     2030
     2031                                values1 = data->mem_mpn_r;
     2032                                linearDependencyTest(values1, A, H, vec_indices_2, vec_values_2,
     2033                                                     vec_length_2, av, n_ps, wi, num_actv_constr,
     2034                                                     wj, num_actv_bnds, data->mem_mpn_i,
     2035                                                     data->mem_for_bounds_int, data->mem_for_bounds_real);
     2036
     2037                                /*  If e_k is part of the linear combination of e_j, then we can
     2038                                    remove e_k since it will be replaced by e_j. */
     2039                                if (fabs(values1[n+m+k]) > 1.0e3 * eq_tol)
     2040                                {
     2041                                    wj[j] = it;
     2042                                    wj[k] = QPP_AT_INACTIVE;
     2043                                }
     2044                                /*  Otherwise e_k is not part of the linear combination and another
     2045                                    constraint must be replaced for e_j. */
     2046                                else
     2047                                {
     2048                                    ldPush(j, -1, it, m, n, wi, &num_actv_constr,
     2049                                           wj, &num_actv_bnds, y, z, values1+n, eq_tol);
     2050                                }
     2051                            }
     2052                            else
     2053                            {
     2054                                wj[j] = it;
     2055                            }
    18592056                        }
    1860                         wj[k] = QPP_AT_INACTIVE;
    1861 
    1862                         if (wi[i] != QPP_AT_INACTIVE)
    1863                             ++sup;
    1864 
    1865                         if (alpha > 0.0)
     2057
     2058                        if (y[i] >= -1.0e3 * eq_tol)
    18662059                        {
    1867                             wi[i] = QPP_AT_LOWER_BOUND;
     2060                            it = QPP_AT_LOWER_BOUND;
    18682061                        }
    18692062                        else
    18702063                        {
    1871                             wi[i] = QPP_AT_UPPER_BOUND;
     2064                            it = QPP_AT_UPPER_BOUND;
    18722065                        }
    18732066
    1874                         /* If constraint i is equality constraint -> no sign restr. */
    1875                         if (y[i] > opt->feas_tol)
     2067                        if (wi[i] == QPP_AT_INACTIVE)
    18762068                        {
    1877                             wi[i] = QPP_AT_LOWER_BOUND;
     2069                            values1 = data->mem_mpn_r;
     2070                            linearDependencyTest(values1, A, H, vec_indices_1, vec_values_1,
     2071                                                 vec_length_1, av, n_ps, wi, num_actv_constr,
     2072                                                 wj, num_actv_bnds, data->mem_mpn_i,
     2073                                                 data->mem_for_bounds_int, data->mem_for_bounds_real);
     2074
     2075                            if (fabs(values1[n+m+k]) > 1.0e3 * eq_tol)
     2076                            {
     2077                                wi[i] = it;
     2078                                ++num_actv_constr;
     2079
     2080                                wj[k] = QPP_AT_INACTIVE;
     2081                                --num_actv_bnds;
     2082                            }
     2083                            else
     2084                            {
     2085                                ldPush(i, 1, it, m, n, wi, &num_actv_constr, wj,
     2086                                       &num_actv_bnds, y, z, values1+n, eq_tol);
     2087                            }
    18782088                        }
    1879 
    1880                         if (y[i] < -opt->feas_tol)
     2089                        else
    18812090                        {
    1882                             wi[i] = QPP_AT_UPPER_BOUND;
     2091                            wi[i] = it;
    18832092                        }
    1884                         printf("1 Number active components / total: %d / %" QPP_PRID "\n", sup, vec_length_1+1);
    18852093                    }
    18862094                }
     
    18892097                for (kk = 0; kk < length2; ++kk)
    18902098                {
    1891                     int sup = 0;
    1892                     int z0 = 0;
    18932099                    k = indices2[kk];
    18942100
     
    19012107
    19022108                        alpha = vec_values_1[it];       /* alpha = a_ik */
     2109
    19032110                        for (jj = 0; jj < vec_length_1; ++jj)
    19042111                        {
    19052112                            j = vec_indices_1[jj];
    19062113
    1907                             if (wj[j] == QPP_AT_INACTIVE)
    1908                                 printf("%" QPP_PRID "  ", j);
    1909 
    1910                             if (wj[j] != QPP_AT_INACTIVE)
    1911                                 ++sup;
    1912 
    1913                             if (fabs(z[j]) <= opt->feas_tol)
    1914                                 ++z0;
    1915 
    1916                             /*if (wj[j] == QPP_AT_INACTIVE)
    1917                             {
    1918                                 getColA(A, j, ac, vec_indices_2, vec_values_2, &vec_length_2);
    1919                                 int found = 0;
    1920                                 for (ii = 0; ii < vec_length_2; ++ii)
    1921                                 {
    1922                                     it = A->irowst[vec_indices_2[ii]];
    1923                                     while (it >= 0)
    1924                                     {
    1925                                         if (av[A->jcn[it]] && (fabs(z[A->jcn[it]]) <= eq_tol) && (wj[A->jcn[it]] != QPP_AT_INACTIVE))
    1926                                         {
    1927                                             wj[A->jcn[it]] = QPP_AT_INACTIVE;
    1928                                             found = 1;
    1929                                             break;
    1930                                         }
    1931                                         it = A->linkrw[it];
    1932                                     }
    1933                                     if (found) break;
    1934                                 }
    1935                             }*/
     2114                            if (j == k) continue;
    19362115
    19372116                            if (alpha > 0.0)
     
    19392118                                if (vec_values_1[jj] > 0.0)
    19402119                                {
    1941                                     wj[j] = QPP_AT_LOWER_BOUND;
     2120                                    it = QPP_AT_LOWER_BOUND;
    19422121                                }
    19432122                                else
    19442123                                {
    1945                                     wj[j] = QPP_AT_UPPER_BOUND;
     2124                                    it = QPP_AT_UPPER_BOUND;
    19462125                                }
    19472126                            }
     
    19502129                                if (vec_values_1[jj] > 0.0)
    19512130                                {
    1952                                     wj[j] = QPP_AT_UPPER_BOUND;
     2131                                    it = QPP_AT_UPPER_BOUND;
    19532132                                }
    19542133                                else
    19552134                                {
    1956                                     wj[j] = QPP_AT_LOWER_BOUND;
     2135                                    it = QPP_AT_LOWER_BOUND;
    19572136                                }
    19582137                            }
     2138
     2139                            if (wj[j] == QPP_AT_INACTIVE)
     2140                            {
     2141                                vec_values_2[0] = 1.0;
     2142                                vec_indices_2[0] = j;
     2143                                vec_length_2 = 1;
     2144
     2145                                values1 = data->mem_mpn_r;
     2146                                linearDependencyTest(values1, A, H, vec_indices_2, vec_values_2,
     2147                                                     vec_length_2, av, n_ps, wi, num_actv_constr,
     2148                                                     wj, num_actv_bnds, data->mem_mpn_i,
     2149                                                     data->mem_for_bounds_int, data->mem_for_bounds_real);
     2150
     2151                                if (fabs(values1[n+m+k]) > 1.0e3 * eq_tol)
     2152                                {
     2153                                    wj[j] = it;
     2154                                    wj[k] = QPP_AT_INACTIVE;
     2155                                }
     2156                                else
     2157                                {
     2158                                    ldPush(j, -1, it, m, n, wi, &num_actv_constr, wj,
     2159                                           &num_actv_bnds, y, z, values1+n, eq_tol);
     2160                                }
     2161                            }
     2162                            else
     2163                            {
     2164                                wj[j] = it;
     2165                            }
    19592166                        }
    1960                         wj[k] = QPP_AT_INACTIVE;
    1961                         printf("%" QPP_PRID "  ", k);
    1962 
    1963                         if (wi[i] != QPP_AT_INACTIVE)
    1964                             ++sup;
    1965 
    1966                         if (fabs(y[i]) <= opt->feas_tol)
    1967                             ++z0;
    1968 
    1969                         if (alpha > 0.0)
     2167
     2168                        if (y[i] >= -1.0e3 * eq_tol)
    19702169                        {
    1971                             wi[i] = QPP_AT_UPPER_BOUND;
     2170                            it = QPP_AT_LOWER_BOUND;
    19722171                        }
    19732172                        else
    19742173                        {
    1975                             wi[i] = QPP_AT_LOWER_BOUND;
     2174                            it = QPP_AT_UPPER_BOUND;
    19762175                        }
    19772176
    1978                         if (y[i] > opt->feas_tol)
     2177                        if (wi[i] == QPP_AT_INACTIVE)
    19792178                        {
    1980                             wi[i] = QPP_AT_LOWER_BOUND;
     2179                            values1 = data->mem_mpn_r;
     2180                            linearDependencyTest(values1, A, H, vec_indices_1, vec_values_1,
     2181                                                 vec_length_1, av, n_ps, wi, num_actv_constr,
     2182                                                 wj, num_actv_bnds, data->mem_mpn_i,
     2183                                                 data->mem_for_bounds_int, data->mem_for_bounds_real);
     2184
     2185                            if (fabs(values1[n+m+k]) > 1.0e3 * eq_tol)
     2186                            {
     2187                                wi[i] = it;
     2188                                ++num_actv_constr;
     2189
     2190                                wj[k] = QPP_AT_INACTIVE;
     2191                                --num_actv_bnds;
     2192                            }
     2193                            else
     2194                            {
     2195                                ldPush(i, 1, it, m, n, wi, &num_actv_constr, wj,
     2196                                       &num_actv_bnds, y, z, values1+n, eq_tol);
     2197                            }
    19812198                        }
    1982 
    1983                         if (y[i] < -opt->feas_tol)
     2199                        else
    19842200                        {
    1985                             wi[i] = QPP_AT_UPPER_BOUND;
     2201                            wi[i] = it;
    19862202                        }
    1987 
    1988                         /*if (vec_length_1-sup > 0)
    1989                         {
    1990                             for (j = 0; j < n; ++j)
    1991                             {
    1992                                 if (av[j] && (fabs(z[j]) <= opt->feas_tol) && (wj[j] != QPP_AT_INACTIVE) )
    1993                                 {
    1994                                     printf("\nIndex = %" QPP_PRID "\n", j);
    1995                                     wj[j] = QPP_AT_INACTIVE;
    1996                                     break;
    1997                                 }
    1998                             }
    1999                         }*/
    2000 
    2001                         printf("\n2 Number active components / total / Zero Multipliers: %d / %" QPP_PRID " / %d\n", sup, vec_length_1+1, z0);
    2002                     }
    2003                 }
    2004                 //wj[378] = QPP_AT_INACTIVE;
    2005                 //wj[835] = QPP_AT_INACTIVE;
    2006             }
    2007 
     2203                    }
     2204                }
     2205            }
    20082206
    20092207            /*free(indices1);
     
    20392237                    wi[i] = QPP_AT_UPPER_BOUND;
    20402238                }
     2239                ++num_actv_constr;
    20412240            }
    20422241
     
    20632262
    20642263            ac[i] = 1;
     2264            ++m_ps;
    20652265            y[i] = new_y;
    20662266
     
    20752275            if ( (wi != NULL) && (wj != NULL) )
    20762276            {
     2277                /*  All variables are already fixed at correct bound due to the postprocessing
     2278                    operation of fixVariable. */
    20772279                if (fabs(y[i]) <= eq_tol)
    20782280                {
    20792281                    wi[i] = QPP_AT_INACTIVE;
     2282                    --num_actv_constr;
    20802283                }
    20812284                else
     
    20862289                    for (jj = 0; jj < length1+length2; ++jj)
    20872290                    {
    2088                         if (fabs(z[indices1[jj]]) <= eq_tol)
     2291                        if (fabs(z[indices1[jj]]) <= 1.0e3 * eq_tol)
    20892292                        {
    20902293                            wj[indices1[jj]] = QPP_AT_INACTIVE;
     
    20942297                    }
    20952298                    assert(it >= 0);
     2299                    --num_actv_bnds;
    20962300                }
    20972301            }
     
    21022306            break;
    21032307
    2104         case QPP_SE_SPARSIFICATION:     /* Multipliers OK, working set not. */
     2308        case QPP_SE_SPARSIFICATION:     /* OK (with linear dependency test) */
    21052309            spf = se->elements.spf;
    21062310            alpha = spf->alpha;
     
    21102314            getRowA(A, i, av, vec_indices_1, vec_values_1, &vec_length_1);
    21112315
     2316            /*  !!! If the entries in A are NOT sorted, then put the following line
     2317                into the next for loop!!! */
     2318            it = A->irowst[k];
    21122319            for (jj = 0; jj < vec_length_1; ++jj)
    21132320            {
    21142321                j = vec_indices_1[jj];
    2115                 it = A->irowst[k];
    2116 
    2117                 /*  The following loop always works, but if the entries are sorted
    2118                     (as they should be) we  do not have to reset it = A->irowst[k]
    2119                     in every iteration */
     2322
    21202323                while (A->jcn[it] != j)
    21212324                {
     
    21392342            if ( (wi != NULL) && (wj != NULL) )
    21402343            {
    2141                 if (fabs(y[k]) <= 1e-8)
    2142                 {
    2143                     continue;
    2144                 }
    2145                 /* Problem if constr. k is active and constraint i is not!!! */
    2146                 /*if ( wi[k] && !wi[i] && (fabs(y[k]) > 1e-8) )
    2147                 {
    2148                     puts("Problem SPF.");
    2149                     printf("z_j[j] = %.6e\n", z[spf->j]);
    2150 
    2151                     if (y[i] <= 0)
     2344                if (wi[i] != QPP_AT_INACTIVE)
     2345                {
     2346                    if (y[i] >= -1.0e3 * eq_tol)
     2347                    {
     2348                        wi[i] = QPP_AT_LOWER_BOUND;
     2349                    }
     2350                    else
    21522351                    {
    21532352                        wi[i] = QPP_AT_UPPER_BOUND;
    21542353                    }
    2155                     else //if (y[i] > opt->feas_tol)
    2156                     {
    2157                         wi[i] = QPP_AT_LOWER_BOUND;
    2158                     }
    2159                 }*/
    2160 
    2161                 if (wi[i] != QPP_AT_INACTIVE)
    2162                 {
    2163                     if (y[i] <= -opt->feas_tol)
    2164                     {
    2165                         wi[i] = QPP_AT_UPPER_BOUND;
    2166                     }
    2167                     else if (y[i] > opt->feas_tol)
    2168                     {
    2169                         wi[i] = QPP_AT_LOWER_BOUND;
    2170                     }
    2171                 }
    2172 
    2173             }
    2174 
     2354                }
     2355
     2356                if ( (wi[i] == QPP_AT_INACTIVE) && (wi[k] != QPP_AT_INACTIVE) )
     2357                {
     2358                    values1 = data->mem_mpn_r;
     2359
     2360                    ret_val = linearDependencyTest(values1, A, H, vec_indices_1, vec_values_1,
     2361                                                   vec_length_1, av, n_ps, wi, num_actv_constr,
     2362                                                   wj, num_actv_bnds, data->mem_mpn_i,
     2363                                                   data->mem_for_bounds_int,
     2364                                                   data->mem_for_bounds_real);
     2365                    if (ret_val != QPP_OK)
     2366                    {
     2367                        return ret_val;
     2368                    }
     2369
     2370                    /* Constraint must be linearly dependent. For debugging we nevertheless check it. */
     2371                    /*value = 0.0;
     2372                    for (jj = 0; jj < n; ++jj)
     2373                    {
     2374                        value += values1[jj] * values1[jj];
     2375                    }
     2376
     2377                    assert(value <= eq_tol);*/
     2378
     2379                    kk = (alpha * y[k] > 1.0e3 * eq_tol ? QPP_AT_UPPER_BOUND : QPP_AT_LOWER_BOUND);
     2380
     2381                    ret_val = ldPush(i, 1, kk, m, n, wi, &num_actv_constr, wj, &num_actv_bnds, y, z,
     2382                                     values1+n, eq_tol);
     2383
     2384                    if (ret_val != QPP_OK)
     2385                    {
     2386                        return ret_val;
     2387                    }
     2388                }
     2389            }
    21752390
    21762391            /*free(spf);*/
     
    22012416                x[j] += sc->au / a_ij;  /* sc->au == sc->al */
    22022417                av[j] = 1;
     2418                ++n_ps;
    22032419                z[j] = -a_ij * y[i];
    22042420                /*++A_nnz_rows[i];*/
     
    22162432                if ( (wi != NULL) && (wj != NULL) )
    22172433                {
     2434                    assert(wi[i] != QPP_AT_EQ_CONSTRAINT);
     2435
    22182436                    /* Here, wi[i] stores the status of the i-th (inequality!) constraint */
    22192437                    if (wi[i] == QPP_AT_LOWER_BOUND)
    22202438                    {
     2439                        ++num_actv_bnds;
    22212440                        if (a_ij > 0.0)
    22222441                        {
     
    22302449                    else if (wi[i] == QPP_AT_UPPER_BOUND)
    22312450                    {
     2451                        ++num_actv_bnds;
    22322452                        if (a_ij > 0.0)
    22332453                        {
     
    22432463                        /* Inequality constraint is inactive and so is x_j. */
    22442464                        wj[j] = QPP_AT_INACTIVE;
    2245                     }
    2246                 }
    2247 
    2248                 /*  The original linear (eq.) constraint must be considered as active
    2249                     since then the dimension of the null space of the constraints stays
    2250                     the same. */
    2251                 if (wi != NULL)
    2252                 {
    2253                     if (y[i] >= 0.0)
     2465
     2466                        /*  i-th linear constraint will become active (if it is not
     2467                            already active). So, if the constraint is not active, we have
     2468                            to increase the number of active constraints. */
     2469                        ++num_actv_constr;
     2470                    }
     2471
     2472                    /* Active (equality) constraint i -> Dimension of null space stays the same! */
     2473                    if (y[i] >= -1.0e3 * eq_tol)
    22542474                    {
    22552475                        wi[i] = QPP_AT_LOWER_BOUND;
     
    22662486                av[j] = 1;
    22672487                ac[i] = 1;
     2488                ++n_ps;
     2489                ++m_ps;
    22682490                z[j] = 0.0;
    22692491
     
    22822504                        {
    22832505                            wi[i] = QPP_AT_LOWER_BOUND;
     2506                            ++num_actv_constr;
    22842507                        }
    22852508                    }
     
    22912514                        {
    22922515                            wi[i] = QPP_AT_UPPER_BOUND;
     2516                            ++num_actv_constr;
    22932517                        }
    22942518                    }
     
    23032527                        {
    23042528                            wi[i] = QPP_AT_UPPER_BOUND;
     2529                            ++num_actv_constr;
    23052530                        }
    23062531                    }
     
    23112536                        {
    23122537                            wi[i] = QPP_AT_LOWER_BOUND;
     2538                            ++num_actv_constr;
    23132539                        }
    23142540                    }
     
    38244050        if (H_nnz_columns[j] == 0)
    38254051        {
    3826             /* Solve LP: min g_j * x s.t. xl_j <= x <= xu_j */
     4052            /* Solve LP: min g_j * x_j s.t. xl_j <= x <= xu_j */
    38274053            if (g_j > QPP_ZERO)
    38284054            {
     
    38754101                if ( isfinite(xl_j) && isfinite(xu_j) )
    38764102                {
    3877                     value = (xu_j + xl_j) / 2;
    3878                     actv_type = QPP_AT_INACTIVE;
    3879                 }
    3880                 else if (isfinite(xl_j))
    3881                 {
     4103                    /* To ensure that the reduced Hessian stays pos. def.: Fix at lower bound. */
    38824104                    value = xl_j;
    38834105                    actv_type = QPP_AT_LOWER_BOUND;
    38844106                }
     4107                else if (isfinite(xl_j))
     4108                {
     4109                    value = xl_j;
     4110                    actv_type = QPP_AT_LOWER_BOUND;
     4111                }
    38854112                else if (isfinite(xu_j))
    38864113                {
     
    38904117                else
    38914118                {
    3892                     value = 0.0;
    3893                     actv_type = QPP_AT_INACTIVE;
     4119                    /*  Here, x_j is a free variable, does neither occur in linear constraints,
     4120                        nor in the objective function. We can not remove it here, as the reduced
     4121                        Hessian might become indefinite (not pos. def.). */
     4122                    continue;
     4123                    /*value = 0.0;
     4124                    actv_type = QPP_AT_INACTIVE;*/
    38944125                }
    38954126            }
     
    39124143            /* Solve QP: min 0.5 * h_jj * x^2 + g_j * x s.t. xl_j <= x <= xu_j */
    39134144            value = -g_j / h_jj;    /* Solution of unconstrained problem */
     4145            actv_type = QPP_AT_INACTIVE;
    39144146
    39154147            if (h_jj > 0.0)
     
    41734405        if (need_to_split == 1)
    41744406        {
     4407            //continue;
    41754408            /* Use split technique only after 2 iterations of the presolving loop. */
    41764409            if ((data->num_iter <= 2) || !qppIsEqual(al[i], au[i], eq_tol))
     
    61146347    return QPP_OK;
    61156348}
     6349
     6350
     6351static cholmod_dense *solveKKTSystem(const qpp_ecrmatrix_t *A,
     6352                                     const qpp_ecrmatrix_t *H,
     6353                                     cholmod_dense *rhs,
     6354                                     const qpp_int_t *av,
     6355                                     const qpp_int_t num_av,
     6356                                     const qpp_int_t *wi,
     6357                                     const qpp_int_t num_actv_constr,
     6358                                     const qpp_int_t *wj,
     6359                                     const qpp_int_t num_actv_bnds,
     6360                                     qpp_int_t *storage,
     6361                                     qpp_int_t *vec_indices,
     6362                                     qpp_real_t *vec_values)
     6363{
     6364    qpp_int_t i, j, k, l, m, n, num_entries, vec_length, dim_it, kkt_dim, status;
     6365    qpp_real_t *Kx, *Ksol_x, *rhs_x;
     6366    double umf_info[UMFPACK_INFO], umf_control[UMFPACK_CONTROL];
     6367    qpp_int_t *Ki, *Kp, *av_offset, *wi_offset;
     6368    cholmod_sparse *K;
     6369    cholmod_dense *Ksol;
     6370    cholmod_common cc;
     6371    void *symbolic, *numeric;
     6372
     6373    m = A->nrow;
     6374    n = A->ncol;
     6375    av_offset = storage;
     6376    qppSetArrayi(av_offset, m + 2*n, 0);
     6377    wi_offset = &av_offset[n];
     6378    K = NULL;
     6379    Ksol = NULL;
     6380
     6381    kkt_dim = num_av + num_actv_constr + num_actv_bnds;
     6382
     6383    /*qpp_int_t n_av = 0;
     6384    qpp_int_t n_wi = 0;
     6385    qpp_int_t n_wj = 0;
     6386
     6387    for (k = 0; k < n; ++k)
     6388    {
     6389        if (av[k])
     6390            ++n_av;
     6391        if (wj[k])
     6392            ++n_wj;
     6393    }
     6394    for (k = 0; k < m; ++k)
     6395    {
     6396        if (wi[k] )
     6397            ++n_wi;
     6398    }
     6399
     6400    printf("REF n_av = %d, n_wi = %d, n_wj = %d\n", (int)(num_av-n_av),
     6401           (int)(num_actv_constr-n_wi), (int)(num_actv_bnds-n_wj));*/
     6402
     6403    if (!av[0])
     6404    {
     6405        av_offset[0] = 1;
     6406    }
     6407    for (k = 1; k < n; ++k)
     6408    {
     6409        av_offset[k] = av_offset[k-1];
     6410        if (!av[k])
     6411        {
     6412            ++av_offset[k];
     6413        }
     6414    }
     6415
     6416    if (!wi[0])
     6417    {
     6418        wi_offset[0] = 1;
     6419    }
     6420    for (k = 1; k < m; ++k)
     6421    {
     6422        wi_offset[k] = wi_offset[k-1];
     6423        if (!wi[k])
     6424        {
     6425            ++wi_offset[k];
     6426        }
     6427    }
     6428
     6429    CM_START(&cc);
     6430    K = CM_ALLOCATE_SPARSE(kkt_dim, kkt_dim, 2 * (H->nnz+A->nnz+num_actv_bnds), 1, 1, 0,
     6431                           CHOLMOD_REAL, &cc);
     6432    CM_FINISH(&cc);
     6433
     6434    if (K == NULL)
     6435    {
     6436        return NULL;
     6437    }
     6438
     6439    num_entries = 0;
     6440    dim_it = 0;
     6441    Ki = K->i;
     6442    Kp = K->p;
     6443    Kx = K->x;
     6444
     6445    Kp[dim_it++] = 0;
     6446    l = 0;
     6447
     6448    /* Build KKT matrix in CSC format */
     6449    for (j = 0; j < n; ++j)
     6450    {
     6451        if (av[j])
     6452        {
     6453            getColH(H, j, av, vec_indices, vec_values, &vec_length);
     6454            for (k = 0; k < vec_length; ++k)
     6455            {
     6456                Ki[num_entries] = vec_indices[k] - av_offset[vec_indices[k]];
     6457                Kx[num_entries++] = vec_values[k];
     6458            }
     6459
     6460            getColA(A, j, wi, vec_indices, vec_values, &vec_length);
     6461            for (k = 0; k < vec_length; ++k)
     6462            {
     6463                Ki[num_entries] = vec_indices[k] - wi_offset[vec_indices[k]] + num_av;
     6464                Kx[num_entries++] = vec_values[k];
     6465            }
     6466
     6467            if (wj[j])
     6468            {
     6469                Ki[num_entries] = l + num_av + num_actv_constr;
     6470                Kx[num_entries++] = 1.0;
     6471                ++l;
     6472            }
     6473            Kp[dim_it++] = num_entries;
     6474        }
     6475    }
     6476
     6477    for (i = 0; i < m; ++i)
     6478    {
     6479        if (wi[i])
     6480        {
     6481            getRowA(A, i, av, vec_indices, vec_values, &vec_length);
     6482            for (k = 0; k < vec_length; ++k)
     6483            {
     6484                Ki[num_entries] = vec_indices[k] - av_offset[vec_indices[k]];
     6485                Kx[num_entries++] = vec_values[k];
     6486            }
     6487            Kp[dim_it++] = num_entries;
     6488        }
     6489    }
     6490
     6491    for (j = 0; j < n; ++j)
     6492    {
     6493        if (wj[j])
     6494        {
     6495            Ki[num_entries] = j - av_offset[j];
     6496            Kx[num_entries++] = 1.0;
     6497            Kp[dim_it++] = num_entries;
     6498        }
     6499    }
     6500
     6501    /*FILE *mat_file = fopen("KKT.mm", "w");
     6502    CM_WRITE_SPARSE(mat_file, K, NULL, NULL, &cc);
     6503    fclose(mat_file);*/
     6504
     6505    /* SOLVE WITH UMFPACK */
     6506    UMFPACK_DEFAULTS(umf_control);
     6507    status = UMFPACK_SYMBOLIC(kkt_dim, kkt_dim, Kp, Ki, Kx, &symbolic, umf_control, umf_info);
     6508    if (status < 0)
     6509    {
     6510        CM_START(&cc);
     6511        CM_FREE_SPARSE(&K, &cc);
     6512        CM_FINISH(&cc);
     6513        return NULL;
     6514    }
     6515
     6516    status = UMFPACK_NUMERIC(Kp, Ki, Kx, symbolic, &numeric, umf_control, umf_info);
     6517    if (status < 0)
     6518    {
     6519        CM_START(&cc);
     6520        CM_FREE_SPARSE(&K, &cc);
     6521        CM_FINISH(&cc);
     6522        return NULL;
     6523    }
     6524
     6525    rhs_x = rhs->x;
     6526    Ksol = CM_ALLOCATE_DENSE(kkt_dim, 1, kkt_dim, CHOLMOD_REAL, &cc);
     6527    Ksol_x = Ksol->x;
     6528
     6529    status = UMFPACK_SOLVE(UMFPACK_A, Kp, Ki, Kx, Ksol_x, rhs_x, numeric, umf_control, umf_info);
     6530    if (status < 0)
     6531    {
     6532        CM_START(&cc);
     6533        CM_FREE_SPARSE(&K, &cc);
     6534        CM_FREE_DENSE(&Ksol, &cc);
     6535        CM_FINISH(&cc);
     6536        return NULL;
     6537    }
     6538
     6539    /*CM_START(&cc);
     6540    Ksol = SPQR_BACKSLASH(K, rhs, &cc);
     6541    CM_FINISH(&cc);*/
     6542
     6543    if (Ksol == NULL)
     6544    {
     6545        puts("Keine Loesung");
     6546        CM_START(&cc);
     6547        CM_FREE_SPARSE(&K, &cc);
     6548        CM_FINISH(&cc);
     6549        return NULL;
     6550    }
     6551
     6552    UMFPACK_FREE_NUMERIC(&numeric);
     6553    UMFPACK_FREE_SYMBOLIC(&symbolic);
     6554
     6555    CM_START(&cc);
     6556    CM_FREE_SPARSE(&K, &cc);
     6557    CM_FINISH(&cc);
     6558
     6559    return Ksol;
     6560}
     6561
     6562
     6563static qpp_return_value_t curvatureTest(qpp_real_t *result,
     6564                                        const qpp_int_t index,
     6565                                        const qpp_ecrmatrix_t *A,
     6566                                        const qpp_ecrmatrix_t *H,
     6567                                        const qpp_int_t *av,
     6568                                        const qpp_int_t num_av,
     6569                                        const qpp_int_t *wi,
     6570                                        const qpp_int_t num_actv_constr,
     6571                                        const qpp_int_t *wj,
     6572                                        const qpp_int_t num_actv_bnds,
     6573                                        qpp_int_t *storage,
     6574                                        qpp_int_t *vec_indices,
     6575                                        qpp_real_t *vec_values)
     6576{
     6577    qpp_int_t m, n, k, kkt_dim, counter;
     6578    qpp_real_t *x;
     6579    cholmod_dense *rhs, *Ksol;
     6580    cholmod_common cc;
     6581
     6582    m = A->nrow;
     6583    n = A->ncol;
     6584    rhs = NULL;
     6585    Ksol = NULL;
     6586
     6587    kkt_dim = num_av + num_actv_constr + num_actv_bnds;
     6588
     6589    CM_START(&cc);
     6590
     6591    rhs = CM_ALLOCATE_DENSE(kkt_dim, 1, kkt_dim, CHOLMOD_REAL, &cc);
     6592
     6593    CM_FINISH(&cc);
     6594
     6595    if (rhs == NULL)
     6596    {
     6597        return QPP_OUT_OF_MEMORY;
     6598    }
     6599
     6600    x = rhs->x;
     6601    qppSetArrayr(x, kkt_dim, 0.0);
     6602
     6603    /*  If index \in {0,...,m-1}, then curvature test with linear constraint <index>,
     6604        else curvature test with bound constraint <index> - <m>. */
     6605    if (index < m)
     6606    {
     6607        counter = 0;
     6608        for (k = 0; k < n; ++k)
     6609        {
     6610            if (!wi[k])
     6611            {
     6612                ++counter;
     6613            }
     6614
     6615            if (k == index)
     6616            {
     6617                x[num_av + k - counter] = 1.0;
     6618                break;
     6619            }
     6620        }
     6621    }
     6622    else
     6623    {
     6624        counter = 0;
     6625        for (k = 0; k < n; ++k)
     6626        {
     6627            if (!wj[k])
     6628            {
     6629                ++counter;
     6630            }
     6631
     6632            if (k == (index-m))
     6633            {
     6634                x[num_av + num_actv_constr + k - counter] = 1.0;
     6635                break;
     6636            }
     6637        }
     6638    }
     6639
     6640    Ksol = solveKKTSystem(A, H, rhs, av, num_av, wi, num_actv_constr, wj, num_actv_bnds,
     6641                          storage, vec_indices, vec_values);
     6642    if (Ksol == NULL)
     6643    {
     6644        CM_START(&cc);
     6645        CM_FREE_DENSE(&rhs, &cc);
     6646        CM_FINISH(&cc);
     6647        return QPP_OUT_OF_MEMORY;
     6648    }
     6649
     6650    x = Ksol->x;
     6651
     6652    qppSetArrayr(result, m + 2*n, 0.0);
     6653    counter = 0;
     6654    for (k = 0; k < n; ++k)
     6655    {
     6656        if (av[k])
     6657        {
     6658            result[k] = x[counter++];
     6659        }
     6660    }
     6661
     6662    assert(counter == num_av);
     6663    counter = num_av;
     6664
     6665    for (k = 0; k < m; ++k)
     6666    {
     6667        if (wi[k])
     6668        {
     6669            result[n+k] = x[counter++];
     6670        }
     6671    }
     6672
     6673    assert(counter == num_av + num_actv_constr);
     6674    counter = num_av + num_actv_constr;
     6675
     6676    for (k = 0; k < n; ++k)
     6677    {
     6678        if (wj[k])
     6679        {
     6680            result[n+m+k] = x[counter++];
     6681        }
     6682    }
     6683
     6684    CM_START(&cc);
     6685
     6686    CM_FREE_DENSE(&rhs, &cc);
     6687    CM_FREE_DENSE(&Ksol, &cc);
     6688
     6689    CM_FINISH(&cc);
     6690
     6691    return QPP_OK;
     6692}
     6693
     6694
     6695static qpp_return_value_t linearDependencyTest(qpp_real_t *result,
     6696                                               const qpp_ecrmatrix_t *A,
     6697                                               const qpp_ecrmatrix_t *H,
     6698                                               const qpp_int_t *rhs_indices,
     6699                                               const qpp_real_t *rhs_values,
     6700                                               const qpp_int_t rhs_length,
     6701                                               const qpp_int_t *av,
     6702                                               const qpp_int_t num_av,
     6703                                               const qpp_int_t *wi,
     6704                                               const qpp_int_t num_actv_constr,
     6705                                               const qpp_int_t *wj,
     6706                                               const qpp_int_t num_actv_bnds,
     6707                                               qpp_int_t *storage,
     6708                                               qpp_int_t *vec_indices,
     6709                                               qpp_real_t *vec_values)
     6710{
     6711    qpp_int_t m, n, counter, k, it, kkt_dim;
     6712    qpp_real_t *x;
     6713    cholmod_dense *rhs, *Ksol;
     6714    cholmod_common cc;
     6715
     6716    m = A->nrow;
     6717    n = A->ncol;
     6718    kkt_dim = num_av + num_actv_constr + num_actv_bnds;
     6719    rhs = NULL;
     6720    Ksol = NULL;
     6721
     6722    CM_START(&cc);
     6723
     6724    rhs = CM_ALLOCATE_DENSE(kkt_dim, 1, kkt_dim, CHOLMOD_REAL, &cc);
     6725
     6726    CM_FINISH(&cc);
     6727
     6728    if (rhs == NULL)
     6729    {
     6730        return QPP_OUT_OF_MEMORY;
     6731    }
     6732
     6733    x = rhs->x;
     6734    qppSetArrayr(x, kkt_dim, 0.0);
     6735
     6736    it = 0;
     6737    counter = 0;
     6738    for (k = 0; k < n; ++k)
     6739    {
     6740        if (!av[k])
     6741        {
     6742            ++counter;
     6743        }
     6744
     6745        if ( (it < rhs_length) && (k == rhs_indices[it]) )
     6746        {
     6747            x[k-counter] = rhs_values[it];
     6748            ++it;
     6749        }
     6750
     6751    }
     6752    assert((n-counter)<=num_av);
     6753
     6754    Ksol = solveKKTSystem(A, H, rhs, av, num_av, wi, num_actv_constr, wj, num_actv_bnds,
     6755                          storage, vec_indices, vec_values);
     6756    if (Ksol == NULL)
     6757    {
     6758        CM_START(&cc);
     6759        CM_FREE_DENSE(&rhs, &cc);
     6760        CM_FINISH(&cc);
     6761        return QPP_OUT_OF_MEMORY;
     6762    }
     6763
     6764    x = Ksol->x;
     6765
     6766    qppSetArrayr(result, m + 2*n, 0.0);
     6767    counter = 0;
     6768    for (k = 0; k < n; ++k)
     6769    {
     6770        if (av[k])
     6771        {
     6772            result[k] = x[counter++];
     6773        }
     6774    }
     6775
     6776    assert(counter == num_av);
     6777    counter = num_av;
     6778
     6779    for (k = 0; k < m; ++k)
     6780    {
     6781        if (wi[k])
     6782        {
     6783            result[n+k] = x[counter++];
     6784        }
     6785    }
     6786
     6787    assert(counter == num_av + num_actv_constr);
     6788    counter = num_av + num_actv_constr;
     6789
     6790    for (k = 0; k < n; ++k)
     6791    {
     6792        if (wj[k])
     6793        {
     6794            result[n+m+k] = x[counter++];
     6795        }
     6796    }
     6797
     6798    CM_START(&cc);
     6799
     6800    CM_FREE_DENSE(&rhs, &cc);
     6801    CM_FREE_DENSE(&Ksol, &cc);
     6802
     6803    CM_FINISH(&cc);
     6804
     6805    return QPP_OK;
     6806}
     6807
     6808
     6809static qpp_return_value_t ldPush(const qpp_int_t index,
     6810                                 const qpp_int_t sign,
     6811                                 const qpp_int_t ws,
     6812                                 qpp_int_t m,
     6813                                 qpp_int_t n,
     6814                                 qpp_int_t *wi,
     6815                                 qpp_int_t *num_wi,
     6816                                 qpp_int_t *wj,
     6817                                 qpp_int_t *num_wj,
     6818                                 qpp_real_t *y,
     6819                                 qpp_real_t *z,
     6820                                 qpp_real_t *scalars,
     6821                                 qpp_real_t tol)
     6822{
     6823    qpp_real_t factor;
     6824    qpp_int_t w;
     6825    qpp_real_t multiplier;
     6826
     6827    qpp_int_t sgn = sign;
     6828    w = index;
     6829
     6830    assert(sign != 0);
     6831
     6832    if (sign > 0)
     6833        multiplier = y[index];
     6834    else
     6835        multiplier = z[index];
     6836
     6837    factor = 1.0;
     6838
     6839    qpp_int_t i, j;
     6840
     6841    for (i = 0; i < m; ++i)
     6842    {
     6843        if ( (wi[i] == QPP_AT_LOWER_BOUND) && (scalars[i]*multiplier < -1.0e3 * tol) )
     6844        {
     6845            if (factor > -y[i] / (scalars[i] * multiplier))
     6846            {
     6847                sgn = 1;
     6848                w = i;
     6849                factor = -y[i] / (scalars[i] * multiplier);
     6850            }
     6851        }
     6852        else if ( (wi[i] == QPP_AT_UPPER_BOUND) && (scalars[i]*multiplier > 1.0e3 * tol) )
     6853        {
     6854            if (factor > -y[i] / (scalars[i] * multiplier))
     6855            {
     6856                sgn = 1;
     6857                w = i;
     6858                factor = -y[i] / (scalars[i] * multiplier);
     6859            }
     6860        }
     6861    }
     6862
     6863    scalars += m;
     6864    for (j = 0; j < n; ++j)
     6865    {
     6866        if ( (wj[j] == QPP_AT_LOWER_BOUND) && (scalars[j]*multiplier < -1.0e3 * tol) )
     6867        {
     6868            if (factor > -z[j] / (scalars[j] * multiplier))
     6869            {
     6870                sgn = -1;
     6871                w = j;
     6872                factor = -z[j] / (scalars[j] * multiplier);
     6873            }
     6874        }
     6875        else if ( (wj[j] == QPP_AT_UPPER_BOUND) && (scalars[j]*multiplier > 1.0e3 * tol) )
     6876        {
     6877            if (factor > -z[j] / (scalars[j] * multiplier))
     6878            {
     6879                sgn = -1;
     6880                w = j;
     6881                factor = -z[j] / (scalars[j] * multiplier);
     6882            }
     6883        }
     6884    }
     6885
     6886    scalars -= m;
     6887
     6888    for (i = 0; i < m; ++i)
     6889    {
     6890        if (wi[i] != QPP_AT_INACTIVE)
     6891        {
     6892            y[i] += factor * multiplier * scalars[i];
     6893        }
     6894    }
     6895
     6896    for (j = 0; j < n; ++j)
     6897    {
     6898        if (wj[j] != QPP_AT_INACTIVE)
     6899        {
     6900            z[j] += factor * multiplier * scalars[m+j];
     6901        }
     6902    }
     6903
     6904    if (sign > 0)
     6905        y[index] *= (1.0 - factor);
     6906    else
     6907        z[index] *= (1.0 - factor);
     6908
     6909    if (sign > 0)
     6910    {
     6911        wi[index] = ws;
     6912        (*num_wi) += 1;
     6913    }
     6914    else
     6915    {
     6916        wj[index] = ws;
     6917        (*num_wj) += 1;
     6918    }
     6919
     6920    assert(sgn != 0);
     6921
     6922    if (sgn > 0)
     6923    {
     6924        wi[w] = QPP_AT_INACTIVE;
     6925        (*num_wi) -= 1;
     6926    }
     6927    else
     6928    {
     6929        wj[w] = QPP_AT_INACTIVE;
     6930        (*num_wj) -= 1;
     6931    }
     6932
     6933    return QPP_OK;
     6934}
  • branches/presolver/src/Makefile

    r260 r261  
    3737
    3838IFLAGS      =  -I. \
    39                -I${IDIR} \
    40                -I${PRESOLVER_IDIR}
     39               -I${IDIR}
     40
     41ifeq (${NO_PRESOLVER}, 0)
     42IFLAGS += \
     43    -I${PRESOLVER_IDIR} \
     44    -I${SUITESPARSE_IDIR}
     45endif
    4146
    4247QPOASES_OBJECTS = \
     
    5459        Matrices.${OBJEXT} \
    5560        MessageHandling.${OBJEXT} \
    56         SparseSolver.${OBJEXT} \
    57         Presolver.${OBJEXT} \
    58         PresolverOptions.${OBJEXT}
     61        SparseSolver.${OBJEXT} 
    5962
     63ifeq (${NO_PRESOLVER}, 0)
     64QPOASES_OBJECTS += \
     65    Presolver.${OBJEXT} \
     66    PresolverOptions.${OBJEXT}
     67endif
    6068
    6169QPOASES_EXTRAS_OBJECTS = \
     
    7987        ${IDIR}/qpOASES/Matrices.hpp \
    8088        ${IDIR}/qpOASES/MessageHandling.hpp \
    81         ${IDIR}/qpOASES/UnitTesting.hpp \
    82         ${IDIR}/qpOASES/Presolver.hpp \
     89        ${IDIR}/qpOASES/UnitTesting.hpp
     90       
     91ifeq (${NO_PRESOLVER}, 0)
     92QPOASES_DEPENDS += \
     93    ${IDIR}/qpOASES/Presolver.hpp \
    8394        ${IDIR}/qpOASES/PresolverOptions.hpp
    84 
     95endif
    8596
    8697##
  • branches/presolver/src/Presolver.cpp

    r260 r261  
    9696        return THROWERROR( convertErrorCode(err) );
    9797    }
     98
    9899    return SUCCESSFUL_RETURN;
    99100}
     
    155156
    156157        for (int_t i = 0; i < *nzH; ++i)
    157             Hcp[Hjcn[i]+1]++;
     158            ++Hcp[Hjcn[i]+1];
    158159
    159160        for (int_t i = 1; i <= nVar; ++i)
     
    167168
    168169        for (int_t i = 0; i < *nzA; ++i)
    169             Acp[Ajcn[i]+1]++;
     170            ++Acp[Ajcn[i]+1];
    170171
    171172        for (int_t i = 1; i <= nVar; ++i)
     
    649650    }
    650651
     652    /* Extract working set regarding bound constraints if necessary. */
    651653    if (bounds != 0)
    652654    {
     
    673675                /* Other status is currently not allowed. */
    674676                delete[] wj;
    675                 return RET_INVALID_ARGUMENTS;
     677                return THROWERROR( RET_INVALID_ARGUMENTS );
    676678            }
    677679        }
    678680    }
    679681
     682    /* Extract working set regarding linear constraints if necessary. */
    680683    if (constraints != 0)
    681684    {
     
    702705                delete[] wi;
    703706                delete[] wj;
    704                 return RET_INVALID_ARGUMENTS;
     707                return THROWERROR( RET_INVALID_ARGUMENTS );
    705708            }
    706709        }
     
    711714        if (err != QPP_OK)
    712715        {
     716            delete[] wi;
     717            delete[] wj;
    713718                return THROWERROR( convertErrorCode(err) );
    714719        }
    715720
     721        /* Create new qpOASES working set regarding bound constraints. */
    716722        if (bounds != 0)
    717723    {
     
    760766    }
    761767
     768    /* Create new qpOASES working set regarding linear constraints. */
    762769    if (constraints != 0)
    763770    {
  • branches/presolver/src/PresolverOptions.cpp

    r260 r261  
    120120        char boundInfo[32];
    121121
    122         if (bound_mode == PBT_MEDIUM)
     122        if (options::bound_mode == PBT_MEDIUM)
    123123    {
    124124        strncpy(boundInfo, "Medium Bounds", 32);
  • branches/presolver/testing/c/Makefile

    r239 r261  
    3131
    3232
     33NO_PRESOLVER = 1
    3334
    3435include ../../make.mk
  • branches/presolver/testing/cpp/Makefile

    r239 r261  
    3030##
    3131
     32NO_PRESOLVER = 1
    3233
    3334include ../../make.mk
Note: See TracChangeset for help on using the changeset viewer.