summaryrefslogtreecommitdiffstats
path: root/include/2geom/numeric/matrix.h
diff options
context:
space:
mode:
Diffstat (limited to 'include/2geom/numeric/matrix.h')
-rw-r--r--include/2geom/numeric/matrix.h603
1 files changed, 603 insertions, 0 deletions
diff --git a/include/2geom/numeric/matrix.h b/include/2geom/numeric/matrix.h
new file mode 100644
index 0000000..02851b4
--- /dev/null
+++ b/include/2geom/numeric/matrix.h
@@ -0,0 +1,603 @@
+/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "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 2008 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_MATRIX_H_
+#define _NL_MATRIX_H_
+
+#include <2geom/exception.h>
+#include <2geom/numeric/vector.h>
+
+#include <cassert>
+#include <utility> // for std::pair
+#include <algorithm> // for std::swap
+#include <sstream>
+#include <string>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseMatrixImpl
+{
+ public:
+ virtual ~BaseMatrixImpl()
+ {
+ }
+
+ ConstVectorView row_const_view(size_t i) const
+ {
+ return ConstVectorView(gsl_matrix_const_row(m_matrix, i));
+ }
+
+ ConstVectorView column_const_view(size_t i) const
+ {
+ return ConstVectorView(gsl_matrix_const_column(m_matrix, i));
+ }
+
+ const double & operator() (size_t i, size_t j) const
+ {
+ return *gsl_matrix_const_ptr(m_matrix, i, j);
+ }
+
+ const gsl_matrix* get_gsl_matrix() const
+ {
+ return m_matrix;
+ }
+
+ bool is_zero() const
+ {
+ return gsl_matrix_isnull(m_matrix);
+ }
+
+ bool is_positive() const
+ {
+ for ( unsigned int i = 0; i < rows(); ++i )
+ {
+ for ( unsigned int j = 0; j < columns(); ++j )
+ {
+ if ( (*this)(i,j) <= 0 ) return false;
+ }
+ }
+ return true;
+ }
+
+ bool is_negative() const
+ {
+ for ( unsigned int i = 0; i < rows(); ++i )
+ {
+ for ( unsigned int j = 0; j < columns(); ++j )
+ {
+ if ( (*this)(i,j) >= 0 ) return false;
+ }
+ }
+ return true;
+ }
+
+ bool is_non_negative() const
+ {
+ for ( unsigned int i = 0; i < rows(); ++i )
+ {
+ for ( unsigned int j = 0; j < columns(); ++j )
+ {
+ if ( (*this)(i,j) < 0 ) return false;
+ }
+ }
+ return true;
+ }
+
+ double max() const
+ {
+ return gsl_matrix_max(m_matrix);
+ }
+
+ double min() const
+ {
+ return gsl_matrix_min(m_matrix);
+ }
+
+ std::pair<size_t, size_t>
+ max_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ std::pair<size_t, size_t>
+ min_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ size_t rows() const
+ {
+ return m_rows;
+ }
+
+ size_t columns() const
+ {
+ return m_columns;
+ }
+
+ std::string str() const;
+
+ protected:
+ size_t m_rows, m_columns;
+ gsl_matrix* m_matrix;
+
+}; // end class BaseMatrixImpl
+
+
+inline
+bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)
+{
+ if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;
+
+ for (size_t i = 0; i < m1.rows(); ++i)
+ for (size_t j = 0; j < m1.columns(); ++j)
+ if (m1(i,j) != m2(i,j)) return false;
+
+ return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)
+{
+ if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;
+
+ os << "[[" << _matrix(0,0);
+ for (size_t j = 1; j < _matrix.columns(); ++j)
+ {
+ os << ", " << _matrix(0,j);
+ }
+ os << "]";
+
+ for (size_t i = 1; i < _matrix.rows(); ++i)
+ {
+ os << ", [" << _matrix(i,0);
+ for (size_t j = 1; j < _matrix.columns(); ++j)
+ {
+ os << ", " << _matrix(i,j);
+ }
+ os << "]";
+ }
+ os << "]";
+ return os;
+}
+
+inline
+std::string BaseMatrixImpl::str() const
+{
+ std::ostringstream oss;
+ oss << (*this);
+ return oss.str();
+}
+
+
+class MatrixImpl : public BaseMatrixImpl
+{
+ public:
+
+ typedef BaseMatrixImpl base_type;
+
+ void set_all( double x )
+ {
+ gsl_matrix_set_all(m_matrix, x);
+ }
+
+ void set_identity()
+ {
+ gsl_matrix_set_identity(m_matrix);
+ }
+
+ using base_type::operator(); // VSC legacy support
+ const double & operator() (size_t i, size_t j) const
+ {
+ return base_type::operator ()(i, j);
+ }
+
+ double & operator() (size_t i, size_t j)
+ {
+ return *gsl_matrix_ptr(m_matrix, i, j);
+ }
+
+ using base_type::get_gsl_matrix;
+
+ gsl_matrix* get_gsl_matrix()
+ {
+ return m_matrix;
+ }
+
+ VectorView row_view(size_t i)
+ {
+ return VectorView(gsl_matrix_row(m_matrix, i));
+ }
+
+ VectorView column_view(size_t i)
+ {
+ return VectorView(gsl_matrix_column(m_matrix, i));
+ }
+
+ void swap_rows(size_t i, size_t j)
+ {
+ gsl_matrix_swap_rows(m_matrix, i, j);
+ }
+
+ void swap_columns(size_t i, size_t j)
+ {
+ gsl_matrix_swap_columns(m_matrix, i, j);
+ }
+
+ MatrixImpl & transpose()
+ {
+ assert(columns() == rows());
+ gsl_matrix_transpose(m_matrix);
+ return (*this);
+ }
+
+ MatrixImpl & scale(double x)
+ {
+ gsl_matrix_scale(m_matrix, x);
+ return (*this);
+ }
+
+ MatrixImpl & translate(double x)
+ {
+ gsl_matrix_add_constant(m_matrix, x);
+ return (*this);
+ }
+
+ MatrixImpl & operator+=(base_type const& _matrix)
+ {
+ gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());
+ return (*this);
+ }
+
+ MatrixImpl & operator-=(base_type const& _matrix)
+ {
+ gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());
+ return (*this);
+ }
+
+}; // end class MatrixImpl
+
+} // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+
+template <size_t N>
+class ConstBaseSymmetricMatrix;
+
+
+class Matrix: public detail::MatrixImpl
+{
+ public:
+ typedef detail::MatrixImpl base_type;
+
+ public:
+ // the matrix is not initialized
+ Matrix(size_t n1, size_t n2)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = gsl_matrix_alloc(n1, n2);
+ }
+
+ Matrix(size_t n1, size_t n2, double x)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = gsl_matrix_alloc(n1, n2);
+ gsl_matrix_set_all(m_matrix, x);
+ }
+
+ Matrix(Matrix const& _matrix)
+ : base_type()
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = gsl_matrix_alloc(rows(), columns());
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ }
+
+ explicit
+ Matrix(base_type::base_type const& _matrix)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = gsl_matrix_alloc(rows(), columns());
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ }
+
+ template <size_t N>
+ explicit
+ Matrix(ConstBaseSymmetricMatrix<N> const& _smatrix)
+ {
+ m_rows = N;
+ m_columns = N;
+ m_matrix = gsl_matrix_alloc(N, N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j < N ; ++j)
+ (*gsl_matrix_ptr(m_matrix, i, j)) = _smatrix(i,j);
+ }
+
+ Matrix & operator=(Matrix const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ Matrix & operator=(base_type::base_type const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ template <size_t N>
+ Matrix & operator=(ConstBaseSymmetricMatrix<N> const& _smatrix)
+ {
+ assert (rows() == N && columns() == N);
+ for (size_t i = 0; i < N; ++i)
+ for (size_t j = 0; j < N ; ++j)
+ (*this)(i,j) = _smatrix(i,j);
+ return *this;
+ }
+
+ ~Matrix() override
+ {
+ gsl_matrix_free(m_matrix);
+ }
+
+ Matrix & transpose()
+ {
+ return static_cast<Matrix &>( base_type::transpose() );
+ }
+
+ Matrix & scale(double x)
+ {
+ return static_cast<Matrix &>( base_type::scale(x) );
+ }
+
+ Matrix & translate(double x)
+ {
+ return static_cast<Matrix &>( base_type::translate(x) );
+ }
+
+ Matrix & operator+=(base_type::base_type const& _matrix)
+ {
+ return static_cast<Matrix &>( base_type::operator+=(_matrix) );
+ }
+
+ Matrix & operator-=(base_type::base_type const& _matrix)
+ {
+ return static_cast<Matrix &>( base_type::operator-=(_matrix) );
+ }
+
+ friend
+ void swap(Matrix & m1, Matrix & m2);
+ friend
+ void swap_any(Matrix & m1, Matrix & m2);
+
+}; // end class Matrix
+
+
+// warning! this operation invalidates any view of the passed matrix objects
+inline
+void swap(Matrix & m1, Matrix & m2)
+{
+ assert(m1.rows() == m2.rows() && m1.columns() == m2.columns());
+ using std::swap;
+ swap(m1.m_matrix, m2.m_matrix);
+}
+
+inline void swap_any(Matrix &m1, Matrix &m2)
+{
+ using std::swap;
+ swap(m1.m_matrix, m2.m_matrix);
+ swap(m1.m_rows, m2.m_rows);
+ swap(m1.m_columns, m2.m_columns);
+}
+
+
+
+class ConstMatrixView : public detail::BaseMatrixImpl
+{
+ public:
+ typedef detail::BaseMatrixImpl base_type;
+
+ public:
+ ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+ : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ ConstMatrixView(const ConstMatrixView & _matrix)
+ : base_type(),
+ m_matrix_view(_matrix.m_matrix_view)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ ConstMatrixView(const base_type & _matrix)
+ : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+ }
+
+ private:
+ gsl_matrix_const_view m_matrix_view;
+
+}; // end class ConstMatrixView
+
+
+
+
+class MatrixView : public detail::MatrixImpl
+{
+ public:
+ typedef detail::MatrixImpl base_type;
+
+ public:
+ MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+ {
+ m_rows = n1;
+ m_columns = n2;
+ m_matrix_view
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView(const MatrixView & _matrix)
+ : base_type()
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix_view = _matrix.m_matrix_view;
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView(Matrix & _matrix)
+ {
+ m_rows = _matrix.rows();
+ m_columns = _matrix.columns();
+ m_matrix_view
+ = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());
+ m_matrix = &(m_matrix_view.matrix);
+ }
+
+ MatrixView & operator=(MatrixView const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
+ return *this;
+ }
+
+ MatrixView & operator=(base_type::base_type const& _matrix)
+ {
+ assert( rows() == _matrix.rows() && columns() == _matrix.columns() );
+ gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+ return *this;
+ }
+
+ MatrixView & transpose()
+ {
+ return static_cast<MatrixView &>( base_type::transpose() );
+ }
+
+ MatrixView & scale(double x)
+ {
+ return static_cast<MatrixView &>( base_type::scale(x) );
+ }
+
+ MatrixView & translate(double x)
+ {
+ return static_cast<MatrixView &>( base_type::translate(x) );
+ }
+
+ MatrixView & operator+=(base_type::base_type const& _matrix)
+ {
+ return static_cast<MatrixView &>( base_type::operator+=(_matrix) );
+ }
+
+ MatrixView & operator-=(base_type::base_type const& _matrix)
+ {
+ return static_cast<MatrixView &>( base_type::operator-=(_matrix) );
+ }
+
+ friend
+ void swap_view(MatrixView & m1, MatrixView & m2);
+
+ private:
+ gsl_matrix_view m_matrix_view;
+
+}; // end class MatrixView
+
+
+inline
+void swap_view(MatrixView & m1, MatrixView & m2)
+{
+ assert(m1.rows() == m2.rows() && m1.columns() == m2.columns());
+ using std::swap;
+ swap(m1.m_matrix_view, m2.m_matrix_view);
+}
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseVectorImpl const& v );
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+ detail::BaseMatrixImpl const& B );
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
+
+double trace (detail::BaseMatrixImpl const& A);
+
+double det (detail::BaseMatrixImpl const& A);
+
+} } // end namespaces
+
+#endif /*_NL_MATRIX_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-21 11:00:23 +0000