2
* LinearSystem class wraps some gsl routines for solving linear systems
5
* Marco Cecchetti <mrcekets at gmail.com>
7
* Copyright 2008 authors
9
* This library is free software; you can redistribute it and/or
10
* modify it either under the terms of the GNU Lesser General Public
11
* License version 2.1 as published by the Free Software Foundation
12
* (the "LGPL") or, at your option, under the terms of the Mozilla
13
* Public License Version 1.1 (the "MPL"). If you do not alter this
14
* notice, a recipient may use your version of this file under either
15
* the MPL or the LGPL.
17
* You should have received a copy of the LGPL along with this library
18
* in the file COPYING-LGPL-2.1; if not, write to the Free Software
19
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
* You should have received a copy of the MPL along with this library
21
* in the file COPYING-MPL-1.1
23
* The contents of this file are subject to the Mozilla Public License
24
* Version 1.1 (the "License"); you may not use this file except in
25
* compliance with the License. You may obtain a copy of the License at
26
* http://www.mozilla.org/MPL/
28
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
29
* OF ANY KIND, either express or implied. See the LGPL or the MPL for
30
* the specific language governing rights and limitations.
34
#ifndef _NL_LINEAR_SYSTEM_H_
35
#define _NL_LINEAR_SYSTEM_H_
40
#include <gsl/gsl_linalg.h>
42
#include <2geom/numeric/matrix.h>
43
#include <2geom/numeric/vector.h>
46
namespace Geom { namespace NL {
52
LinearSystem(MatrixView & _matrix, VectorView & _vector)
53
: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
57
LinearSystem(Matrix & _matrix, Vector & _vector)
58
: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
62
const Vector & LU_solve()
64
assert( matrix().rows() == matrix().columns()
65
&& matrix().rows() == vector().size() );
67
gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
68
gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
69
gsl_linalg_LU_solve( matrix().get_gsl_matrix(),
71
vector().get_gsl_vector(),
72
m_solution.get_gsl_vector()
74
gsl_permutation_free(p);
78
const Vector & SV_solve()
80
assert( matrix().rows() >= matrix().columns()
81
&& matrix().rows() == vector().size() );
83
gsl_matrix* U = matrix().get_gsl_matrix();
84
gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
85
gsl_vector* S = gsl_vector_alloc(matrix().columns());
86
gsl_vector* work = gsl_vector_alloc(matrix().columns());
88
gsl_linalg_SV_decomp( U, V, S, work );
90
gsl_vector* b = vector().get_gsl_vector();
91
gsl_vector* x = m_solution.get_gsl_vector();
93
gsl_linalg_SV_solve( U, V, S, b, x);
97
gsl_vector_free(work);
102
MatrixView & matrix()
107
VectorView & vector()
112
const Vector & solution() const
124
} } // end namespaces
127
#endif /*_NL_LINEAR_SYSTEM_H_*/
132
c-file-style:"stroustrup"
133
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
138
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :