gauss.hpp

Go to the documentation of this file.

/*****************************************************************************
    
    gauss.hpp -- Some algorithms of linear algebra.

    This file is part of Arageli library.

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

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

#ifndef _ARAGELI_GAUSS_HPP_
#define _ARAGELI_GAUSS_HPP_


#include "config.hpp"

#include <cmath>
#include <iostream>
#include <iomanip>
#include <functional>

#include "type_opers.hpp"
#include "type_traits.hpp"
#include "factory.hpp"
#include "exception.hpp"
#include "powerest.hpp"
#include "vector.hpp"
#include "matrix.hpp"

#include "std_import.hpp"

namespace Arageli
{



template <typename M>
void gauss_field_row_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col
);



template <typename M>
void gauss_field_col_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col
);



template <typename M, typename T>
void gauss_bareiss_row_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col,
    T& det, T& det_denom
);



template <typename M, typename T>
void gauss_bareiss_col_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col,
    T& det, T& det_denom
);



template <typename M>
void gauss_bareiss_row_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col
)
{
    typename M::element_type det, det_denom;
    gauss_bareiss_row_iter(a, row, col, det, det_denom);
}



template <typename M>
void gauss_bareiss_col_iter
(
    M& a,
    typename M::size_type row,
    typename M::size_type col
)
{
    typename M::element_type det, det_denom;
    gauss_bareiss_col_iter(a, row, col, det, det_denom);
}


namespace ctrl
{

struct rref_gauss_field_idler
{
    class abort : public ctrl::abort {};

    template <typename A>
    void preamble (const A&) const {}

    template <typename B, typename Q, typename Det, typename Basis>
    void conclusion (const B&, const Q&, const Det&, const Basis&) const {}

    template <typename B, typename Q, typename Det, typename Basis>
    void before_iter (const B&, const Q&, const Det&, const Basis&) const {}

    template <typename B, typename Q, typename Det, typename Basis>
    void after_iter (const B&, const Q&, const Det&, const Basis&) const {}

    template <typename J>
    void find_biggest_in_col (const J&) const {}

    template <typename J>
    void negligible_col (const J&) const {}

    template <typename I, typename J>
    void pivot_item (const I&, const J& j) const {}

    template <typename I1, typename I2, typename B, typename Q, typename Det>
    void swap_rows (const I1& i1, const I2& i2, const B&, const Q&, const Det&) const {}

    template <typename J>
    void eliminate_col (const J&) const {}

    template <typename B, typename Q, typename Det>
    void after_elimination (const B&, const Q&, const Det&) const {}
};


struct rref_gauss_bareiss_idler : public rref_gauss_field_idler
{
    class abort : public ctrl::abort {};

    template
    <
        typename B, typename Q, typename Det,
        typename Basis, typename Det_denom
    >
    void before_iter
    (const B&, const Q&, const Det&, const Basis&, const Det_denom&) const {}

    template
    <
        typename B, typename Q, typename Det,
        typename Basis, typename Det_denom
    >
    void after_iter
    (const B&, const Q&, const Det&, const Basis&, const Det_denom&) const {}

    template
    <
        typename I1, typename I2, typename B,
        typename Q, typename Det, typename Det_denom
    >
    void swap_rows
    (const I1&, const I2&, const B&, const Q&, const Det&, const Det_denom&) const {}

    template <typename B, typename Q, typename Det, typename Det_denom>
    void after_elimination
    (const B&, const Q&, const Det&, const Det_denom&) const {}

    template <typename J>
    void find_nonzero_in_col (const J&) const {}
};


struct rref_field_idler
{
    class abort : public ctrl::abort {};

    rref_gauss_field_idler ctrl_rref_gauss_field () const
    { return rref_gauss_field_idler(); }
};


struct rref_int_idler
{
    class abort : public ctrl::abort {};

    rref_gauss_bareiss_idler ctrl_rref_gauss_bareiss () const
    { return rref_gauss_bareiss_idler(); }
};


struct rref_idler
{
    class abort : public ctrl::abort {};

    rref_field_idler ctrl_rref_field () const
    { return rref_field_idler(); }

    rref_int_idler ctrl_rref_int () const
    { return rref_int_idler(); }
};


} // namespace ctrl



template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
void rref_gauss_field
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler
);


template
<
    typename A, typename B, typename Q, typename Basis, typename T,
    typename Ctrler // control visitor
>
inline void rref_gauss_field
(const A& a, B& b, Q& q, Basis& basis, T& det, Ctrler ctrler)
{ rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrler); }


template <typename A, typename B, typename Q, typename Basis, typename T>
inline void rref_gauss_field (const A& a, B& b, Q& q, Basis& basis, T& det)
{ rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_field_idler()); }


template <typename A, typename B, typename Q, typename Basis>
inline void rref_gauss_field (const A& a, B& b, Q& q, Basis& basis)
{
    typedef typename B::element_type T;
    T det;
    rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_field_idler());
}


template <typename A, typename B, typename Q>
inline void rref_gauss_field (const A& a, B& b, Q& q)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_field_idler());
}


template <typename A, typename B>
inline void rref_gauss_field (const A& a, B& b)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    B q;
    rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_field_idler());
}


template <typename A>
inline A rref_gauss_field (const A& a)
{
    typedef typename A::element_type T;
    T det;
    vector<typename A::size_type, false> basis;
    A b, q;
    rref_gauss_field(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_field_idler());
    return b;
}

// ----------------------------------------


template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
inline void rref_field
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler
)
{ rref_gauss_field(a, b, q, basis, det, tfctr, ctrler.ctrl_rref_gauss_field()); }


template
<
    typename A, typename B, typename Q, typename Basis, typename T,
    typename Ctrler // control visitor
>
inline void rref_field
(const A& a, B& b, Q& q, Basis& basis, T& det, Ctrler ctrler)
{ rref_field(a, b, q, basis, det, factory<T>(), ctrler); }


template <typename A, typename B, typename Q, typename Basis, typename T>
inline void rref_field (const A& a, B& b, Q& q, Basis& basis, T& det)
{ rref_field(a, b, q, basis, det, factory<T>(), ctrl::rref_field_idler()); }


template <typename A, typename B, typename Q, typename Basis>
inline void rref_field (const A& a, B& b, Q& q, Basis& basis)
{
    typedef typename B::element_type T;
    T det;
    rref_field(a, b, q, basis, det, factory<T>(), ctrl::rref_field_idler());
}


template <typename A, typename B, typename Q>
inline void rref_field (const A& a, B& b, Q& q)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    rref_field(a, b, q, basis, det, factory<T>(), ctrl::rref_field_idler());
}


template <typename A, typename B>
inline void rref_field (const A& a, B& b)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    B q;
    rref_field(a, b, q, basis, det, factory<T>(), ctrl::rref_field_idler());
}


template <typename A>
inline A rref_field (const A& a)
{
    typedef typename A::element_type T;
    T det;
    vector<typename A::size_type> basis;
    A b, q;
    rref_field(a, b, q, basis, det, factory<T>(), ctrl::rref_field_idler());
    return b;
}


// ----------------------------------------

template <typename A1, typename A2, typename I>
void mat_col_copy_order
(const A1& a1, A2& a2, const I& idx);

template <typename A1, typename A2, typename I>
void vec_copy_order
(const A1& a1, A2& a2, const I& idx);


template <typename A, typename B>
void vec_inverse_permutation (const A& a, B& b);



// TODO Make this function with an algorithm controller.
template
<
    typename A, typename B, typename Q,
    typename Basis_in, typename Basis_out,
    typename Det
>
void rref_order
(
    const A& a, B& b, Q& q,
    const Basis_in& basis_in, Basis_out& basis_out,
    Det& det
);



template <typename T, bool REFCNT>
matrix<T, REFCNT> inverse (const matrix<T, REFCNT>& A);



template <typename T, bool REFCNT>
T det (const matrix<T, REFCNT>& A);


template <typename T, bool REFCNT>
typename matrix<T, REFCNT>::size_type rank
(const matrix<T, REFCNT>& A);




template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
void rref_gauss_bareiss
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler
);


template
<
    typename A, typename B, typename Q, typename Basis, typename T,
    typename Ctrler // control visitor
>
inline void rref_gauss_bareiss
(const A& a, B& b, Q& q, Basis& basis, T& det, Ctrler ctrler)
{ rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrler); }


template <typename A, typename B, typename Q, typename Basis, typename T>
inline void rref_gauss_bareiss (const A& a, B& b, Q& q, Basis& basis, T& det)
{ rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_bareiss_idler()); }


template <typename A, typename B, typename Q, typename Basis>
inline void rref_gauss_bareiss (const A& a, B& b, Q& q, Basis& basis)
{
    typedef typename B::element_type T;
    T det;
    rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_bareiss_idler());
}


template <typename A, typename B, typename Q>
inline void rref_gauss_bareiss (const A& a, B& b, Q& q)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_bareiss_idler());
}


template <typename A, typename B>
inline void rref_gauss_bareiss (const A& a, B& b)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    B q;
    rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_bareiss_idler());
}


template <typename A>
inline A rref_gauss_bareiss (const A& a)
{
    typedef typename A::element_type T;
    T det;
    vector<typename A::size_type, false> basis;
    A b, q;
    rref_gauss_bareiss(a, b, q, basis, det, factory<T>(), ctrl::rref_gauss_bareiss_idler());
    return b;
}


// ----------------------------------------


template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
inline void rref_int
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler
)
{ rref_gauss_bareiss(a, b, q, basis, det, tfctr, ctrler.ctrl_rref_gauss_bareiss()); }


template
<
    typename A, typename B, typename Q, typename Basis, typename T,
    typename Ctrler // control visitor
>
inline void rref_int
(const A& a, B& b, Q& q, Basis& basis, T& det, Ctrler ctrler)
{ rref_int(a, b, q, basis, det, factory<T>(), ctrler); }


template <typename A, typename B, typename Q, typename Basis, typename T>
inline void rref_int (const A& a, B& b, Q& q, Basis& basis, T& det)
{ rref_int(a, b, q, basis, det, factory<T>(), ctrl::rref_int_idler()); }


template <typename A, typename B, typename Q, typename Basis>
inline void rref_int (const A& a, B& b, Q& q, Basis& basis)
{
    typedef typename B::element_type T;
    T det;
    rref_int(a, b, q, basis, det, factory<T>(), ctrl::rref_int_idler());
}


template <typename A, typename B, typename Q>
inline void rref_int (const A& a, B& b, Q& q)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    rref_int(a, b, q, basis, det, factory<T>(), ctrl::rref_int_idler());
}


template <typename A, typename B>
inline void rref_int (const A& a, B& b)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    B q;
    rref_int(a, b, q, basis, det, factory<T>(), ctrl::rref_int_idler());
}


template <typename A>
inline A rref_int (const A& a)
{
    typedef typename A::element_type T;
    T det;
    vector<typename A::size_type, false> basis;
    A b, q;
    rref_int(a, b, q, basis, det, factory<T>(), ctrl::rref_int_idler());
    return b;
}


// ----------------------------------------

namespace _Internal
{

template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
inline void rref_helper_1
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler,
    const true_type&
)
{ rref_field(a, b, q, basis, det, tfctr, ctrler.ctrl_rref_field()); }


template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
inline void rref_helper_1
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler,
    const false_type&
)
{ rref_int(a, b, q, basis, det, tfctr, ctrler.ctrl_rref_int()); }


}



template
<
    typename A, typename B, typename Q, typename Basis,
    typename T, typename T_factory,
    typename Ctrler // control visitor
>
inline void rref
(
    const A& a, B& b, Q& q, Basis& basis,
    T& det,
    const T_factory& tfctr,
    Ctrler ctrler
)
{
    _Internal::rref_helper_1
    (
        a, b, q, basis, det, tfctr, ctrler,
        bool_type<type_traits<typename B::element_type>::is_field>::value
    );
}


template
<
    typename A, typename B, typename Q, typename Basis, typename T,
    typename Ctrler // control visitor
>
inline void rref
(const A& a, B& b, Q& q, Basis& basis, T& det, Ctrler ctrler)
{ rref(a, b, q, basis, det, factory<T>(), ctrler); }


template <typename A, typename B, typename Q, typename Basis, typename T>
inline void rref (const A& a, B& b, Q& q, Basis& basis, T& det)
{ rref(a, b, q, basis, det, factory<T>(), ctrl::rref_idler()); }


template <typename A, typename B, typename Q, typename Basis>
void rref (const A& a, B& b, Q& q, Basis& basis)
{
    typedef typename B::element_type T;
    T det;
    rref(a, b, q, basis, det, factory<T>(), ctrl::rref_idler());
}


template <typename A, typename B, typename Q>
inline void rref (const A& a, B& b, Q& q)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    rref(a, b, q, basis, det, factory<T>(), ctrl::rref_idler());
}


template <typename A, typename B>
inline void rref (const A& a, B& b)
{
    typedef typename B::element_type T;
    T det;
    vector<typename B::size_type, false> basis;
    B q;
    rref(a, b, q, basis, det, factory<T>(), ctrl::rref_idler());
}


template <typename A>
inline A rref (const A& a)
{
    typedef typename A::element_type T;
    T det;
    vector<typename A::size_type, false> basis;
    A b, q;
    rref(a, b, q, basis, det, factory<T>(), ctrl::rref_idler());
    return b;
}


// ----------------------------------------

template <typename T, bool REFCNT>
T det_int (const matrix<T, REFCNT>& A);


template <typename T, bool REFCNT>
typename matrix<T, REFCNT>::size_type rank_int
(const matrix<T, REFCNT>& A);


// The following solve_linsys function have migrated to solve_linsys.{hpp, cpp} files.
//template
//<
//  typename T_A, bool REFCNT_A,
//  typename T_B, bool REFCNT_B,
//  typename T_res, bool REFCNT_res
//>
//inline void solve_linsys
//(
//  const matrix<T_A, REFCNT_A>& A,
//  const vector<T_B, REFCNT_B>& B,
//  vector<T_res, REFCNT_res>& res
//)
//{
//  // WARNING! Slow implementation.
//  res = inverse(A)*B;
//}
//
//
//template
//<
//  typename T_A, bool REFCNT_A,
//  typename T_B, bool REFCNT_B
//>
//inline vector<T_B, REFCNT_B> solve_linsys
//(
//  const matrix<T_A, REFCNT_A>& A,
//  const vector<T_B, REFCNT_B>& B
//)
//{
//  vector<T_B, REFCNT_B> res(B.size());
//  solve_linsys(A, B, res);
//  return res;
//}


}


#ifdef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE
    #define ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_GAUSS
    #include "gauss.cpp"
    #undef  ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE_GAUSS
#endif


#endif

---