QLDDirect.h

Go to the documentation of this file.

/*
 * Copyright 2012-2020 CNRS-UM LIRMM, CNRS-AIST JRL
 */
 
#pragma once
 
// includes
// std
#include <cassert>
 
// Eigen
#include <Eigen/Core>
 
// eigen-qld
#include "eigen_qld_api.h"
 
namespace Eigen
{
 
class QLDDirect
{
public:
  EIGEN_QLD_API QLDDirect();
 
  EIGEN_QLD_API QLDDirect(int nrvar, int nreq, int nrineq, int ldq = -1, int lda = -1, bool verbose = false);
 
  EIGEN_QLD_API void fdOut(int fd);
  EIGEN_QLD_API int fdOut() const;
 
  EIGEN_QLD_API void verbose(bool v);
  EIGEN_QLD_API bool verbose() const;
 
  EIGEN_QLD_API int fail() const;
 
  EIGEN_QLD_API void problem(int nrvar, int nreq, int nrineq, int ldq = -1, int lda = -1);
 
  EIGEN_QLD_API const VectorXd & result() const;
 
  EIGEN_QLD_API const VectorXd & multipliers() const;
 
  template<typename MatObj, typename VecObj, typename MatConstr, typename VecConstr, typename VecVar>
  bool solve(const MatrixBase<MatObj> & Q,
             const MatrixBase<VecObj> & c,
             const MatrixBase<MatConstr> & A,
             const MatrixBase<VecConstr> & b,
             const MatrixBase<VecVar> & xl,
             const MatrixBase<VecVar> & xu,
             int nreq,
             bool isDecomp = false,
             double eps = 1e-12);
 
private:
  EIGEN_QLD_API int fortran_ql(const int * m,
                               const int * me,
                               const int * mmax,
                               const int * n,
                               const int * nmax,
                               const int * mnn,
                               const double * c,
                               const double * d,
                               const double * a,
                               const double * b,
                               const double * xl,
                               const double * xu,
                               double * x,
                               double * u,
                               const double * eps,
                               const int * mode,
                               const int * iout,
                               int * ifail,
                               const int * iprint,
                               double * war,
                               int * lwar,
                               int * iwar,
                               int * liwar);
 
  VectorXd X_;
  int fdOut_;
  int verbose_;
  int fail_;
  VectorXd U_;
  VectorXd WAR_;
  VectorXi IWAR_;
};
 
template<typename MatObj, typename VecObj, typename MatConstr, typename VecConstr, typename VecVar>
inline bool QLDDirect::solve(const MatrixBase<MatObj> & Q,
                             const MatrixBase<VecObj> & c,
                             const MatrixBase<MatConstr> & A,
                             const MatrixBase<VecConstr> & b,
                             const MatrixBase<VecVar> & xl,
                             const MatrixBase<VecVar> & xu,
                             int nreq,
                             bool isDecomp,
                             double eps)
{
  assert(A.rows() == b.rows()); // check constraint size
  assert(A.cols() == X_.rows());
  assert(Q.rows() == Q.cols()); // check Q is square
  assert(Q.cols() == X_.rows()); // check Q has the good number of variable
  assert(c.rows() == X_.rows()); // check C size
  assert(xl.rows() == X_.rows()); // check XL size
  assert(xu.rows() == X_.rows()); // check XU size
 
  int mode = isDecomp ? 0 : 1;
 
  int M = int(A.rows());
  int N = int(X_.rows());
 
  int MMAX = std::max(int(A.stride()), 1);
  int NMAX = std::max(int(Q.stride()), 1);
 
  int NMN = M + 2 * N;
  int LWAR = int(WAR_.rows());
  int LIWAR = int(IWAR_.rows());
 
  assert(LWAR >= (3. * NMAX * NMAX) / 2. + 10. * NMAX + MMAX + M + 1.
         && "Please call QLD::problem with the correct dimensions.");
  assert(LIWAR >= N && "Please call QLD::problem with the correct dimensions.");
 
  fortran_ql(&M, &nreq, &MMAX, &N, &NMAX, &NMN, Q.derived().data(), c.derived().data(), A.derived().data(),
             b.derived().data(), xl.derived().data(), xu.derived().data(), X_.data(), U_.data(), &eps, &mode, &fdOut_,
             &fail_, &verbose_, WAR_.data(), &LWAR, IWAR_.data(), &LIWAR);
 
  return fail_ == 0;
}
 
} // namespace Eigen