triangulation.cpp

Go to the documentation of this file.

/*****************************************************************************
    
    triangulation.cpp -- See declarations in triangulation.hpp.

    This file is a part of the Arageli library.

    Copyright (C) Sergey S. Lyalin, 2006
    University of Nizhni Novgorod, Russia, 2005--2006

*****************************************************************************/

#include "config.hpp"

#if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || \
    defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_triangulation)

#include <stack>

#include "submatrix/indexed.hpp"
#include "vector.hpp"

#include "triangulation.hpp"


namespace Arageli
{

template
<
    typename Q,
    typename Dim1,
    typename TR,
    typename Dim2,
    typename Ctrler
>
void triangulate_simple_1
(
    const Q& q,
    const Dim1& dim,
    TR& tr,
    const Dim2& subspdim,
    Ctrler ctrler
)
{
    typedef typename Q::size_type Idx;
    typedef indexed_submatrix<const Q> SubQ;

    SubQ curq(&q, fromall);
    vector<Idx, false> cursimplex;
    Idx curineq = 0;
    vector<Idx, false> erasedrows;
    Idx curdim = dim;
    
    // stack of pairs (number of inequation, corresponding q matrix)
    typedef std::pair<Idx, SubQ> IdxAndSubQ;
    std::stack<IdxAndSubQ> stackq;
    
    for(;;)
    {
        if(curq.is_empty())
        {
            // Add subsimplexes from a subspace (if any).
            for(Idx i = 0; i <= subspdim; ++i)
            {
                for(Idx j = 0; j <= subspdim; ++j)
                    if(i != j)cursimplex.push_back(j + q.nrows());
                if(tr.ncols() < cursimplex.size())
                    tr.assign();
                if(tr.ncols() <= cursimplex.size())
                    tr.insert_row(tr.nrows(), cursimplex);
                cursimplex.erase(cursimplex.size() - subspdim, subspdim);
            }

            curineq = curq.ncols();
        }

        for(; curineq < curq.ncols() && is_null(curq(0, curineq)); ++curineq);

        if(curineq < curq.ncols())
        {
            erasedrows.push_back(0);
            
            for(Idx i = 1; i < curq.nrows(); ++i)
                if(!is_null(curq(i, curineq)))
                    erasedrows.push_back(i);

            if(erasedrows.size() < curq.nrows() || curq.nrows() == 1)
            {
                // Reduction of the matrix q to the next facet
                // that is not incidented to the first ray.

                cursimplex.push_back(curq.row_index()[0]);
                stackq.push(IdxAndSubQ(curineq, curq));
                curq.row_index().erase_subvector(erasedrows);
                curq.col_index().erase(curineq);
                curineq = 0;
                --curdim;

                for(Idx j = 0; j < curq.ncols();)
                {
                    Idx zeros = 0;
                    for(Idx i = 0; i < curq.nrows(); ++i)
                        zeros += is_null(curq(i, j));

                    if(zeros < curdim-1)curq.col_index().erase(j);
                    else ++j;
                } 
            }
            else
            {
                ++curineq;
            }

            erasedrows.assign();
        }
        else
        {
            if(stackq.empty())break;    // end of triangulation process
            curineq = stackq.top().first + 1;   // move to the next inequation
            curq = stackq.top().second;
            stackq.pop();
            ARAGELI_ASSERT_1(!cursimplex.is_empty());
            cursimplex.pop_back();
            ++curdim;
        }
    }
}



}


#else   // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...


namespace Arageli
{

// PLACE ALL NOT TEMPLATE NOT INLINE IMPLEMENTATIONS HERE

}


#endif  // #if !defined(ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE) || ...

---