linear_complementarity_problem.h

/* Copyright (C) 2004-2012 MBSim Development Team
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA
 *
 * Contact: martin.o.foerg@googlemail.com
 */
 
#ifndef NUMERICS_LINEAR_COMPLEMENTARITY_PROBLEM_H_
#define NUMERICS_LINEAR_COMPLEMENTARITY_PROBLEM_H_
 
#include <fmatvec/fmatvec.h>
 
#include <mbsim/numerics/nonlinear_algebra/multi_dimensional_newton_method.h>
#include <mbsim/numerics/nonlinear_algebra/multi_dimensional_fixpoint_solver.h>
#include <mbsim/numerics/linear_complementarity_problem/lemke_algorithm.h>
 
#include <mbsim/numerics/functions/lcp_reformulation_functions.h>
 
namespace MBSim {
 
  class LinearComplementarityProblem : public virtual fmatvec::Atom {
    public:
      enum LCPSolvingStrategy {
        Standard,      //trying to solve the LCP in a standard way (first some steps of Lemke algorithm, then trying iterative schemes in reformulated system, then trying Lemke again as fallback with all steps)
        ReformulatedStandard,  //trying to use the iterative schemes at first (as fallback use Lemke)
        ReformulatedNewtonOnly,
        ReformulatedFixpointOnly,
        LemkeOnly
          //only use LemkeAgorithm
      };
 
      enum JacobianType {
        Numerical,  //use only a numerical Jacobian Matrix
        LCPSpecial
          //use a analytical Jacobian Matrix almost everywhere and choose a derivation at the kink (=LinearComplementarityJacobianFunction) //TODO: find a better name
      };
 
      enum CriteriaType {
        Global,
        Local
      };
 
      LinearComplementarityProblem(const fmatvec::SymMat & M_, const fmatvec::Vec & q_, const LCPSolvingStrategy & strategy_ = Standard, const JacobianType & jacobianType_ = LCPSpecial);
 
      virtual ~LinearComplementarityProblem();
 
      /*GETTER / SETTER*/
      void setSystem(const fmatvec::SymMat & M_, const fmatvec::Vec & q_);
      void setStrategy(const LCPSolvingStrategy & strategy_) {
        strategy = strategy_;
        oldStrategy = strategy_;
      }
      void setJacobianType(const JacobianType & jacobianType_) {
        jacobianType = jacobianType_;
      }
      void setNewtonCriteriaFunction(CriteriaFunction * criteriaFunction_) {
        newtonSolver->setCriteriaFunction(criteriaFunction_);
      }
      void setFixpointCriteriaFunction(CriteriaFunction * criteriaFunction_) {
        fixpointSolver->setCriteriaFunction(criteriaFunction_);
      }
      /****************/
 
      fmatvec::Vec solve(const fmatvec::Vec & initialSolution = fmatvec::Vec(0, fmatvec::NONINIT));
 
      static double computeMediumEigVal(const fmatvec::SqrMat & M);
 
      static fmatvec::Vec createInitialSolution(const fmatvec::SymMat & M, const fmatvec::Vec & q, double mediumEigVal = 0);
 
      static fmatvec::Vec createInitialSolution(const fmatvec::SqrMat & M, const fmatvec::Vec & q, double mediumEigVal = 0);
 
    protected:
      void useNewton(fmatvec::Vec & solution, bool & solved);
 
      void useFixpoint(fmatvec::Vec & solution, bool & solved);
 
      fmatvec::SqrMat M;
 
      fmatvec::Vec q;
 
      LCPSolvingStrategy strategy;
 
      LCPSolvingStrategy oldStrategy;
 
      double mediumEigenValue;
 
      JacobianType jacobianType;
 
      LemkeAlgorithm lemkeSolver;
 
      MBSim::MultiDimensionalNewtonMethod * newtonSolver;
 
      LCPNewtonReformulationFunction * newtonFunction;
 
      NewtonJacobianFunction * jacobianFunction;
 
      CriteriaFunction * criteriaNewton;
 
      MultiDimensionalFixpointSolver * fixpointSolver;
 
      LCPFixpointReformulationFunction * fixpointFunction;
 
      CriteriaFunction * criteriaFixedpoint;
  };
 
  std::ostream& operator <<(std::ostream & o, const LinearComplementarityProblem::LCPSolvingStrategy &strategy);
 
  std::ostream& operator <<(std::ostream & o, const LinearComplementarityProblem::JacobianType &jacobianType);
 
}
 
#endif /*__LINEAR_COMPLEMENTARITY_PROBLEM_H_*/