basis_inverse.h

CGAL is a collaborative effort of several sites in Europe and Israel. The goal is to make the most i

// Copyright (c) 1997-2001  ETH Zurich (Switzerland).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $Source: /CVSROOT/CGAL/Packages/_QP_solver/include/CGAL/_QP_solver/Basis_inverse.h,v $
// $Revision: 1.7 $ $Date: 2004/09/03 17:41:06 $
// $Name:  $
//
// Author(s)     : Sven Schoenherr <[e-mail]>
                                                                               

#ifndef CGAL_BASIS_INVERSE_H
#define CGAL_BASIS_INVERSE_H

// includes
#include <CGAL/Optimisation/basic.h>
#include <vector>
#include <iterator>
#include <CGAL/IO/Verbose_ostream.h>


CGAL_BEGIN_NAMESPACE
                    

// Class declaration
// =================
template < class ET_, class Is_lp_ >
class Basis_inverse;
                    

template < class ET, class Is_lp >
class Basis_inverse__entry_iterator;


// Class interface
// ===============
template < class ET_, class Is_lp_ >
class Basis_inverse {
  public:
    // self
    typedef  ET_                        ET;
    typedef  Is_lp_                     Is_lp;
    typedef  Basis_inverse<ET,Is_lp>    Self;

  private:
    
    // private types
    typedef  std::vector<ET>            Row;
    typedef  std::vector<Row>           Matrix;
    
    
    typedef  CGAL::Tag_true             Tag_true;
    typedef  CGAL::Tag_false            Tag_false;
    
    
    // friends
    friend  class Basis_inverse__entry_iterator<ET,Is_lp>;
    
    

  public:
    
    // types
    typedef  Basis_inverse__entry_iterator<ET,Is_lp>
                                        Entry_iterator;
                                                       

    
    // creation
    Basis_inverse( CGAL::Verbose_ostream&  verbose_out)
        : et_0( 0), et_1( 1), vout( verbose_out) { }
    
    
    // initialization
    template < class InputIterator1, class InputIterator2 >
    void
    init( unsigned int qp_m, unsigned int l,
          InputIterator1 u_it, InputIterator2 w_it,
          unsigned int max_size = 0)
        {
            CGAL_optimisation_precondition( qp_m > 0);
            m = qp_m;
            d = et_1;
            init( l, u_it, w_it, max_size, Is_lp());
            
            CGAL_optimisation_debug {
                if ( vout.verbose()) {
                    for ( unsigned int row = 0; row < k; ++row) {
                        std::copy( M[ row].begin(), M[ row].end(),
                                   std::ostream_iterator<ET>( vout.out()," "));
                        vout.out() << std::endl;
                    }
                    vout.out() << "denominator = " << d << std::endl;
                }
            }
             
        }
    
    
    // access
    const ET&  denominator( ) const { return d; }
    
    
    Entry_iterator  column_begin( unsigned int j) const
        { return Entry_iterator( M, 0, j); }
    Entry_iterator  column_end  ( unsigned int j) const
        { return Entry_iterator( M, k, j); }
    
    
    // multiplication functions
    template < class ForwardIterator, class OutputIterator >
    void
    multiply( ForwardIterator z_l, ForwardIterator z_x,
               OutputIterator y_l,  OutputIterator y_x) const
        {
            multiply( z_l, z_x, y_l, y_x, Is_lp());
        }
    
    
    // special multiplication functions for LPs
    template < class ForwardIterator, class OutputIterator >
    void
    multiply_l( ForwardIterator z_x, OutputIterator y_l) const
        {
            multiply_l( z_x, y_l, Is_lp());
        }
    
    template < class ForwardIterator, class OutputIterator >
    void
    multiply_x( ForwardIterator z_l, OutputIterator y_x) const
        {
            multiply_x( z_l, y_x, Is_lp());
        }
    
    
    // update functions
    /*
    template < class ForwardIterator > inline
    void  append( ForwardIterator q_l, ForwardIterator q_x, const ET& nu);
    
    
    inline
    void  remove( unsigned int i);
    
    
    template < class RandomAccessIterator > inline
    void  replace( unsigned int i, RandomAccessIterator q_x);
    */
    
