hermite.cpp

Go to the documentation of this file.

/* /////////////////////////////////////////////////////////////////////////////
//
//  File:       hermite.cpp
//  Created:    2005/10/18    23:24
//
//  Author: Andrey Somsikov
*/

#include "config.hpp"

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

#include "hermite.hpp"
#include <iostream>

namespace Arageli
{

template <class T>
matrix<T> hermite_old (const matrix<T> &m)
{
    typedef matrix<T> M;
    typedef typename M::value_type MV;
    typedef typename M::size_type size_type;

    M res(m);

    size_type i = 0, j = 0;
    while(i < res.nrows()-1 && j < res.ncols()-1)
    {
/*
output_aligned(std::cout, res, "|| ", " ||", "  ");
std::cout << "\n";
*/
        // ensure that el(i,j)>0
        if (is_null(res.el(i, j)))
        {
            size_type k;
            // find non zero element in column j
            for (k = i+1; k < res.nrows(); k++)
                if (!is_null(res.el(k, j)))
                    break;
            if (k == res.nrows())
            {
                // skip column if not found
                j++;
                continue;
            }
            res.swap_rows(i, k);
/*
            // ???
            if (res.el(i, j) < factory<MV>::null())
                res.mult_row(i, factory<MV>::opposite_unit());
*/
        }

        // find pivot element
        size_type minValueIndex = i;
        bool isNonZero = false;
        for (size_type k = i+1; k < res.ncols(); k++)
        {
            if (!is_null(res.el(k, j)))
            {   
                isNonZero = true;
                if (res.el(k, j).degree() < res.el(minValueIndex, j).degree())
                    minValueIndex = k;
            }
        }
        if (!isNonZero)
        {
            // nothing to do with this element
            i++;
            j++;
            continue;
        }
        res.swap_rows(i, minValueIndex);
/*
        if (res.el(i, j) < factory<MV>::null())
            res.mult_row(i, factory<MV>::opposite_unit());
*/

        for (size_type k = 0; k < i; k++)
            res.addmult_rows(k, i, factory<MV>::opposite_unit()*(res.el(k, j) / res.el(i, j)));
        for (size_type k = i+1; k < res.ncols(); k++)
            res.addmult_rows(k, i, factory<MV>::opposite_unit()*(res.el(k, j) / res.el(i, j)));
    }

    return res;
}

} // namespace Arageli

#endif // #ifndef ARAGELI_INCLUDE_CPP_WITH_EXPORT_TEMPLATE

/* End of file hermite.cpp */

---