linear_algebra.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 1,075 行 · 第 1/3 页

// Boost.Units - A C++ library for zero-overhead dimensional analysis and 
// unit/quantity manipulation and conversion
//
// Copyright (C) 2003-2008 Matthias Christian Schabel
// Copyright (C) 2008 Steven Watanabe
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP
#define BOOST_UNITS_DETAIL_LINEAR_ALGEBRA_HPP

#include <boost/units/static_rational.hpp>
#include <boost/mpl/next.hpp>
#include <boost/mpl/arithmetic.hpp>
#include <boost/mpl/and.hpp>
#include <boost/mpl/assert.hpp>

#include <boost/units/dim.hpp>
#include <boost/units/dimensionless_type.hpp>
#include <boost/units/static_rational.hpp>
#include <boost/units/detail/dimension_list.hpp>
#include <boost/units/detail/sort.hpp>

namespace boost {

namespace units {

namespace detail {

// typedef list<rational> equation;

template<int N>
struct eliminate_from_pair_of_equations_impl;

template<class E1, class E2>
struct eliminate_from_pair_of_equations;

template<int N>
struct elimination_impl;

template<bool is_zero, bool element_is_last>
struct elimination_skip_leading_zeros_impl;

template<class Equation, class Vars>
struct substitute;

template<int N>
struct substitute_impl;

template<bool is_end>
struct solve_impl;

template<class T>
struct solve;

template<int N>
struct check_extra_equations_impl;

template<int N>
struct normalize_units_impl;

struct inconsistent {};

// generally useful utilies.

template<int N>
struct divide_equation {
    template<class Begin, class Divisor>
    struct apply {
        typedef list<typename mpl::divides<typename Begin::item, Divisor>::type, typename divide_equation<N - 1>::template apply<typename Begin::next, Divisor>::type> type;
    };
};

template<>
struct divide_equation<0> {
    template<class Begin, class Divisor>
    struct apply {
        typedef dimensionless_type type;
    };
};

// eliminate_from_pair_of_equations takes a pair of
// equations and eliminates the first variable.
//
// equation eliminate_from_pair_of_equations(equation l1, equation l2) {
//     rational x1 = l1.front();
//     rational x2 = l2.front();
//     return(transform(pop_front(l1), pop_front(l2), _1 * x2 - _2 * x1));
// }

template<int N>
struct eliminate_from_pair_of_equations_impl {
    template<class Begin1, class Begin2, class X1, class X2>
    struct apply {
        typedef list<
            typename mpl::minus<
                typename mpl::times<typename Begin1::item, X2>::type,
                typename mpl::times<typename Begin2::item, X1>::type
            >::type,
            typename eliminate_from_pair_of_equations_impl<N - 1>::template apply<
                typename Begin1::next,
                typename Begin2::next,
                X1,
                X2
            >::type
        > type;
    };
};

template<>
struct eliminate_from_pair_of_equations_impl<0> {
    template<class Begin1, class Begin2, class X1, class X2>
    struct apply {
        typedef dimensionless_type type;
    };
};

template<class E1, class E2>
struct eliminate_from_pair_of_equations {
    typedef E1 begin1;
    typedef E2 begin2;
    typedef typename eliminate_from_pair_of_equations_impl<(E1::size::value - 1)>::template apply<
        typename begin1::next,
        typename begin2::next,
        typename begin1::item,
        typename begin2::item
    >::type type;
};



// Stage 1.  Determine which dimensions should
// have dummy base units.  For this purpose
// row reduce the matrix.

template<int N>
struct make_zero_vector {
    typedef list<static_rational<0>, typename make_zero_vector<N - 1>::type> type;
};
template<>
struct make_zero_vector<0> {
    typedef dimensionless_type type;
};

template<int Column, int TotalColumns>
struct create_row_of_identity {
    typedef list<static_rational<0>, typename create_row_of_identity<Column - 1, TotalColumns - 1>::type> type;
};
template<int TotalColumns>
struct create_row_of_identity<0, TotalColumns> {
    typedef list<static_rational<1>, typename make_zero_vector<TotalColumns - 1>::type> type;
};
template<int Column>
struct create_row_of_identity<Column, 0> {
    // error
};

template<int RemainingRows>
struct determine_extra_equations_impl;

template<bool first_is_zero, bool is_last>
struct determine_extra_equations_skip_zeros_impl;

// not the last row and not zero.
template<>
struct determine_extra_equations_skip_zeros_impl<false, false> {
    template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
    struct apply {
        // remove the equation being eliminated against from the set of equations.
        typedef typename determine_extra_equations_impl<RemainingRows - 1>::template apply<typename RowsBegin::next, typename RowsBegin::item>::type next_equations;
        // since this column was present, strip it out.
        typedef Result type;
    };
};

// the last row but not zero.
template<>
struct determine_extra_equations_skip_zeros_impl<false, true> {
    template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
    struct apply {
        // remove this equation.
        typedef dimensionless_type next_equations;
        // since this column was present, strip it out.
        typedef Result type;
    };
};


// the first columns is zero but it is not the last column.
// continue with the same loop.
template<>
struct determine_extra_equations_skip_zeros_impl<true, false> {
    template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class Result>
    struct apply {
        typedef typename RowsBegin::item current_row;
        typedef typename determine_extra_equations_skip_zeros_impl<
            current_row::item::Numerator == 0,
            RemainingRows == 2  // the next one will be the last.
        >::template apply<
            typename RowsBegin::next,
            RemainingRows - 1,
            CurrentColumn,
            TotalColumns,
            Result
        > next;
        typedef list<typename RowsBegin::item::next, typename next::next_equations> next_equations;
        typedef typename next::type type;
    };
};

// all the elements in this column are zero.
template<>
struct determine_extra_equations_skip_zeros_impl<true, true> {
    template<class RowsBegin, int RemainingRows, int CurrentColumn, int TotalColumns, class MatrixWithFirstColumnStripped, class Result>
    struct apply {
        typedef list<typename RowsBegin::item::next, dimensionless_type> next_equations;
        typedef list<typename create_row_of_identity<CurrentColumn, TotalColumns>::type, Result> type;
    };
};

template<int RemainingRows>
struct determine_extra_equations_impl {
    template<class RowsBegin, class EliminateAgainst>
    struct apply {
        typedef list<
            typename eliminate_from_pair_of_equations<typename RowsBegin::item, EliminateAgainst>::type,
            typename determine_extra_equations_impl<RemainingRows-1>::template apply<typename RowsBegin::next, EliminateAgainst>::type
        > type;
    };
};

template<>
struct determine_extra_equations_impl<0> {
    template<class RowsBegin, class EliminateAgainst>
    struct apply {
        typedef dimensionless_type type;
    };
};

template<int RemainingColumns, bool is_done>
struct determine_extra_equations {
    template<class RowsBegin, int TotalColumns, class Result>
    struct apply {
        typedef typename RowsBegin::item top_row;
        typedef typename determine_extra_equations_skip_zeros_impl<
            top_row::item::Numerator == 0,
            RowsBegin::item::size::value == 1
        >::template apply<
            RowsBegin,
            RowsBegin::size::value,
            TotalColumns - RemainingColumns,
            TotalColumns,
            Result
        > column_info;
        typedef typename determine_extra_equations<
            RemainingColumns - 1,
            column_info::next_equations::size::value == 0
        >::template apply<
            typename column_info::next_equations,
            TotalColumns,
            typename column_info::type
        >::type type;
    };
};

template<int RemainingColumns>
struct determine_extra_equations<RemainingColumns, true> {
    template<class RowsBegin, int TotalColumns, class Result>
    struct apply {
        typedef typename determine_extra_equations<RemainingColumns - 1, true>::template apply<
            RowsBegin,
            TotalColumns,
            list<typename create_row_of_identity<TotalColumns - RemainingColumns, TotalColumns>::type, Result>
        >::type type;
    };
};

template<>
struct determine_extra_equations<0, true> {
    template<class RowsBegin, int TotalColumns, class Result>
    struct apply {
        typedef Result type;
    };
};

// Stage 2
// invert the matrix using Gauss-Jordan elimination


template<bool is_zero, bool is_last>
struct invert_strip_leading_zeroes;

template<int N>
struct invert_handle_after_pivot_row;

// When processing column N, none of the first N rows 
// can be the pivot column.
template<int N>
struct invert_handle_inital_rows {
    template<class RowsBegin, class IdentityBegin>
    struct apply {
        typedef typename invert_handle_inital_rows<N - 1>::template apply<
            typename RowsBegin::next,
            typename IdentityBegin::next
        > next;
        typedef typename RowsBegin::item current_row;
        typedef typename IdentityBegin::item current_identity_row;
        typedef typename next::pivot_row pivot_row;
        typedef typename next::identity_pivot_row identity_pivot_row;
        typedef list<
            typename eliminate_from_pair_of_equations_impl<(current_row::size::value) - 1>::template apply<
                typename current_row::next,
                pivot_row,
                typename current_row::item,
                static_rational<1>
            >::type,
            typename next::new_matrix
        > new_matrix;
        typedef list<
            typename eliminate_from_pair_of_equations_impl<(current_identity_row::size::value)>::template apply<
                current_identity_row,
                identity_pivot_row,
                typename current_row::item,
                static_rational<1>
            >::type,
            typename next::identity_result
        > identity_result;
    };
};

// This handles the switch to searching for a pivot column.
// The pivot row will be propagated up in the typedefs
// pivot_row and identity_pivot_row.  It is inserted here.
template<>
struct invert_handle_inital_rows<0> {
    template<class RowsBegin, class IdentityBegin>
    struct apply {
        typedef typename RowsBegin::item current_row;
        typedef typename invert_strip_leading_zeroes<
            (current_row::item::Numerator == 0),
            (RowsBegin::size::value == 1)
        >::template apply<
            RowsBegin,
            IdentityBegin
        > next;
        // results
        typedef list<typename next::pivot_row, typename next::new_matrix> new_matrix;
        typedef list<typename next::identity_pivot_row, typename next::identity_result> identity_result;
        typedef typename next::pivot_row pivot_row;
        typedef typename next::identity_pivot_row identity_pivot_row;
    };
};

// The first internal element which is not zero.
template<>
struct invert_strip_leading_zeroes<false, false> {
    template<class RowsBegin, class IdentityBegin>