    // swap function
    /*
    inline
    void  swap( unsigned int, unsigned int);
    */
    

  private:
    
    // some constants
    const ET                 et_0, et_1;
                                        

    
    // data members
    unsigned int             m;         // number of constraints
    unsigned int             k;         // size of matrix
    Matrix                   M;         // basis inverse, stored row-wise
    ET                       d;         // denominator
    
    
    CGAL::Verbose_ostream&   vout;      // used for verbose output
    
    

    
    // initialization
    /*
    template < class InputIterator1, class InputIterator2 >
    inline  void  init( unsigned int l,
                        InputIterator1 u_it, InputIterator2 w_it,
                        unsigned int max_size, Tag_false);
    template < class InputIterator1, class InputIterator2 >
    inline  void  init( unsigned int l,
                        InputIterator1 u_it, InputIterator2 w_it,
                        unsigned int max_size, Tag_true );
    */
    
    // multiplication functions
    /*
    template < class ForIt, class OutIt > inline                        // QP
    void  multiply( ForIt z_l, ForIt z_x,
                    OutIt y_l, OutIt y_x, Tag_false) const;
    template < class ForIt, class OutIt > inline                        // LP
    void  multiply( ForIt z_l, ForIt z_x,
                    OutIt y_l, OutIt y_x, Tag_true ) const;
    */
    
