SOCLaw.hpp

/*
 * This file is part of So-bogus, a C++ sparse block matrix library and
 * Second Order Cone solver.
 *
 * Copyright 2013 Gilles Daviet <gdaviet@gmail.com>
 *
 * So-bogus is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.

 * So-bogus 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 General Public License for more details.

 * You should have received a copy of the GNU General Public License
 * along with So-bogus.  If not, see <http://www.gnu.org/licenses/>.
*/


#ifndef BOGUS_SOCLAW_HPP
#define BOGUS_SOCLAW_HPP

#include "FischerBurmeister.hpp"
#include "LocalSOCSolver.hpp"

#include <vector>

namespace bogus
{


template < DenseIndexType Dimension, typename Scalar, bool DeSaxceCOV,
           local_soc_solver::Strategy Strat >
class SOCLaw
{
public:
    typedef LocalProblemTraits< Dimension, Scalar > Traits ;
    enum{ dimension = Dimension } ;


    SOCLaw( const unsigned n, const double * mu ) ;

    Scalar eval( const unsigned problemIndex,
                 const typename Traits::Vector &x,
                 const typename Traits::Vector &y ) const
    {
        typedef FischerBurmeister< Traits::dimension, typename Traits::Scalar, DeSaxceCOV > FBFunction ;

        if(m_mu[problemIndex] < 0) return y.squaredNorm() ;

        typename Traits::Vector fb( x.rows() ) ;
        FBFunction::compute( m_mu[problemIndex], x, y, fb ) ;

        return fb.squaredNorm() ;
    }


    bool solveLocal(
            const unsigned problemIndex,
            const typename Traits::Matrix &A,
            const typename Traits::Vector &b,
            typename Traits::Vector &x,
            const Scalar scaling
            ) const ;

    void projectOnConstraint( const unsigned problemIndex, typename Traits::Vector &x ) const ;


    template< typename Segment >
    void dualityCOV( const unsigned problemIndex, const Segment& y,
                     typename Traits::Vector& s ) const
    {
        if( DeSaxceCOV ) {
            Traits::np( s ) = m_mu[problemIndex] * Traits::tp(y).norm() ;
            Traits::tp( s ).setZero() ;
        } else  s.setZero() ;
    }

private:

    const double * m_mu ;
    const unsigned m_n ;
    Scalar m_localTol ;

} ;

typedef SOCLaw< 2, double,  true > Coulomb2D ;
typedef SOCLaw< 3, double,  true > Coulomb3D ;
typedef SOCLaw< 2, double, false > SOC2D ;
typedef SOCLaw< 3, double, false > SOC3D ;

}

#endif // SOCLAW_HPP