update to latest 2geom. this adds gsl dependency, doesn't seem to make inskape executable bigger

(bzr r5649)

3 files changed, 480 insertions, 0 deletions

diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
new file mode 100644
index 000000000..2b4963998
--- /dev/null
+++ b/src/2geom/numeric/linear_system.h
@@ -0,0 +1,89 @@
+#ifndef _NL_LINEAR_SYSTEM_H_

+#define _NL_LINEAR_SYSTEM_H_

+

+

+#include <cassert>

+

+#include <gsl/gsl_linalg.h>

+

+#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
new file mode 100644
index 000000000..bdfab75ef
--- /dev/null
+++ b/src/2geom/numeric/matrix.h
@@ -0,0 +1,207 @@
+#ifndef _NL_MATRIX_H_

+#define _NL_MATRIX_H_

+

+

+#include <cassert>

+#include <utility>

+#include <algorithm>

+

+#include <gsl/gsl_matrix.h>

+

+

+

+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<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;

+	}

+	

+	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
new file mode 100644
index 000000000..136b2cc3f
--- /dev/null
+++ b/src/2geom/numeric/vector.h
@@ -0,0 +1,184 @@
+#ifndef _NL_VECTOR_H_

+#define _NL_VECTOR_H_

+

+#include <cassert>

+#include <utility>

+

+#include <gsl/gsl_vector.h>

+

+

+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_*/