    // special multiplication functions for LPs
    /*
    template < class ForIt, class OutIt > inline                        // QP
    void  multiply_l( ForIt z_x, OutIt y_l, Tag_false) const;
    template < class ForIt, class OutIt > inline                        // LP
    void  multiply_l( ForIt z_x, OutIt y_l, Tag_true ) const;
    
    template < class ForIt, class OutIt > inline                        // QP
    void  multiply_x( ForIt z_l, OutIt y_x, Tag_false) const;
    template < class ForIt, class OutIt > inline                        // LP
    void  multiply_x( ForIt z_l, OutIt y_x, Tag_true ) const;
    */
  

  

// ============================================================================
                                                                               

// Class Implementation
// ====================

// initialization
// --------------

// initialization (QP)
template < class InputIterator1, class InputIterator2 > inline
void
init( unsigned int l,
      InputIterator1 u_it, InputIterator2 w_it,
      unsigned int max_size, Tag_false)
{
    k = 2*m;
    M.erase( M.begin(), M.end());
    M.reserve( max_size > m ? m+max_size : k+1);
    unsigned int  i;
    for ( i = 0; i < k; ) M.push_back( Row( ++i, et_0));
    for ( i = 0; i < m; ++i, ++u_it, ++w_it) {
        M[ m+i][ l] = *w_it;
        M[ m+i][ i] = *u_it;
    }
}

// initialization (LP)
template < class InputIterator1, class InputIterator2 > inline
void
init( unsigned int l,
      InputIterator1 u_it, InputIterator2 w_it,
      unsigned int, Tag_true)
{
    k = m;
    M = Matrix( m, Row( m, et_0));
    for ( unsigned int i = 0; i < m; ++i, ++u_it, ++w_it) {
        M[ i][ l] = *w_it;
        M[ i][ i] = *u_it;
    }
}


// multiplication functions
// ------------------------

// multiply (QP)
template < class ForIt, class OutIt > inline
void
multiply( ForIt z_l, ForIt z_x, OutIt y_l, OutIt y_x, Tag_false) const
{
    typename Matrix::const_iterator  matrix_it = M.begin();
    typename Row   ::const_iterator     row_it;     // left  of diagonal
    typename Matrix::const_iterator  column_it;     // right of diagonal
    ForIt                                 z_it;

    unsigned int  row, count;
    ET            sum;

    // compute  P z_l + Q^T z_x
    for ( row = 0; row < m; ++row,                                ++y_l) {
        sum = et_0;

        // P: left of diagonal (inclusive)
        for (   row_it =  matrix_it->begin(),                z_it = z_l;
                row_it != matrix_it->end();
              ++row_it,                                    ++z_it      ) {
            sum += *row_it * *z_it;
        }

        // P: right of diagonal (exclusive)
        for (   count = row+1,   column_it = ++matrix_it;
                count < m;
              ++count,         ++column_it,                ++z_it      ) {
            sum += (*column_it)[ row] * *z_it;
        }

        // Q^T:
        for (                                                z_it = z_x;
                count < k;
              ++count,         ++column_it,                ++z_it      ) {
            sum += (*column_it)[ row] * *z_it;
        }

        // store result
        *y_l = sum;
    }

    // compute  Q z_l + R z_x
    for ( ; row < k; ++row,                                       ++y_x) {
        sum = et_0;

        // Q:
        for (   count = 0,   row_it =  matrix_it->begin(),   z_it = z_l;
                count < m;
              ++count,     ++row_it,                       ++z_it      ) {
            sum += *row_it * *z_it;
        }

        // R: left of diagonal (inclusive)
        for (                                                z_it = z_x;
                             row_it != matrix_it->end();
                           ++row_it,                       ++z_it      ) {
            sum += *row_it * *z_it;
        }

        // R: right of diagonal (exclusive)
        for (   count = row+1,   column_it = ++matrix_it;
                count < k;
              ++count,         ++column_it,                ++z_it      ) {
            sum += (*column_it)[ row] * *z_it;
        }

        // store result
        *y_x = sum;
    }
}

// multiply (LP)
template < class ForIt, class OutIt > inline
void
multiply( ForIt z_l, ForIt z_x, OutIt y_l, OutIt y_x, Tag_true) const
{
    multiply_l( z_x, y_l);
    multiply_x( z_l, y_x);
}

// multiply_l (QP)
template < class ForIt, class OutIt > inline
void
multiply_l( ForIt z_x, OutIt y_l, Tag_false) const
{
    typename Matrix::const_iterator  matrix_it = M.begin()+m;
    typename Matrix::const_iterator  column_it;
    ForIt                                 z_it;

    unsigned int  row, count;
    ET            sum;

    // compute  Q^T z_x
    for ( row = 0; row < m; ++row,                                ++y_l) {
        sum = et_0;

        for (   count = 0,   column_it = matrix_it,   z_it = z_x;
                count < m;
              ++count,     ++column_it,             ++z_it      ) {
            sum += (*column_it)[ row] * *z_it;
        }

        *y_l = sum;
    }
}

// multiply_x (QP)
template < class ForIt, class OutIt > inline
void
multiply_x( ForIt z_l, OutIt y_x, Tag_false) const
{
    typename Matrix::const_iterator  matrix_it = M.begin()+m;
    typename Row   ::const_iterator     row_it;
    ForIt                                 z_it;

    unsigned int  row, count;
    ET            sum;

    // compute  Q z_l
    for ( row = 0; row < m; ++row,  ++matrix_it,                  ++y_x) {
        sum = et_0;

        for (   count = 0,   row_it = matrix_it->begin(),   z_it = z_l;
                count < m;
              ++count,     ++row_it,                      ++z_it      ) {
            sum += *row_it * *z_it;
        }

        *y_x = sum;
    }
}

// multiply_l (LP)
template < class ForIt, class OutIt > inline
void
multiply_l( ForIt z_x, OutIt y_l, Tag_true) const
{
    typename Matrix::const_iterator  matrix_it;
    ForIt                                 z_it;

    unsigned int  row;
    ET            sum;

    // compute  Q^T z_x
    for ( row = 0; row < m; ++row,                   ++y_l) {
        sum = et_0;

        for (   matrix_it =  M.begin(),   z_it = z_x;
                matrix_it != M.end();
              ++matrix_it,              ++z_it      ) {