Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

8 files changed, 3680 insertions, 0 deletions

diff --git a/include/2geom/numeric/fitting-model.h b/include/2geom/numeric/fitting-model.h
new file mode 100644
index 0000000..0316f57
--- /dev/null
+++ b/include/2geom/numeric/fitting-model.h
@@ -0,0 +1,521 @@
+/*
+ * Fitting Models for Geom Types
+ *
+ * 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_FITTING_MODEL_H_
+#define _NL_FITTING_MODEL_H_
+
+
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
+#include <2geom/bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/polynomial.h>
+#include <2geom/ellipse.h>
+#include <2geom/circle.h>
+#include <2geom/utils.h>
+#include <2geom/conicsec.h>
+
+
+namespace Geom { namespace NL {
+
+/*
+ * A model is an abstraction for an expression dependent from a parameter where
+ * the coefficients of this expression are the unknowns of the fitting problem.
+ * For a ceratain number of parameter values we know the related values
+ * the expression evaluates to: from each parameter value we get a row of
+ * the matrix of the fitting problem, from each expression value we get
+ * the related constant term.
+ * Example: given the model a*x^2 + b*x + c = 0; from x = 1 we get
+ * the equation a + b + c = 0, in this example the constant term is always
+ * the same for each parameter value.
+ *
+ * A model is required to implement 3 methods:
+ *
+ *  - size : returns the number of unknown coefficients that appear in
+ *           the expression of the fitting problem;
+ *  - feed : its input is a parameter value and the related expression value,
+ *           it generates a matrix row and a new entry of the constant vector
+ *           of the fitting problem;
+ *  - instance : it has an input parameter represented by the raw vector
+ *               solution of the fitting problem and an output parameter
+ *               of type InstanceType that return a specific object that is
+ *               generated using the fitting problem solution, in the example
+ *               above the object could be a Poly type.
+ */
+
+/*
+ *   completely unknown models must inherit from this template class;
+ *   example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;
+ *   example: the model A(t) = known_sample_value_at(t) to be solved wrt
+ *       the coefficients of the curve A(t) expressed in S-Basis form;
+ *   parameter type: the type of x and t variable in the examples above;
+ *   value type:     the type of the known sample values (in the first example
+ *                   is constant )
+ *   instance type:  the type of the objects produced by using
+ *                   the fitting raw data solution
+ */
+
+
+
+
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModel
+{
+  public:
+    typedef ParameterType       parameter_type;
+    typedef ValueType           value_type;
+    typedef InstanceType        instance_type;
+
+    static const bool WITH_FIXED_TERMS = false;
+
+    /*
+     * a LinearFittingModel must implement the following methods:
+     *
+     * void feed( VectorView & vector,
+     *            parameter_type const& sample_parameter ) const;
+     *
+     * size_t size() const;
+     *
+     * void instance(instance_type &, raw_type const& raw_data) const;
+     *
+     */
+};
+
+
+/*
+ *   partially known models must inherit from this template class
+ *   example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c
+ */
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModelWithFixedTerms
+{
+  public:
+    typedef ParameterType       parameter_type;
+    typedef ValueType           value_type;
+    typedef InstanceType        instance_type;
+
+    static const bool WITH_FIXED_TERMS = true;
+
+    /*
+     * a LinearFittingModelWithFixedTerms must implement the following methods:
+     *
+     * void feed( VectorView & vector,
+     *            value_type & fixed_term,
+     *            parameter_type const& sample_parameter ) const;
+     *
+     * size_t size() const;
+     *
+     * void instance(instance_type &, raw_type const& raw_data) const;
+     *
+     */
+
+
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// represented in standard power basis
+template< typename InstanceType >
+class LFMPowerBasis
+    : public LinearFittingModel<double, double, InstanceType>
+{
+  public:
+    LFMPowerBasis(size_t degree)
+        : m_size(degree + 1)
+    {
+    }
+
+    void feed( VectorView & coeff, double sample_parameter ) const
+    {
+        coeff[0] = 1;
+        double x_i = 1;
+        for (size_t i = 1; i < coeff.size(); ++i)
+        {
+          x_i *= sample_parameter;
+          coeff[i] = x_i;
+        }
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+  private:
+    size_t m_size;
+};
+
+
+// this model generates Geom::Poly objects
+class LFMPoly
+    : public LFMPowerBasis<Poly>
+{
+  public:
+    LFMPoly(size_t degree)
+        : LFMPowerBasis<Poly>(degree)
+    {
+    }
+
+    void instance(Poly & poly, ConstVectorView const& raw_data) const
+    {
+        poly.clear();
+        poly.resize(size());
+        for (size_t i = 0; i < raw_data.size(); ++i)
+        {
+            poly[i] =  raw_data[i];
+        }
+    }
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// represented in standard power basis with leading term coefficient equal to 1
+template< typename InstanceType >
+class LFMNormalizedPowerBasis
+    : public LinearFittingModelWithFixedTerms<double, double, InstanceType>
+{
+  public:
+    LFMNormalizedPowerBasis(size_t _degree)
+        : m_model( _degree - 1)
+    {
+        assert(_degree > 0);
+    }
+
+
+    void feed( VectorView & coeff,
+               double & known_term,
+               double sample_parameter ) const
+    {
+        m_model.feed(coeff, sample_parameter);
+        known_term = coeff[m_model.size()-1] * sample_parameter;
+    }
+
+    size_t size() const
+    {
+        return m_model.size();
+    }
+
+  private:
+    LFMPowerBasis<InstanceType> m_model;
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of an ellipse curve
+//template< typename InstanceType >
+//class LFMEllipseEquation
+//    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+//{
+//  public:
+//    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+//    {
+//        coeff[0] = p[X] * p[Y];
+//        coeff[1] = p[Y] * p[Y];
+//        coeff[2] = p[X];
+//        coeff[3] = p[Y];
+//        coeff[4] = 1;
+//        fixed_term = p[X] * p[X];
+//    }
+//
+//    size_t size() const
+//    {
+//        return 5;
+//    }
+//};
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of a conic section
+template< typename InstanceType >
+class LFMConicEquation
+    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+{
+  public:
+    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+    {
+        coeff[0] = p[X] * p[Y];
+        coeff[1] = p[Y] * p[Y];
+        coeff[2] = p[X];
+        coeff[3] = p[Y];
+        coeff[4] = 1;
+        fixed_term = p[X] * p[X];
+    }
+
+    size_t size() const
+    {
+        return 5;
+    }
+};
+
+// this model generates Ellipse curves
+class LFMConicSection
+    : public LFMConicEquation<xAx>
+{
+  public:
+    void instance(xAx & c, ConstVectorView const& coeff) const
+    {
+        c.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+    }
+};
+
+// this model generates Ellipse curves
+class LFMEllipse
+    : public LFMConicEquation<Ellipse>
+{
+  public:
+    void instance(Ellipse & e, ConstVectorView const& coeff) const
+    {
+        e.setCoefficients(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+    }
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of a circle curve
+template< typename InstanceType >
+class LFMCircleEquation
+    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+{
+  public:
+    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+    {
+        coeff[0] = p[X];
+        coeff[1] = p[Y];
+        coeff[2] = 1;
+        fixed_term = p[X] * p[X] + p[Y] * p[Y];
+    }
+
+    size_t size() const
+    {
+        return 3;
+    }
+};
+
+
+// this model generates Ellipse curves
+class LFMCircle
+    : public LFMCircleEquation<Circle>
+{
+  public:
+    void instance(Circle & c, ConstVectorView const& coeff) const
+    {
+        c.setCoefficients(1, coeff[0], coeff[1], coeff[2]);
+    }
+};
+
+
+// this model generates SBasis objects
+class LFMSBasis
+    : public LinearFittingModel<double, double, SBasis>
+{
+  public:
+    LFMSBasis( size_t _order )
+        : m_size( 2*(_order+1) ),
+          m_order(_order)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        double u0 = 1-t;
+        double u1 = t;
+        double s = u0 * u1;
+        coeff[0] = u0;
+        coeff[1] = u1;
+        for (size_t i = 2; i < size(); i+=2)
+        {
+            u0 *= s;
+            u1 *= s;
+            coeff[i] = u0;
+            coeff[i+1] = u1;
+        }
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+    void instance(SBasis & sb, ConstVectorView const& raw_data) const
+    {
+        sb.resize(m_order+1);
+        for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)
+        {
+            sb[k][0] = raw_data[i];
+            sb[k][1] = raw_data[i+1];
+        }
+    }
+
+  private:
+    size_t m_size;
+    size_t m_order;
+};
+
+
+// this model generates D2<SBasis> objects
+class LFMD2SBasis
+    : public LinearFittingModel< double, Point, D2<SBasis> >
+{
+  public:
+    LFMD2SBasis( size_t _order )
+        : mosb(_order)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        mosb.feed(coeff, t);
+    }
+
+    size_t size() const
+    {
+        return mosb.size();
+    }
+
+    void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const
+    {
+        mosb.instance(d2sb[X], raw_data.column_const_view(X));
+        mosb.instance(d2sb[Y], raw_data.column_const_view(Y));
+    }
+
+  private:
+    LFMSBasis mosb;
+};
+
+
+// this model generates Bezier objects
+class LFMBezier
+    : public LinearFittingModel<double, double, Bezier>
+{
+  public:
+    LFMBezier( size_t _order )
+        : m_size(_order + 1),
+          m_order(_order)
+    {
+        binomial_coefficients(m_bc, m_order);
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        double s = 1;
+        for (size_t i = 0; i < size(); ++i)
+        {
+            coeff[i] = s * m_bc[i];
+            s *= t;
+        }
+        double u = 1-t;
+        s = 1;
+        for (size_t i = size()-1; i > 0; --i)
+        {
+            coeff[i] *= s;
+            s *= u;
+        }
+        coeff[0] *= s;
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+    void instance(Bezier & b, ConstVectorView const& raw_data) const
+    {
+        assert(b.size() == raw_data.size());
+        for (unsigned int i = 0; i < raw_data.size(); ++i)
+        {
+            b[i] = raw_data[i];
+        }
+    }
+
+  private:
+    size_t m_size;
+    size_t m_order;
+    std::vector<size_t> m_bc;
+};
+
+
+// this model generates Bezier curves
+template <unsigned degree>
+class LFMBezierCurveN
+    : public LinearFittingModel< double, Point, BezierCurveN<degree> >
+{
+  public:
+    LFMBezierCurveN()
+        : mob(degree+1)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        mob.feed(coeff, t);
+    }
+
+    size_t size() const
+    {
+        return mob.size();
+    }
+
+    void instance(BezierCurveN<degree> & bc, ConstMatrixView const& raw_data) const
+    {
+        Bezier bx(degree);
+        Bezier by(degree);
+        mob.instance(bx, raw_data.column_const_view(X));
+        mob.instance(by, raw_data.column_const_view(Y));
+        bc = BezierCurveN<degree>(bx, by);
+    }
+
+  private:
+    LFMBezier mob;
+};
+
+}  // end namespace NL
+}  // end namespace Geom
+
+
+#endif // _NL_FITTING_MODEL_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 :
diff --git a/include/2geom/numeric/fitting-tool.h b/include/2geom/numeric/fitting-tool.h
new file mode 100644
index 0000000..78d66ca
--- /dev/null
+++ b/include/2geom/numeric/fitting-tool.h
@@ -0,0 +1,562 @@
+/*
+ * Fitting Tools
+ *
+ * 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_FITTING_TOOL_H_
+#define _NL_FITTING_TOOL_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+
+#include <2geom/point.h>
+
+#include <vector>
+
+
+/*
+ *  The least_square_fitter class represents a tool for solving a fitting
+ *  problem with respect to a given model that represents an expression
+ *  dependent from a parameter where the coefficients of this expression
+ *  are the unknowns of the fitting problem.
+ *  The minimizing solution is found by computing the pseudo-inverse
+ *  of the problem matrix
+ */
+
+
+namespace Geom { namespace NL {
+
+namespace detail {
+
+template< typename ModelT>
+class lsf_base
+{
+  public:
+    typedef ModelT                                   model_type;
+    typedef typename model_type::parameter_type      parameter_type;
+    typedef typename model_type::value_type          value_type;
+
+    lsf_base( model_type const& _model, size_t forecasted_samples )
+        : m_model(_model),
+          m_total_samples(0),
+          m_matrix(forecasted_samples, m_model.size()),
+          m_psdinv_matrix(NULL)
+    {
+    }
+
+    // compute pseudo inverse
+    void update()
+    {
+        if (total_samples() == 0) return;
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+        MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
+        m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
+        assert(m_psdinv_matrix != NULL);
+    }
+
+    size_t total_samples() const
+    {
+        return m_total_samples;
+    }
+
+    bool is_full() const
+    {
+        return (total_samples() == m_matrix.rows());
+    }
+
+    void clear()
+    {
+        m_total_samples = 0;
+    }
+
+    virtual
+    ~lsf_base()
+    {
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+    }
+
+  protected:
+    const model_type &      m_model;
+    size_t                  m_total_samples;
+    Matrix                  m_matrix;
+    Matrix*                 m_psdinv_matrix;
+
+};  // end class lsf_base
+
+
+
+
+template< typename ModelT, typename ValueType = typename ModelT::value_type>
+class lsf_solution
+{
+};
+
+// a fitting process on samples with value of type double
+// produces a solution of type Vector
+template< typename ModelT>
+class lsf_solution<ModelT, double>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Vector                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, double>( model_type const& _model,
+                                  size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size())
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        ConstVectorView sv(sample_values);
+        m_solution = (*m_psdinv_matrix) * sv;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can be performed between two result invocations
+    template< typename VectorT >
+    solution_type& result( VectorT const& old_sample_values,
+                           VectorT const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Vector diff(total_samples());
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            diff[i] = new_sample_values[i] - old_sample_values[i];
+        }
+        Vector column(m_model.size());
+        Vector delta(m_model.size(), 0.0);
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            if (diff[i] != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff[i]);
+                delta += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, double>
+
+
+// a fitting process on samples with value of type Point
+// produces a solution of type Matrix (with 2 columns)
+template< typename ModelT>
+class lsf_solution<ModelT, Point>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Matrix                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, Point>( model_type const& _model,
+                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size(), 2)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        Matrix svm(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            svm(i, X) = sample_values[i][X];
+            svm(i, Y) = sample_values[i][Y];
+        }
+        m_solution = (*m_psdinv_matrix) * svm;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can to be performed between two result invocations
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Matrix diff(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            diff(i, X) =  new_sample_values[i][X] - old_sample_values[i][X];
+            diff(i, Y) =  new_sample_values[i][Y] - old_sample_values[i][Y];
+        }
+        Vector column(m_model.size());
+        Matrix delta(m_model.size(), 2, 0.0);
+        VectorView deltax = delta.column_view(X);
+        VectorView deltay = delta.column_view(Y);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            if (diff(i, X) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, X));
+                deltax += column;
+            }
+            if (diff(i, Y) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, Y));
+                deltay += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, Point>
+
+
+
+
+template< typename ModelT,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class lsf_with_fixed_terms
+{
+};
+
+
+// fitting tool for completely unknown models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, false>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+      lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
+                                           size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        ++m_total_samples;
+    }
+
+}; // end class lsf_with_fixed_terms<ModelT, false>
+
+
+// fitting tool for partially known models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, true>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+    lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
+                                        size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_vector(forecasted_samples),
+          m_vector_view(NULL)
+    {
+    }
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, m_vector[total_samples()], sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        m_vector[total_samples()] = m_vector[sample_index];
+        ++m_total_samples;
+    }
+
+    void update()
+    {
+        base_type::update();
+        if (total_samples() == 0) return;
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+        m_vector_view = new VectorView(m_vector, base_type::total_samples());
+        assert(m_vector_view != NULL);
+    }
+
+    
+    ~lsf_with_fixed_terms<model_type, true>() override
+    {
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+    }
+
+  protected:
+    Vector          m_vector;
+    VectorView*     m_vector_view;
+
+}; // end class lsf_with_fixed_terms<ModelT, true>
+
+
+} // end namespace detail
+
+
+
+
+template< typename ModelT,
+          typename ValueType = typename ModelT::value_type,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class least_squeares_fitter
+{
+};
+
+
+template< typename ModelT, typename ValueType >
+class least_squeares_fitter<ModelT, ValueType, false>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+  public:
+    least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,
+                                                     size_t forecasted_samples )
+          : base_type(_model, forecasted_samples)
+  {
+  }
+}; // end class least_squeares_fitter<ModelT, ValueType, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, double, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    //using base_type::result; // VSC legacy support
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        return base_type::result(old_sample_values, new_sample_values);
+    }
+ 
+    solution_type& result()
+    {
+        return base_type::result();
+    }
+
+  public:
+    least_squeares_fitter<ModelT, double, true>( model_type const& _model,
+                                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        Vector sv(sample_values.size());
+        for (size_t i = 0; i < sv.size(); ++i)
+            sv[i] = sample_values[i] - (*m_vector_view)[i];
+        return base_type::result(sv);
+    }
+
+}; // end class least_squeares_fitter<ModelT, double, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, Point, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    //using base_type::result; // VCS legacy support
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        return base_type::result(old_sample_values, new_sample_values);
+    }
+ 
+    solution_type& result()
+    {
+        return base_type::result();
+    }
+
+
+  public:
+    least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
+                                                size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        NL::Matrix sv(sample_values.size(), 2);
+        for (size_t i = 0; i < sample_values.size(); ++i)
+        {
+            sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
+            sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
+        }
+        return base_type::result(sv);
+    }
+
+};  // end class least_squeares_fitter<ModelT, Point, true>
+
+
+}  // end namespace NL
+}  // end namespace Geom
+
+
+
+#endif  // _NL_FITTING_TOOL_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 :
diff --git a/include/2geom/numeric/linear_system.h b/include/2geom/numeric/linear_system.h
new file mode 100644
index 0000000..f793e20
--- /dev/null
+++ b/include/2geom/numeric/linear_system.h
@@ -0,0 +1,138 @@
+/*
+ * LinearSystem class wraps some gsl routines for solving linear systems
+ *
+ * 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_LINEAR_SYSTEM_H_
+#define _NL_LINEAR_SYSTEM_H_
+
+
+#include <cassert>
+
+#include <gsl/gsl_linalg.h>
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/vector.h>
+
+
+namespace Geom { namespace NL {
+
+
+class LinearSystem
+{
+public:
+	LinearSystem(MatrixView & _matrix, VectorView & _vector)
+		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
+	{
+	}
+	
+	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();			  
+	}
+	
+	MatrixView & matrix()
+	{
+		return m_matrix;
+	}
+	
+	VectorView & vector()
+	{
+		return m_vector;
+	}
+	
+	const Vector & solution() const
+	{
+		return m_solution;
+	}
+	
+private:
+	MatrixView m_matrix;
+	VectorView m_vector;
+	Vector m_solution;
+};
+
+
+} } // end namespaces
+
+
+#endif /*_NL_LINEAR_SYSTEM_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 :
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 :
diff --git a/include/2geom/numeric/symmetric-matrix-fs-operation.h b/include/2geom/numeric/symmetric-matrix-fs-operation.h
new file mode 100644
index 0000000..c5aaa72
--- /dev/null
+++ b/include/2geom/numeric/symmetric-matrix-fs-operation.h
@@ -0,0 +1,102 @@
+/*
+ * SymmetricMatrix basic operation
+ *
+ * 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_OPERATION_H_
+#define _NL_SYMMETRIC_MATRIX_FS_OPERATION_H_
+
+
+#include <2geom/numeric/symmetric-matrix-fs.h>
+#include <2geom/numeric/symmetric-matrix-fs-trace.h>
+
+
+
+
+namespace Geom { namespace NL {
+
+template <size_t N>
+SymmetricMatrix<N> adj(const ConstBaseSymmetricMatrix<N> & S);
+
+template <>
+inline
+SymmetricMatrix<2> adj(const ConstBaseSymmetricMatrix<2> & S)
+{
+    SymmetricMatrix<2> result;
+    result.get<0,0>() = S.get<1,1>();
+    result.get<1,0>() = -S.get<1,0>();
+    result.get<1,1>() = S.get<0,0>();
+    return result;
+}
+
+template <>
+inline
+SymmetricMatrix<3> adj(const ConstBaseSymmetricMatrix<3> & S)
+{
+    SymmetricMatrix<3> result;
+
+    result.get<0,0>() = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+    result.get<1,0>() = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+    result.get<1,1>() = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+    result.get<2,0>() = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+    result.get<2,1>() = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+    result.get<2,2>() = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    return result;
+}
+
+template <size_t N>
+inline
+SymmetricMatrix<N> inverse(const ConstBaseSymmetricMatrix<N> & S)
+{
+    SymmetricMatrix<N> result = adj(S);
+    double d = det(S);
+    assert (d != 0);
+    result.scale (1/d);
+    return result;
+}
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+#endif // _NL_SYMMETRIC_MATRIX_FS_OPERATION_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 :
diff --git a/include/2geom/numeric/symmetric-matrix-fs-trace.h b/include/2geom/numeric/symmetric-matrix-fs-trace.h
new file mode 100644
index 0000000..0e7a28c
--- /dev/null
+++ b/include/2geom/numeric/symmetric-matrix-fs-trace.h
@@ -0,0 +1,427 @@
+/*
+ * SymmetricMatrix trace
+ *
+ * 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_TRACE_H_
+#define _NL_TRACE_H_
+
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/symmetric-matrix-fs.h>
+
+
+
+
+
+namespace Geom { namespace NL {
+
+
+namespace detail
+{
+
+/*
+ *  helper routines
+ */
+
+inline
+int sgn_prod (int x, int y)
+{
+    if (x == 0 || y == 0)  return 0;
+    if (x == y) return 1;
+    return -1;
+}
+
+inline
+bool abs_less (double x, double y)
+{
+    return (std::fabs(x) < std::fabs(y));
+}
+
+
+/*
+ * trace K-th of symmetric matrix S of order N
+ */
+template <size_t K, size_t N>
+struct trace
+{
+    static double evaluate(const ConstBaseSymmetricMatrix<N> &S);
+};
+
+template <size_t N>
+struct trace<1,N>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        double t = 0;
+        for (size_t i = 0; i < N; ++i)
+        {
+            t += S(i,i);
+        }
+        return t;
+    }
+};
+
+template <size_t N>
+struct trace<N,N>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        Matrix M(S);
+        return det(M);
+    }
+};
+
+/*
+ *  trace for symmetric matrix of order 2
+ */
+template <>
+struct trace<1,2>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        return (S.get<0,0>() + S.get<1,1>());
+    }
+};
+
+template <>
+struct trace<2,2>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        return (S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>());
+    }
+};
+
+
+/*
+ *  trace for symmetric matrix of order 3
+ */
+template <>
+struct trace<1,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        return (S.get<0,0>() + S.get<1,1>() + S.get<2,2>());
+    }
+};
+
+template <>
+struct trace<2,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+        double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+        double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+        return (a00 + a11 + a22);
+    }
+};
+
+template <>
+struct trace<3,3>
+{
+    static
+    double evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double d = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+        d += (2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>());
+        d -= (S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+        d -= (S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+        d -= (S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+        return d;
+    }
+};
+
+
+/*
+ *  sign of trace K-th
+ */
+template <size_t K, size_t N>
+struct trace_sgn
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<N> & S)
+    {
+        double d = trace<K, N>::evaluate(S);
+        return sgn(d);
+    }
+};
+
+
+/*
+ *  sign of trace for symmetric matrix of order 2
+ */
+template <>
+struct trace_sgn<2,2>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<2> & S)
+    {
+        double m00 = S.get<0,0>();
+        double m10 = S.get<1,0>();
+        double m11 = S.get<1,1>();
+
+        int sm00 = sgn (m00);
+        int sm10 = sgn (m10);
+        int sm11 = sgn (m11);
+
+        if (sm10 == 0)
+        {
+            return sgn_prod (sm00, sm11);
+        }
+        else
+        {
+            int sm00m11 = sgn_prod (sm00, sm11);
+            if (sm00m11 == 1)
+            {
+                int e00, e10, e11;
+                double f00 = std::frexp (m00, &e00);
+                double f10 = std::frexp (m10, &e10);
+                double f11 = std::frexp (m11, &e11);
+
+                int e0011 = e00 + e11;
+                int e1010 = e10 << 1;
+                int ed = e0011 - e1010;
+
+                if (ed > 1)
+                {
+                    return 1;
+                }
+                else if (ed < -1)
+                {
+                    return -1;
+                }
+                else
+                {
+                    double d = std::ldexp (f00 * f11, ed) - f10 * f10;
+                    //std::cout << "trace_sgn<2,2>: det = " << d << std::endl;
+                    double eps = std::ldexp (1, -50);
+                    if (std::fabs(d) < eps)  return 0;
+                    return sgn (d);
+                }
+            }
+            return -1;
+        }
+    }
+};
+
+
+/*
+ *  sign of trace for symmetric matrix of order 3
+ */
+template <>
+struct trace_sgn<2,3>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+        double eps = std::ldexp (1, -50);
+        double t[6];
+
+        t[0] = S.get<1,1>() * S.get<2,2>();
+        t[1] = - S.get<1,2>() * S.get<2,1>();
+        t[2] = S.get<0,0>() * S.get<2,2>();
+        t[3] = - S.get<0,2>() * S.get<2,0>();
+        t[4] = S.get<0,0>() * S.get<1,1>();
+        t[5] = - S.get<0,1>() * S.get<1,0>();
+
+
+        double* maxp = std::max_element (t, t+6, abs_less);
+        int em;
+        std::frexp(*maxp, &em);
+        double d = 0;
+        for (double i : t)
+        {
+            d += i;
+        }
+        double r = std::fabs (std::ldexp (d, -em));  // relative error
+        //std::cout << "trace_sgn<2,3>: d = " << d << std::endl;
+        //std::cout << "trace_sgn<2,3>: r = " << r << std::endl;
+        if (r < eps) return 0;
+        if (d > 0) return 1;
+        return -1;
+    }
+};
+
+template <>
+struct trace_sgn<3,3>
+{
+    static
+    int evaluate (const ConstBaseSymmetricMatrix<3> & S)
+    {
+
+        double eps = std::ldexp (1, -48);
+        double t[5];
+
+        t[0] = S.get<0,0>() * S.get<1,1>() * S.get<2,2>();
+        t[1] = 2 * S.get<1,0>() * S.get<2,0>() * S.get<2,1>();
+        t[2] = -(S.get<0,0>() * S.get<2,1>() * S.get<2,1>());
+        t[3] = -(S.get<1,1>() * S.get<2,0>() * S.get<2,0>());
+        t[4] = -(S.get<2,2>() * S.get<1,0>() * S.get<1,0>());
+
+        double* maxp = std::max_element (t, t+5, abs_less);
+        int em;
+        std::frexp(*maxp, &em);
+        double d = 0;
+        for (double i : t)
+        {
+            d += i;
+        }
+        //std::cout << "trace_sgn<3,3>: d = " << d << std::endl;
+        double r = std::fabs (std::ldexp (d, -em));  // relative error
+        //std::cout << "trace_sgn<3,3>: r = " << r << std::endl;
+
+        if (r < eps)  return 0;
+        if (d > 0) return 1;
+        return -1;
+    }
+}; // end struct trace_sgn<3,3>
+
+} // end namespace detail
+
+
+template <size_t K, size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double trace (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+double det (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace<N, N>::evaluate(_matrix);
+}
+
+
+template <size_t K, size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<K, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int trace_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<1, N>::evaluate(_matrix);
+}
+
+template <size_t N>
+inline
+int det_sgn (const ConstBaseSymmetricMatrix<N> & _matrix)
+{
+    return detail::trace_sgn<N, N>::evaluate(_matrix);
+}
+
+/*
+template <size_t N>
+inline
+size_t rank (const ConstBaseSymmetricMatrix<N> & S)
+{
+    THROW_NOTIMPLEMENTED();
+    return 0;
+}
+
+template <>
+inline
+size_t rank<2> (const ConstBaseSymmetricMatrix<2> & S)
+{
+    if (S.is_zero())  return 0;
+    double d = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    if (d != 0)  return 2;
+    return 1;
+}
+
+template <>
+inline
+size_t rank<3> (const ConstBaseSymmetricMatrix<3> & S)
+{
+    if (S.is_zero())  return 0;
+
+    double a20 = S.get<0,1>() * S.get<1,2>() - S.get<0,2>() * S.get<1,1>();
+    double a21 = S.get<0,2>() * S.get<1,0>() - S.get<0,0>() * S.get<1,2>();
+    double a22 = S.get<0,0>() * S.get<1,1>() - S.get<0,1>() * S.get<1,0>();
+    double d = a20 * S.get<2,0>() + a21 * S.get<2,1>() + a22 * S.get<2,2>();
+
+    if (d != 0)  return 3;
+
+    if (a20 != 0 || a21 != 0 || a22 != 0)  return 2;
+
+    double a00 = S.get<1,1>() * S.get<2,2>() - S.get<1,2>() * S.get<2,1>();
+    if (a00 != 0)  return 2;
+
+    double a10 = S.get<0,2>() * S.get<2,1>() - S.get<0,1>() * S.get<2,2>();
+    if (a10 != 0)  return 2;
+
+    double a11 = S.get<0,0>() * S.get<2,2>() - S.get<0,2>() * S.get<2,0>();
+    if (a11 != 0)  return 2;
+
+    return 1;
+}
+*/
+
+} /* end namespace NL*/ } /* end namespace Geom*/
+
+
+
+
+#endif // _NL_TRACE_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 :
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 :
diff --git a/include/2geom/numeric/vector.h b/include/2geom/numeric/vector.h
new file mode 100644
index 0000000..b66115b
--- /dev/null
+++ b/include/2geom/numeric/vector.h
@@ -0,0 +1,594 @@
+/*
+ * Vector, VectorView, ConstVectorView classes wrap the gsl vector 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_VECTOR_H_
+#define _NL_VECTOR_H_
+
+#include <cassert>
+#include <algorithm> // for std::swap
+#include <vector>
+#include <sstream>
+#include <string>
+
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseVectorImpl
+{
+  public:
+	double const& operator[](size_t i) const
+	{
+		return *gsl_vector_const_ptr(m_vector, i);
+	}
+
+	const gsl_vector* get_gsl_vector() const
+	{
+		return m_vector;
+	}
+	bool is_zero() const
+	{
+		return gsl_vector_isnull(m_vector);
+	}
+
+	bool is_positive() const
+	{
+	    for ( size_t i = 0; i < size(); ++i )
+	    {
+	        if ( (*this)[i] <= 0 ) return false;
+	    }
+		return true;
+	}
+
+	bool is_negative() const
+	{
+        for ( size_t i = 0; i < size(); ++i )
+        {
+            if ( (*this)[i] >= 0 ) return false;
+        }
+        return true;
+	}
+
+	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);
+	}
+
+	size_t size() const
+	{
+		return m_size;
+	}
+
+	std::string str() const;
+
+	virtual ~BaseVectorImpl()
+	{
+	}
+
+  protected:
+	size_t m_size;
+	gsl_vector* m_vector;
+
+};  // end class BaseVectorImpl
+
+
+inline
+bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+	if (v1.size() != v2.size())	return false;
+
+	for (size_t i = 0; i < v1.size(); ++i)
+	{
+		if (v1[i] != v2[i])  return false;
+	}
+	return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)
+{
+	if (_vector.size() == 0 ) return os;
+	os << "[" << _vector[0];
+	for (unsigned int i = 1; i < _vector.size(); ++i)
+	{
+		os << ", " << _vector[i];
+	}
+	os << "]";
+	return os;
+}
+
+inline
+std::string BaseVectorImpl::str() const
+{
+	std::ostringstream oss;
+	oss << (*this);
+	return oss.str();
+}
+
+inline
+double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+    double result;
+    gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);
+    return result;
+}
+
+
+class VectorImpl : public BaseVectorImpl
+{
+  public:
+	typedef BaseVectorImpl base_type;
+
+  public:
+	void set_all(double x)
+	{
+		gsl_vector_set_all(m_vector, x);
+	}
+
+	void set_basis(size_t i)
+	{
+		gsl_vector_set_basis(m_vector, i);
+	}
+
+	using base_type::operator[];
+
+	double & operator[](size_t i)
+	{
+		return *gsl_vector_ptr(m_vector, i);
+	}
+
+	using base_type::get_gsl_vector;
+
+	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);
+	}
+
+	VectorImpl & scale(double x)
+	{
+		gsl_vector_scale(m_vector, x);
+		return (*this);
+	}
+
+	VectorImpl & translate(double x)
+	{
+		gsl_vector_add_constant(m_vector, x);
+		return (*this);
+	}
+
+	VectorImpl & operator+=(base_type const& _vector)
+	{
+		gsl_vector_add(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorImpl & operator-=(base_type const& _vector)
+	{
+		gsl_vector_sub(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+};  // end class VectorImpl
+
+}  // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+class Vector : public detail::VectorImpl
+{
+  public:
+	typedef detail::VectorImpl base_type;
+
+  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)
+        : base_type()
+	{
+		m_size = _vector.size();
+		m_vector = gsl_vector_alloc(size());
+		gsl_vector_memcpy(m_vector, _vector.m_vector);
+	}
+
+	explicit
+	Vector(base_type::base_type const& _vector)
+	{
+		m_size = _vector.size();
+		m_vector = gsl_vector_alloc(size());
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+	}
+
+	~Vector() override
+	{
+		gsl_vector_free(m_vector);
+	}
+
+
+	Vector & operator=(Vector const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.m_vector);
+		return (*this);
+	}
+
+	Vector & operator=(base_type::base_type const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	Vector & scale(double x)
+	{
+		return static_cast<Vector&>( base_type::scale(x) );
+	}
+
+	Vector & translate(double x)
+	{
+		return static_cast<Vector&>( base_type::translate(x) );
+	}
+
+	Vector & operator+=(base_type::base_type const& _vector)
+	{
+		return static_cast<Vector&>( base_type::operator+=(_vector) );
+	}
+
+	Vector & operator-=(base_type::base_type const& _vector)
+	{
+		return static_cast<Vector&>( base_type::operator-=(_vector) );
+	}
+
+	friend
+	void swap(Vector & v1, Vector & v2);
+	friend
+	void swap_any(Vector & v1, Vector & v2);
+
+}; // end class Vector
+
+
+// warning! these operations invalidate any view of the passed vector objects
+inline
+void swap(Vector & v1, Vector & v2)
+{
+    assert(v1.size() == v2.size());
+    using std::swap;
+    swap(v1.m_vector, v2.m_vector);
+}
+
+inline
+void swap_any(Vector & v1, Vector & v2)
+{
+	using std::swap;
+    swap(v1.m_vector, v2.m_vector);
+    swap(v1.m_size, v2.m_size);
+}
+
+
+class ConstVectorView : public detail::BaseVectorImpl
+{
+  public:
+	typedef detail::BaseVectorImpl base_type;
+
+  public:
+	ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)
+		: m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )
+	{
+		m_size = n;
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+	ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)
+		: m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )
+	{
+		m_size = n;
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+    ConstVectorView(const double* _vector, size_t n, size_t offset = 0)
+        : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )
+    {
+        m_size = n;
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+    }
+
+    ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)
+        : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )
+    {
+        m_size = n;
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+    }
+
+	explicit
+	ConstVectorView(gsl_vector_const_view  _gsl_vector_view)
+		: m_vector_view(_gsl_vector_view)
+	{
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+		m_size = m_vector->size;
+	}
+
+    explicit
+    ConstVectorView(const std::vector<double>&  _vector)
+        : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )
+    {
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+        m_size = _vector.size();
+    }
+
+	ConstVectorView(const ConstVectorView & _vector)
+		: base_type(),
+		  m_vector_view(_vector.m_vector_view)
+	{
+		m_size = _vector.size();
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+	ConstVectorView(const base_type & _vector)
+		: m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))
+	{
+		m_size = _vector.size();
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+  private:
+	gsl_vector_const_view m_vector_view;
+
+}; // end class  ConstVectorView
+
+
+
+
+class VectorView : public detail::VectorImpl
+{
+  public:
+	typedef detail::VectorImpl base_type;
+
+  public:
+	VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)
+	{
+		m_size = n;
+		if (stride == 1)
+		{
+			m_vector_view
+				= gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);
+			m_vector = &(m_vector_view.vector);
+		}
+		else
+		{
+			m_vector_view
+				= gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);
+			m_vector = &(m_vector_view.vector);
+		}
+	}
+
+    VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)
+    {
+        m_size = n;
+        if (stride == 1)
+        {
+            m_vector_view
+                = gsl_vector_view_array(_vector + offset, n);
+            m_vector = &(m_vector_view.vector);
+        }
+        else
+        {
+            m_vector_view
+                = gsl_vector_view_array_with_stride(_vector + offset, stride, n);
+            m_vector = &(m_vector_view.vector);
+        }
+
+    }
+
+	VectorView(const VectorView & _vector)
+        : base_type()
+	{
+		m_size = _vector.size();
+		m_vector_view = _vector.m_vector_view;
+		m_vector = &(m_vector_view.vector);
+	}
+
+	VectorView(Vector & _vector)
+	{
+		m_size = _vector.size();
+		m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());
+		m_vector = &(m_vector_view.vector);
+	}
+
+	explicit
+	VectorView(gsl_vector_view _gsl_vector_view)
+		: m_vector_view(_gsl_vector_view)
+	{
+		m_vector = &(m_vector_view.vector);
+		m_size = m_vector->size;
+	}
+
+	explicit
+	VectorView(std::vector<double> & _vector)
+	{
+	    m_size = _vector.size();
+	    m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());
+	    m_vector = &(m_vector_view.vector);
+	}
+
+	VectorView & operator=(VectorView const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorView & operator=(base_type::base_type const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorView & scale(double x)
+	{
+		return static_cast<VectorView&>( base_type::scale(x) );
+	}
+
+	VectorView & translate(double x)
+	{
+		return static_cast<VectorView&>( base_type::translate(x) );
+	}
+
+	VectorView & operator+=(base_type::base_type const& _vector)
+	{
+		return static_cast<VectorView&>( base_type::operator+=(_vector) );
+	}
+
+	VectorView & operator-=(base_type::base_type const& _vector)
+	{
+		return static_cast<VectorView&>( base_type::operator-=(_vector) );
+	}
+
+	friend
+	void swap_view(VectorView & v1, VectorView & v2);
+
+  private:
+	gsl_vector_view m_vector_view;
+
+}; // end class VectorView
+
+
+inline
+void swap_view(VectorView & v1, VectorView & v2)
+{
+	assert( v1.size() == v2.size() );
+	using std::swap;
+	swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
+}
+
+inline
+const VectorView & const_vector_view_cast (const ConstVectorView & view)
+{
+    const detail::BaseVectorImpl & bvi
+        = static_cast<const detail::BaseVectorImpl &>(view);
+    const VectorView & vv = reinterpret_cast<const VectorView &>(bvi);
+    return vv;
+}
+
+inline
+VectorView & const_vector_view_cast (ConstVectorView & view)
+{
+    detail::BaseVectorImpl & bvi
+        = static_cast<detail::BaseVectorImpl &>(view);
+    VectorView & vv = reinterpret_cast<VectorView &>(bvi);
+    return vv;
+}
+
+
+} } // end namespaces
+
+
+#endif /*_NL_VECTOR_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 :