From 3cd345ae277f34e13420e4f7849f4e030b2437d6 Mon Sep 17 00:00:00 2001 From: mcecchetti Date: Tue, 20 May 2008 22:29:23 +0000 Subject: synchronization with 2geom library (bzr r5723) --- src/2geom/numeric/linear_system.h | 178 ++++++++-------- src/2geom/numeric/matrix.h | 414 +++++++++++++++++++------------------- src/2geom/numeric/vector.h | 369 ++++++++++++++++----------------- 3 files changed, 481 insertions(+), 480 deletions(-) (limited to 'src/2geom/numeric') diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h index 2b4963998..9cb521eb2 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,89 +1,89 @@ -#ifndef _NL_LINEAR_SYSTEM_H_ -#define _NL_LINEAR_SYSTEM_H_ - - -#include - -#include - -#include "numeric/matrix.h" -#include "numeric/vector.h" - - -namespace Geom { namespace NL { - - -class LinearSystem -{ -public: - LinearSystem(Matrix & _matrix, Vector & _vector) - : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) - { - } - - const Vector & LU_solve() - { - assert( matrix().rows() == matrix().columns() - && matrix().rows() == vector().size() ); - int s; - gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); - gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); - gsl_linalg_LU_solve( matrix().get_gsl_matrix(), - p, - vector().get_gsl_vector(), - m_solution.get_gsl_vector() - ); - gsl_permutation_free(p); - return solution(); - } - - const Vector & SV_solve() - { - assert( matrix().rows() >= matrix().columns() - && matrix().rows() == vector().size() ); - - gsl_matrix* U = matrix().get_gsl_matrix(); - gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); - gsl_vector* S = gsl_vector_alloc(matrix().columns()); - gsl_vector* work = gsl_vector_alloc(matrix().columns()); - - gsl_linalg_SV_decomp( U, V, S, work ); - - gsl_vector* b = vector().get_gsl_vector(); - gsl_vector* x = m_solution.get_gsl_vector(); - - gsl_linalg_SV_solve( U, V, S, b, x); - - gsl_matrix_free(V); - gsl_vector_free(S); - gsl_vector_free(work); - - return solution(); - } - - Matrix & matrix() - { - return m_matrix; - } - - Vector & vector() - { - return m_vector; - } - - const Vector & solution() const - { - return m_solution; - } - -private: - Matrix & m_matrix; - Vector & m_vector; - Vector m_solution; -}; - - -} } // end namespaces - - -#endif /*_NL_LINEAR_SYSTEM_H_*/ +#ifndef _NL_LINEAR_SYSTEM_H_ +#define _NL_LINEAR_SYSTEM_H_ + + +#include + +#include + +#include "numeric/matrix.h" +#include "numeric/vector.h" + + +namespace Geom { namespace NL { + + +class LinearSystem +{ +public: + LinearSystem(Matrix & _matrix, Vector & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + const Vector & LU_solve() + { + assert( matrix().rows() == matrix().columns() + && matrix().rows() == vector().size() ); + int s; + gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); + gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); + gsl_linalg_LU_solve( matrix().get_gsl_matrix(), + p, + vector().get_gsl_vector(), + m_solution.get_gsl_vector() + ); + gsl_permutation_free(p); + return solution(); + } + + const Vector & SV_solve() + { + assert( matrix().rows() >= matrix().columns() + && matrix().rows() == vector().size() ); + + gsl_matrix* U = matrix().get_gsl_matrix(); + gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); + gsl_vector* S = gsl_vector_alloc(matrix().columns()); + gsl_vector* work = gsl_vector_alloc(matrix().columns()); + + gsl_linalg_SV_decomp( U, V, S, work ); + + gsl_vector* b = vector().get_gsl_vector(); + gsl_vector* x = m_solution.get_gsl_vector(); + + gsl_linalg_SV_solve( U, V, S, b, x); + + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + + return solution(); + } + + Matrix & matrix() + { + return m_matrix; + } + + Vector & vector() + { + return m_vector; + } + + const Vector & solution() const + { + return m_solution; + } + +private: + Matrix & m_matrix; + Vector & m_vector; + Vector m_solution; +}; + + +} } // end namespaces + + +#endif /*_NL_LINEAR_SYSTEM_H_*/ diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h index bdfab75ef..bdd647e53 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,207 +1,207 @@ -#ifndef _NL_MATRIX_H_ -#define _NL_MATRIX_H_ - - -#include -#include -#include - -#include - - - -namespace Geom { namespace NL { - -class Matrix; -void swap(Matrix & m1, Matrix & m2); - - -class Matrix -{ -public: - // the matrix is not inizialized - 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 ) - : m_rows(_matrix.rows()), m_columns(_matrix.columns()) - { - m_matrix = gsl_matrix_alloc(rows(), columns()); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - } - -// explicit -// Matrix( gsl_matrix* m, size_t n1, size_t n2) -// : m_rows(n1), m_columns(n2), m_matrix(m) -// { -// } - - Matrix & operator=(Matrix const& _matrix) - { - assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - - return *this; - } - - virtual ~Matrix() - { - gsl_matrix_free(m_matrix); - } - - void set_all( double x = 0) - { - gsl_matrix_set_all(m_matrix, x); - } - - void set_identity() - { - gsl_matrix_set_identity(m_matrix); - } - - const double & operator() (size_t i, size_t j) const - { - return *gsl_matrix_const_ptr(m_matrix, i, j); - } - - double & operator() (size_t i, size_t j) - { - return *gsl_matrix_ptr(m_matrix, i, j); - } - - gsl_matrix* get_gsl_matrix() - { - return m_matrix; - } - - 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); - } - - Matrix & transpose() - { - gsl_matrix_transpose(m_matrix); - return (*this); - } - - Matrix & scale(double x) - { - gsl_matrix_scale(m_matrix, x); - return (*this); - } - - Matrix & translate(double x) - { - gsl_matrix_add_constant(m_matrix, x); - return (*this); - } - - Matrix & operator+=(Matrix const& _matrix) - { - gsl_matrix_add(m_matrix, _matrix.m_matrix); - return (*this); - } - - Matrix & operator-=(Matrix const& _matrix) - { - gsl_matrix_sub(m_matrix, _matrix.m_matrix); - return (*this); - } - - bool is_zero() const - { - return gsl_matrix_isnull(m_matrix); - } - - bool is_positive() const - { - return gsl_matrix_ispos(m_matrix); - } - - bool is_negative() const - { - return gsl_matrix_isneg(m_matrix); - } - - 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 - max_index() const - { - std::pair indices; - gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - std::pair - min_index() const - { - std::pair indices; - gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - friend - void swap(Matrix & m1, Matrix & m2); - - size_t rows() const - { - return m_rows; - } - - size_t columns() const - { - return m_columns; - } - -private: - size_t m_rows, m_columns; - gsl_matrix* m_matrix; -}; - -void swap(Matrix & m1, Matrix & m2) -{ - assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); - std::swap(m1.m_matrix, m2.m_matrix); -} - - -} } // end namespaces - -#endif /*_NL_MATRIX_H_*/ +#ifndef _NL_MATRIX_H_ +#define _NL_MATRIX_H_ + + +#include +#include +#include + +#include + + + +namespace Geom { namespace NL { + +class Matrix; +void swap(Matrix & m1, Matrix & m2); + + +class Matrix +{ +public: + // the matrix is not inizialized + 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 ) + : m_rows(_matrix.rows()), m_columns(_matrix.columns()) + { + m_matrix = gsl_matrix_alloc(rows(), columns()); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + } + +// explicit +// Matrix( gsl_matrix* m, size_t n1, size_t n2) +// : m_rows(n1), m_columns(n2), m_matrix(m) +// { +// } + + Matrix & operator=(Matrix const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + + return *this; + } + + virtual ~Matrix() + { + gsl_matrix_free(m_matrix); + } + + void set_all( double x = 0) + { + gsl_matrix_set_all(m_matrix, x); + } + + void set_identity() + { + gsl_matrix_set_identity(m_matrix); + } + + const double & operator() (size_t i, size_t j) const + { + return *gsl_matrix_const_ptr(m_matrix, i, j); + } + + double & operator() (size_t i, size_t j) + { + return *gsl_matrix_ptr(m_matrix, i, j); + } + + gsl_matrix* get_gsl_matrix() + { + return m_matrix; + } + + 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); + } + + Matrix & transpose() + { + gsl_matrix_transpose(m_matrix); + return (*this); + } + + Matrix & scale(double x) + { + gsl_matrix_scale(m_matrix, x); + return (*this); + } + + Matrix & translate(double x) + { + gsl_matrix_add_constant(m_matrix, x); + return (*this); + } + + Matrix & operator+=(Matrix const& _matrix) + { + gsl_matrix_add(m_matrix, _matrix.m_matrix); + return (*this); + } + + Matrix & operator-=(Matrix const& _matrix) + { + gsl_matrix_sub(m_matrix, _matrix.m_matrix); + return (*this); + } + + bool is_zero() const + { + return gsl_matrix_isnull(m_matrix); + } + + bool is_positive() const + { + return gsl_matrix_ispos(m_matrix); + } + + bool is_negative() const + { + return gsl_matrix_isneg(m_matrix); + } + + 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 + max_index() const + { + std::pair indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair + min_index() const + { + std::pair indices; + gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + friend + void swap(Matrix & m1, Matrix & m2); + + size_t rows() const + { + return m_rows; + } + + size_t columns() const + { + return m_columns; + } + +private: + size_t m_rows, m_columns; + gsl_matrix* m_matrix; +}; + +void swap(Matrix & m1, Matrix & m2) +{ + assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); + std::swap(m1.m_matrix, m2.m_matrix); +} + + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h index 136b2cc3f..b25861e22 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,184 +1,185 @@ -#ifndef _NL_VECTOR_H_ -#define _NL_VECTOR_H_ - -#include -#include - -#include - - -namespace Geom { namespace NL { - -class Vector; -void swap(Vector & v1, Vector & v2); - - -class Vector -{ -public: - Vector(size_t n) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - } - - Vector(size_t n, double x) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_all(m_vector, x); - } - - // create a vector with n elements all set to zero - // but the i-th that is set to 1 - Vector(size_t n, size_t i) - : m_size(n) - { - m_vector = gsl_vector_alloc(n); - gsl_vector_set_basis(m_vector, i); - } - - Vector(Vector const& _vector) - : m_size(_vector.size()) - { - m_vector = gsl_vector_alloc(size()); - gsl_vector_memcpy(m_vector, _vector.m_vector); - } - - virtual ~Vector() - { - gsl_vector_free(m_vector); - } - - Vector & operator=(Vector const& _vector) - { - assert( size() == _vector.size() ); - gsl_vector_memcpy(m_vector, _vector.m_vector); - } - - void set_all(double x = 0) - { - gsl_vector_set_all(m_vector, x); - } - - void set_basis(size_t i) - { - gsl_vector_set_basis(m_vector, i); - } - - double const& operator[](size_t i) const - { - return *gsl_vector_const_ptr(m_vector, i); - } - - double & operator[](size_t i) - { - return *gsl_vector_ptr(m_vector, i); - } - - gsl_vector* get_gsl_vector() - { - return m_vector; - } - - void swap_elements(size_t i, size_t j) - { - gsl_vector_swap_elements(m_vector, i, j); - } - - void reverse() - { - gsl_vector_reverse(m_vector); - } - - Vector & scale(double x) - { - gsl_vector_scale(m_vector, x); - return (*this); - } - - Vector & translate(double x) - { - gsl_vector_add_constant(m_vector, x); - return (*this); - } - - Vector & operator+=(Vector const& _vector) - { - gsl_vector_add(m_vector, _vector.m_vector); - return (*this); - } - - Vector & operator-=(Vector const& _vector) - { - gsl_vector_sub(m_vector, _vector.m_vector); - return (*this); - } - - bool is_zero() const - { - return gsl_vector_isnull(m_vector); - } - - bool is_positive() const - { - return gsl_vector_ispos(m_vector); - } - - bool is_negative() const - { - return gsl_vector_isneg(m_vector); - } - - bool is_non_negative() const - { - for ( size_t i = 0; i < size(); ++i ) - { - if ( (*this)[i] < 0 ) return false; - } - return true; - } - - double max() const - { - return gsl_vector_max(m_vector); - } - - double min() const - { - return gsl_vector_min(m_vector); - } - - size_t max_index() const - { - return gsl_vector_max_index(m_vector); - } - - size_t min_index() const - { - return gsl_vector_min_index(m_vector); - } - - friend - void swap(Vector & v1, Vector & v2); - - size_t size() const - { - return m_size; - } - -private: - size_t m_size; - gsl_vector* m_vector; -}; - -void swap(Vector & v1, Vector & v2) -{ - assert( v1.size() == v2.size() ); - std::swap(v1.m_vector, v2.m_vector); -} - -} } // end namespaces - - -#endif /*_NL_VECTOR_H_*/ +#ifndef _NL_VECTOR_H_ +#define _NL_VECTOR_H_ + +#include +#include + +#include + + +namespace Geom { namespace NL { + +class Vector; +void swap(Vector & v1, Vector & v2); + + +class Vector +{ +public: + Vector(size_t n) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + } + + Vector(size_t n, double x) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + gsl_vector_set_all(m_vector, x); + } + + // create a vector with n elements all set to zero + // but the i-th that is set to 1 + Vector(size_t n, size_t i) + : m_size(n) + { + m_vector = gsl_vector_alloc(n); + gsl_vector_set_basis(m_vector, i); + } + + Vector(Vector const& _vector) + : m_size(_vector.size()) + { + m_vector = gsl_vector_alloc(size()); + gsl_vector_memcpy(m_vector, _vector.m_vector); + } + + virtual ~Vector() + { + gsl_vector_free(m_vector); + } + + Vector & operator=(Vector const& _vector) + { + assert( size() == _vector.size() ); + gsl_vector_memcpy(m_vector, _vector.m_vector); + return (*this); + } + + void set_all(double x = 0) + { + gsl_vector_set_all(m_vector, x); + } + + void set_basis(size_t i) + { + gsl_vector_set_basis(m_vector, i); + } + + double const& operator[](size_t i) const + { + return *gsl_vector_const_ptr(m_vector, i); + } + + double & operator[](size_t i) + { + return *gsl_vector_ptr(m_vector, i); + } + + gsl_vector* get_gsl_vector() + { + return m_vector; + } + + void swap_elements(size_t i, size_t j) + { + gsl_vector_swap_elements(m_vector, i, j); + } + + void reverse() + { + gsl_vector_reverse(m_vector); + } + + Vector & scale(double x) + { + gsl_vector_scale(m_vector, x); + return (*this); + } + + Vector & translate(double x) + { + gsl_vector_add_constant(m_vector, x); + return (*this); + } + + Vector & operator+=(Vector const& _vector) + { + gsl_vector_add(m_vector, _vector.m_vector); + return (*this); + } + + Vector & operator-=(Vector const& _vector) + { + gsl_vector_sub(m_vector, _vector.m_vector); + return (*this); + } + + bool is_zero() const + { + return gsl_vector_isnull(m_vector); + } + + bool is_positive() const + { + return gsl_vector_ispos(m_vector); + } + + bool is_negative() const + { + return gsl_vector_isneg(m_vector); + } + + bool is_non_negative() const + { + for ( size_t i = 0; i < size(); ++i ) + { + if ( (*this)[i] < 0 ) return false; + } + return true; + } + + double max() const + { + return gsl_vector_max(m_vector); + } + + double min() const + { + return gsl_vector_min(m_vector); + } + + size_t max_index() const + { + return gsl_vector_max_index(m_vector); + } + + size_t min_index() const + { + return gsl_vector_min_index(m_vector); + } + + friend + void swap(Vector & v1, Vector & v2); + + size_t size() const + { + return m_size; + } + +private: + size_t m_size; + gsl_vector* m_vector; +}; + +void swap(Vector & v1, Vector & v2) +{ + assert( v1.size() == v2.size() ); + std::swap(v1.m_vector, v2.m_vector); +} + +} } // end namespaces + + +#endif /*_NL_VECTOR_H_*/ -- cgit v1.2.3