summaryrefslogtreecommitdiffstats
path: root/include/2geom/numeric/symmetric-matrix-fs.h
diff options
context:
space:
mode:
Diffstat (limited to '')
-rw-r--r--include/2geom/numeric/symmetric-matrix-fs.h733
1 files changed, 733 insertions, 0 deletions
diff --git a/include/2geom/numeric/symmetric-matrix-fs.h b/include/2geom/numeric/symmetric-matrix-fs.h
new file mode 100644
index 0000000..2fadd69
--- /dev/null
+++ b/include/2geom/numeric/symmetric-matrix-fs.h
@@ -0,0 +1,733 @@
+/*
+ * SymmetricMatrix, ConstSymmetricMatrixView, SymmetricMatrixView template
+ * classes implement fixed size symmetric matrix; "views" mimic the semantic
+ * of C++ references: any operation performed on a "view" is actually performed
+ * on the "viewed object"
+ *
+ * Authors:
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2009 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_SYMMETRIC_MATRIX_FS_H_
+#define _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+#include <2geom/utils.h>
+#include <2geom/exception.h>
+
+#include <boost/static_assert.hpp>
+
+#include <cassert>
+#include <utility> // for std::pair
+#include <algorithm> // for std::swap, std::copy
+#include <sstream>
+#include <string>
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+template <size_t I, size_t J, bool B = (I < J)>
+struct index
+{
+ static const size_t K = index<J, I>::K;
+};
+
+template <size_t I, size_t J>
+struct index<I, J, false>
+{
+ static const size_t K = (((I+1) * I) >> 1) + J;
+};
+
+} // end namespace detail
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix;
+
+template <size_t N>
+class BaseSymmetricMatrix;
+
+template <size_t N>
+class SymmetricMatrix;
+
+template <size_t N>
+class ConstSymmetricMatrixView;
+
+template <size_t N>
+class SymmetricMatrixView;
+
+
+
+// declaration needed for friend clause
+template <size_t N>
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+ ConstBaseSymmetricMatrix<N> const& _smatrix2);
+
+
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix
+{
+ public:
+ const static size_t DIM = N;
+ const static size_t DATA_SIZE = ((DIM+1) * DIM) / 2;
+
+ public:
+
+ ConstBaseSymmetricMatrix (VectorView const& _data)
+ : m_data(_data)
+ {
+ }
+
+ double operator() (size_t i, size_t j) const
+ {
+ return m_data[get_index(i,j)];
+ }
+
+ template <size_t I, size_t J>
+ double get() const
+ {
+ BOOST_STATIC_ASSERT ((I < N && J < N));
+ return m_data[detail::index<I, J>::K];
+ }
+
+
+ size_t rows() const
+ {
+ return DIM;
+ }
+
+ size_t columns() const
+ {
+ return DIM;
+ }
+
+ bool is_zero() const
+ {
+ return m_data.is_zero();
+ }
+
+ bool is_positive() const
+ {
+ return m_data.is_positive();
+ }
+
+ bool is_negative() const
+ {
+ return m_data.is_negative();
+ }
+
+ bool is_non_negative() const
+ {
+ return m_data.is_non_negative();
+ }
+
+ double min() const
+ {
+ return m_data.min();
+ }
+
+ double max() const
+ {
+ return m_data.max();
+ }
+
+ std::pair<size_t, size_t>
+ min_index() const
+ {
+ std::pair<size_t, size_t> indices(0,0);
+ double min_value = m_data[0];
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ for (size_t j = 0; j <= i; ++j)
+ {
+ if (min_value > (*this)(i,j))
+ {
+ min_value = (*this)(i,j);
+ indices.first = i;
+ indices.second = j;
+ }
+ }
+ }
+ return indices;
+ }
+
+ std::pair<size_t, size_t>
+ max_index() const
+ {
+ std::pair<size_t, size_t> indices(0,0);
+ double max_value = m_data[0];
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ for (size_t j = 0; j <= i; ++j)
+ {
+ if (max_value < (*this)(i,j))
+ {
+ max_value = (*this)(i,j);
+ indices.first = i;
+ indices.second = j;
+ }
+ }
+ }
+ return indices;
+ }
+
+ size_t min_on_row_index (size_t i) const
+ {
+ size_t idx = 0;
+ double min_value = (*this)(i,0);
+ for (size_t j = 1; j < DIM; ++j)
+ {
+ if (min_value > (*this)(i,j))
+ {
+ min_value = (*this)(i,j);
+ idx = j;
+ }
+ }
+ return idx;
+ }
+
+ size_t max_on_row_index (size_t i) const
+ {
+ size_t idx = 0;
+ double max_value = (*this)(i,0);
+ for (size_t j = 1; j < DIM; ++j)
+ {
+ if (max_value < (*this)(i,j))
+ {
+ max_value = (*this)(i,j);
+ idx = j;
+ }
+ }
+ return idx;
+ }
+
+ size_t min_on_column_index (size_t j) const
+ {
+ return min_on_row_index(j);
+ }
+
+ size_t max_on_column_index (size_t j) const
+ {
+ return max_on_row_index(j);
+ }
+
+ size_t min_on_diag_index () const
+ {
+ size_t idx = 0;
+ double min_value = (*this)(0,0);
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ if (min_value > (*this)(i,i))
+ {
+ min_value = (*this)(i,i);
+ idx = i;
+ }
+ }
+ return idx;
+ }
+
+ size_t max_on_diag_index () const
+ {
+ size_t idx = 0;
+ double max_value = (*this)(0,0);
+ for (size_t i = 1; i < DIM; ++i)
+ {
+ if (max_value < (*this)(i,i))
+ {
+ max_value = (*this)(i,i);
+ idx = i;
+ }
+ }
+ return idx;
+ }
+
+ std::string str() const;
+
+ ConstSymmetricMatrixView<N-1> main_minor_const_view() const;
+
+ SymmetricMatrix<N> operator- () const;
+
+ Vector operator* (ConstVectorView _vector) const
+ {
+ assert (_vector.size() == DIM);
+ Vector result(DIM, 0.0);
+
+ for (size_t i = 0; i < DIM; ++i)
+ {
+ for (size_t j = 0; j < DIM; ++j)
+ {
+ result[i] += (*this)(i,j) * _vector[j];
+ }
+ }
+ return result;
+ }
+
+ protected:
+ static size_t get_index (size_t i, size_t j)
+ {
+ if (i < j) return get_index (j, i);
+ size_t k = (i+1) * i;
+ k >>= 1;
+ k += j;
+ return k;
+ }
+
+ protected:
+ ConstVectorView get_data() const
+ {
+ return m_data;
+ }
+
+ friend
+ bool operator==<N> (ConstBaseSymmetricMatrix const& _smatrix1,
+ ConstBaseSymmetricMatrix const& _smatrix2);
+
+ protected:
+ VectorView m_data;
+
+}; //end ConstBaseSymmetricMatrix
+
+
+template <size_t N>
+class BaseSymmetricMatrix : public ConstBaseSymmetricMatrix<N>
+{
+ public:
+ typedef ConstBaseSymmetricMatrix<N> base_type;
+
+
+ public:
+
+ BaseSymmetricMatrix (VectorView const& _data)
+ : base_type(_data)
+ {
+ }
+
+ double operator() (size_t i, size_t j) const
+ {
+ return m_data[base_type::get_index(i,j)];
+ }
+
+ double& operator() (size_t i, size_t j)
+ {
+ return m_data[base_type::get_index(i,j)];
+ }
+
+ template <size_t I, size_t J>
+ double& get()
+ {
+ BOOST_STATIC_ASSERT ((I < N && J < N));
+ return m_data[detail::index<I, J>::K];
+ }
+
+ void set_all (double x)
+ {
+ m_data.set_all(x);
+ }
+
+ SymmetricMatrixView<N-1> main_minor_view();
+
+ BaseSymmetricMatrix& transpose() const
+ {
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& translate (double c)
+ {
+ m_data.translate(c);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& scale (double c)
+ {
+ m_data.scale(c);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& operator+= (base_type const& _smatrix)
+ {
+ m_data += (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+ return (*this);
+ }
+
+ BaseSymmetricMatrix& operator-= (base_type const& _smatrix)
+ {
+ m_data -= (static_cast<const BaseSymmetricMatrix &>(_smatrix).m_data);
+ return (*this);
+ }
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+ using base_type::operator-;
+ using base_type::operator*;
+
+}; //end BaseSymmetricMatrix
+
+
+template <size_t N>
+class SymmetricMatrix : public BaseSymmetricMatrix<N>
+{
+ public:
+ typedef BaseSymmetricMatrix<N> base_type;
+ typedef typename base_type::base_type base_base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+ public:
+ SymmetricMatrix ()
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ }
+
+ explicit
+ SymmetricMatrix (ConstVectorView _data)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ assert (_data.size() == DATA_SIZE);
+ m_data = _data;
+ }
+
+ explicit
+ SymmetricMatrix (const double _data[DATA_SIZE])
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ std::copy (_data, _data + DATA_SIZE, m_adata);
+ }
+
+ SymmetricMatrix (SymmetricMatrix const& _smatrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ m_data = _smatrix.m_data;
+ }
+
+ explicit
+ SymmetricMatrix (base_base_type const& _smatrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ }
+
+ explicit
+ SymmetricMatrix (ConstMatrixView const& _matrix)
+ : base_type (VectorView(m_adata, DATA_SIZE))
+ {
+ assert (_matrix.rows() == N && _matrix.columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j <= i ; ++j)
+ (*this)(i,j) = _matrix(i,j);
+ }
+
+ SymmetricMatrix& operator= (SymmetricMatrix const& _smatrix)
+ {
+ m_data = _smatrix.m_data;
+ return (*this);
+ }
+
+ SymmetricMatrix& operator= (base_base_type const& _smatrix)
+ {
+
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ return (*this);
+ }
+
+ SymmetricMatrix& operator= (ConstMatrixView const& _matrix)
+ {
+ assert (_matrix.rows() == N && _matrix.columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j <= i ; ++j)
+ (*this)(i,j) = _matrix(i,j);
+
+ return (*this);
+ }
+
+ // needed for accessing m_adata
+ friend class ConstSymmetricMatrixView<DIM>;
+ friend class SymmetricMatrixView<DIM>;
+ private:
+ double m_adata[DATA_SIZE];
+}; //end SymmetricMatrix
+
+
+template <size_t N>
+class ConstSymmetricMatrixView : public ConstBaseSymmetricMatrix<N>
+{
+ public:
+ typedef ConstBaseSymmetricMatrix<N> base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+
+ public:
+
+ explicit
+ ConstSymmetricMatrixView (ConstVectorView _data)
+ : base_type (const_vector_view_cast(_data))
+ {
+ assert (_data.size() == DATA_SIZE);
+ }
+
+ explicit
+ ConstSymmetricMatrixView (const double _data[DATA_SIZE])
+ : base_type (const_vector_view_cast (ConstVectorView (_data, DATA_SIZE)))
+ {
+ }
+
+ ConstSymmetricMatrixView (const ConstSymmetricMatrixView & _smatrix)
+ : base_type (_smatrix.m_data)
+ {
+ }
+
+ ConstSymmetricMatrixView (const base_type & _smatrix)
+ : base_type (static_cast<const ConstSymmetricMatrixView &>(_smatrix).m_data)
+ {
+ }
+
+}; //end SymmetricMatrix
+
+
+// declaration needed for friend clause
+template <size_t N>
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2);
+
+
+template <size_t N>
+class SymmetricMatrixView : public BaseSymmetricMatrix<N>
+{
+ public:
+ typedef BaseSymmetricMatrix<N> base_type;
+ typedef typename base_type::base_type base_base_type;
+
+ using base_type::DIM;
+ using base_type::DATA_SIZE;
+ using base_type::m_data;
+
+ public:
+
+ explicit
+ SymmetricMatrixView (VectorView _data)
+ : base_type (_data)
+ {
+ assert (_data.size() == DATA_SIZE);
+ }
+
+ explicit
+ SymmetricMatrixView (double _data[DATA_SIZE])
+ : base_type (VectorView (_data, DATA_SIZE))
+ {
+ }
+
+ SymmetricMatrixView (const SymmetricMatrixView & _smatrix)
+ : base_type (_smatrix.m_data)
+ {
+ }
+
+ SymmetricMatrixView (SymmetricMatrix<DIM> & _smatrix)
+ : base_type (VectorView (_smatrix.m_adata, DATA_SIZE))
+ {
+ }
+
+ SymmetricMatrixView& operator= (const SymmetricMatrixView & _smatrix)
+ {
+ m_data = _smatrix.m_data;
+ return (*this);
+ }
+
+ SymmetricMatrixView& operator= (const base_base_type & _smatrix)
+ {
+ m_data = static_cast<const ConstSymmetricMatrixView<DIM> &>(_smatrix).m_data;
+ return (*this);
+ }
+
+ friend
+ void swap_view<N>(SymmetricMatrixView & m1, SymmetricMatrixView & m2);
+
+}; //end SymmetricMatrix
+
+
+
+
+/*
+ * class ConstBaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+std::string ConstBaseSymmetricMatrix<N>::str() const
+{
+ std::ostringstream oss;
+ oss << (*this);
+ return oss.str();
+}
+
+template <size_t N>
+inline
+ConstSymmetricMatrixView<N-1>
+ConstBaseSymmetricMatrix<N>::main_minor_const_view() const
+{
+ ConstVectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+ ConstSymmetricMatrixView<N-1> mm(data);
+ return mm;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> ConstBaseSymmetricMatrix<N>::operator- () const
+{
+ SymmetricMatrix<N> result;
+ for (size_t i = 0; i < DATA_SIZE; ++i)
+ {
+ result.m_data[i] = -m_data[i];
+ }
+ return result;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix friend free functions
+ */
+
+template <size_t N>
+inline
+bool operator== (ConstBaseSymmetricMatrix<N> const& _smatrix1,
+ ConstBaseSymmetricMatrix<N> const& _smatrix2)
+{
+ return (_smatrix1.m_data == _smatrix2.m_data);
+}
+
+/*
+ * class ConstBaseSymmetricMatrix related free functions
+ */
+
+template< size_t N, class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os,
+ const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+ os << "[[" << _matrix(0,0);
+ for (size_t j = 1; j < N; ++j)
+ {
+ os << ", " << _matrix(0,j);
+ }
+ os << "]";
+ for (size_t i = 1; i < N; ++i)
+ {
+ os << "\n [" << _matrix(i,0);
+ for (size_t j = 1; j < N; ++j)
+ {
+ os << ", " << _matrix(i,j);
+ }
+ os << "]";
+ }
+ os << "]";
+ return os;
+}
+
+
+/*
+ * class ConstBaseSymmetricMatrix specialized methods
+ */
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<2>::get_index (size_t i, size_t j)
+{
+ return (i+j);
+}
+
+template<>
+inline
+size_t ConstBaseSymmetricMatrix<3>::get_index (size_t i, size_t j)
+{
+ size_t k = i + j;
+ if (i == 2 || j == 2) ++k;
+ return k;
+}
+
+
+/*
+ * class BaseSymmetricMatrix methods
+ */
+
+template <size_t N>
+inline
+SymmetricMatrixView<N-1> BaseSymmetricMatrix<N>::main_minor_view()
+{
+ VectorView data(m_data.get_gsl_vector()->data, DATA_SIZE - DIM);
+ SymmetricMatrixView<N-1> mm(data);
+ return mm;
+}
+
+
+/*
+ * class SymmetricMatrixView friend free functions
+ */
+
+template <size_t N>
+inline
+void swap_view(SymmetricMatrixView<N> & m1, SymmetricMatrixView<N> & m2)
+{
+ swap_view(m1.m_data, m2.m_data);
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif // _NL_SYMMETRIC_MATRIX_FS_H_
+
+
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-24 07:49:49 +0000